How did Gauss sum Eisenstein series? Entry 61 in Gauss's diary states (this is a translation from latin):

From the integer powers of $$\int_0^1 \frac{dx}{\sqrt{1-x^4}}$$
depends $$\sum_{m,n}(\frac{m^4 - 6m^2n^2 +n^4}{(m^2+n^2)^4})^k$$.

The expression inside the summation is: $\frac{m^4 - 6m^2n^2 +n^4}{(m^2+n^2)^4} = \frac{1}{(m+n\sqrt{-1})^4}+\frac{1}{(m-n\sqrt{-1})^4} $, so this entry apparently connects Eisenstein series with integral powers of half the "lemniscate constant" $\varpi = 2\int_0^1 \frac{dx}{\sqrt{1-x^4}}$. Note also of the analogy with "Basel problem"; here the summation of the fourth powers of inverse "Gaussian integers" rationally depends on the fourth power not of $\pi = 2\int_0^1 \frac{dx}{\sqrt{1-x^2}}$ but rather on the fourth power of its lemniscatic analogue.
Felix klein and Ludwig Schlesinger, the mathematicians who commented on this diary entry, interpreted Gauss's entry as the statement that the sum of Eisenstein series of weight $4k$ is a rational number times $\varpi^{4k}$. I also found additional verification of the conjecture that Gauss apparently summed up Eisenstein series in several other sources (mentioned below).
Schlesinger and Klein offer a reconstruction of the way Gauss arrived at the statement in entry 61 in p. 516 of volume 10.1 of his works. They say he began with the infinite product expansion of $\mathbb{sinlemn}(z) = \frac{M(z)}{N(z)}$, and then arrived somehow at an infinite series for the logarithm of the numerator $M(z)$ and the denominator $N(z)$. The numerator $M(z)$ has zeros at Gaussian integers multiples of $\varpi$ (that is: $M((m+ni)\varpi) = 0$, where $m,n$ are integers), while $N(z)$ has zeros at Gaussian half-integers multiples of $\varpi$, and therefore:
$$M(z) = z\cdot \prod_{(m,n)\ne (0,0)}(1-\frac{z}{(m+ni)\varpi})$$
$$N(z) = \prod (1-\frac{z}{((m+\frac{1}{2})+(n+\frac{1}{2})i)\varpi})$$
Since the logarithm of an infinite product equals an infinite sum of logarithms, one gets:
$$(*)\mathbb{log}M(z) = \mathbb{log}(z)+\sum_{(m,n)\ne (0,0)}\mathbb{log}(1-\frac{z}{(m+ni)\varpi})$$
$$(**)\mathbb{log}N(z) = \sum \mathbb{log}(1-\frac{z}{((m+\frac{1}{2})+(n+\frac{1}{2})i)\varpi})$$.
In addition, on p.159 of volume 10.1 in Gauss's werke, Gauss somehow gives the following infinite series:
$$(***)\mathbb{log}M(z) = \mathbb{log}(z) - \frac{1}{60}z^4 - \frac{1}{4200}z^8 -\frac{1}{321750}z^{12} - ...$$
$$(****)\mathbb{log}N(z) =\frac{1}{12}z^4 - \frac{1}{280}z^8 +\frac{1}{4950}z^{12} - ...$$.
Comparing the coefficients of the series for $\mathbb{log}M(z)$ with the coefficients that result from (*) apparently enabled Gauss to sum particular cases of Eisenstein series, and to come up with the statement in his diary entry.
Now, I have two misunderstandings concerning this conjecture:

*

*The sum $$\sum_{m,n} (\frac{1}{(m+n\sqrt{-1})^4}+\frac{1}{(m-n\sqrt{-1})^4})^k $$ is not exactly   the canonical form of an Eisenstein series because of the second fraction inside the basis of the $k$ power. Therefore, how to show the equivalence of the two sums?

*I don't understand how the evaluation of the logarithm of the numerator $\mathbb{log}M(z)$ amounts to summation of Eisenstein series. I simply need to see the algebra written more clearly. So can someone explain to me how this reasoning amounts to summing Eisenstein series?

Update (13 February, 2022)
An additional point of interest that I noted now is that Gauss's series for $\mathbb{log}N(z)$ corresponds to a kind of "generalized Eisenstein series" in which the lattice is shifted by $(\frac{1}{2}+\frac{1}{2}i)\varpi$. Since  such shifted lattice cannot be generated by some action of the modular group on the lattice of Gaussian integers (if it was, one could use the modularity of the Eisenstein series and deduce the series for $\mathbb{log}N(z)$ from that of $\mathbb{log}M(z)$), I wondered what tools enable to sum such series and what was Gauss's original method in this case. I tried to make a Google search about "modular forms defined on shifted lattices", but without success.
Sources

*

*Chapter 8 of the book "Groups Acting on Hyperbolic Space: Harmonic Analysis and Number Theory" (p.388-390). On p. 390 there it is mentioned that entry 61 in Gauss's diary suggests that the summation of Eisenstein series of weight 4 (to $\frac{\varpi^4}{15}$) was already known to Gauss in March 1797.

*"Abel's Theorem on the Lemniscate", article by Michael Rosen (1981). In this article there is also (on p.8) a reference to entry 61.

*Gauss's collected works, in particular Klein's and Schlesinger's comments on Gauss's works on the lemniscatic elliptic functions (I have linked to Gauss's diary entry 61 in the beginning of my post). Their comments are necessary to gain insight into Gauss's work since Gauss himself isn't particulary clear on this point and doesn't explain what he is doing.

 A: $(\log M(z)'' = - \varpi^{-2} \wp_i(z/\varpi)$ so Gauss was into the theory elliptic functions.
There is an easy way to show that for $k\ge 4$ $$G_k(i)=\sum_{m,n\ne (0,0)}\frac1{(m+in)^k} \text{ is in }  G_4(i)^{k/4} \Bbb{Q}$$
Note that $G_k(i)=0$ and $G_{2k}(i)=0$ for $k$ odd.
Let
$$\wp_i(z)=\sum_{n,m} \frac1{(m+in+z)^2}-\frac{1_{(m,n)\ne (0,0)}}{(m+in)^2}$$
$\wp_i(z)$ is even, meromorphic, $\Bbb{Z}+i\Bbb{Z}$ periodic. Expand in Laurent series around $0$
$$\wp_i(z) =z^{-2}+\sum_{k\ge 1} \frac{(2k+1)! G_{2k+2}(i) }{(2k)!}z^{2k}= z^{-2}+3 G_4(i) z^2+O(z^6)$$
$$\wp_i''(z)= 6 z^{-4}+6G_4(i)+O(z^2), \qquad \wp_i(z)^2 = z^{-4}+6G_4(i)+O(z^2)$$
So $\wp_i''(z G_4(i))-6\wp_i(z)^2+30 G_4(i)$ is entire, it vanishes at $0$, and since it is $\Bbb{Z}+i\Bbb{Z}$ periodic it must be constant $=0$.

Let $G_4(i)=\frac{\varpi^4}{15}$ and  $f(z)= \varpi^{-2}\wp_i(z/\varpi)$ then
$$f''(z) =6 f(z)^2-2 $$

Equating the Laurent series, you'll find by induction that $\color{red}{\text{the Laurent coefficients of f are all rational}}$. If Gauss had a rigorous proof then it was for sure equivalent to thisone, obtaining a differential equation for something closely related to $\wp_i(z)$.
The relation between $\varpi$ and $G_4(i)$ is obtained from
$$\int_0^1 \frac{dw}{\sqrt{1-w^4}}=\int_0^1 \frac{du}{2 \sqrt{u(1-u^2)}}=\frac12 \int_0^\infty \frac{du}{2 \sqrt{u(1-u^2)}}$$ On the other hand with $g_2(i)=60 G_4(i)$ and $\wp_i'(z)^2=4\wp_i(z)^3-g_2(i)\wp(z)$ (obtained similarly to the previous differential equation)
$$-\frac{1+i}{2} =\int_{(1+i)/2}^0 \frac{d\wp_i(z)}{\wp_i'(z)}=
\int_{(1+i)/2}^0 \frac{d\wp_i(z)}{\sqrt{4\wp_i(z)^4-g_2(i)\wp_i(z)}}$$ $$=
 \int_0^\infty \frac{x}{\sqrt{4x^3-g_2(i)x}}$$
A: Despite not knowing how Gauss arrived at the infinite series for $\mathbb{log}M(z)$, I finally understood how it is connected with summation of Eisenstein series.
First recall the taylor series for $\mathbb{log}(1-z)$ around $z=0$:
$$\mathbb{log}(1-z) = -z - \frac{z^2}{2}-\frac{z^3}{3}-\frac{z^4}{4}-...$$
If we use this for relation (*) in my question, then we get:
$$\mathbb{log}M(z) = \mathbb{log}(z)-\frac{1}{4\varpi^4}\sum\frac{1}{(m+ni)^4}z^4-\frac{1}{8\varpi^8}\sum\frac{1}{(m+ni)^8}z^8-...$$
where coefficients of powers of $z$ not divisible by 4 cancel out because of the modularity of the Eisenstein series. That is, $G_n(\tau)=\sum\frac{1}{(m+n\tau)^n}$ is a modular form of weight $n$ in variable $\tau$, so transforming $\tau = i$ into itself by substituting $\tau' =\frac{-1}{\tau}$ (this transformation belongs to $SL(2,\mathbb{Z})$) transforms the Eisenstein series into $G_n(\tau') = \tau^n G_n(\tau)$. Since $\tau = \tau' = i$, any value of $n$ which is not divisible by 4 leads to the equality $G_n(\tau) = C\cdot G_n(\tau)$ where $C\ne 1$, and this implies $G_n(i) = 0$ for $n\ne 0 \pmod{4}$.
Comparing the previous equality with the infinite series obtained by Gauss in (***), we get:
$$\sum\frac{1}{(m+ni)^4} = \frac{4\varpi^4}{60} = \frac{\varpi^4}{15}, \sum\frac{1}{(m+ni)^8} = \frac{8\varpi^8}{4200} = \frac{\varpi^8}{525}, \sum\frac{1}{(m+ni)^{12}} = \frac{12\varpi^{12}}{321750} = \frac{2\varpi^{12}}{53625}...$$,
in exact accordance with the first known values of sums of Eisenstein series. This partially answers my second question (only partially because i don't know how Gauss arrived at the series (***)).
Update (15 February, 2022)
In order to bring this post into some form of conclusion, the calculation of the infinite series for $\mathbb{log}N(z)$ can be done in the following manner.
Thanks to the answer of user @reuns to my question: How to approach the problem of summation of Eisenstein series on shifted lattices? , which concludes that the infinite series for $\mathbb{log}N(z)$ can be deduced from the identity: $$M(z/2)N(z/2) = M(z/(1+i))$$, one can now take log of both sides and conclude:
$$\mathbb{log}N(z/2) = \mathbb{log}M(z/(1+i))-\mathbb{log}M(z/2)$$,
and since the series for $\mathbb{log}M(z)$ is known, one just needs to make a calculation:
$$\mathbb{log}M(z/(1+i)) = ...-\frac{z^4}{60(1+i)^4}-\frac{z^8}{4200(1+i)^8}-... = ...\frac{z^4}{240}-\frac{z^8}{67200}+...$$
$$\mathbb{log}M(z/2) = ..-\frac{z^4}{960}-\frac{z^8}{1075200}-...$$
so therefore:
$$\mathbb{log}N(z/2) = \frac{z^4}{192}-\frac{z^8}{71680}+...\implies \mathbb{log}N(z) = \frac{z^4\cdot 2^4}{192}-\frac{z^8\cdot 2^8}{71680}+...= \frac{z^4}{12}-\frac{z^8}{280}+...$$
in exact accordance with Gauss's series.
