[proofplan]
We first justify that the defining lattice series for $G_k$ is absolutely convergent, so that the terms may be grouped by rows indexed by $m \in \mathbb{Z}$. The row $m=0$ gives the constant term $2\zeta(k)$, while the positive rows are evaluated using the classical Fourier expansion of the one-dimensional row sum. The negative rows equal the positive rows because $k$ is even. Finally, we group the resulting double series by the product $N=mr$ and divide by $2\zeta(k)$, using [Euler's formula](/theorems/2014) for $\zeta(k)$ in terms of the Bernoulli number $B_k$.
[/proofplan]
[step:Justify grouping the Eisenstein lattice series by horizontal rows]
Fix $z \in \mathbb{H}$, and write $z = x + iy$ with $x,y \in \mathbb{R}$ and $y > 0$. Define the real-[linear map](/page/Linear%20Map)
\begin{align*}
T_z: \mathbb{R}^2 &\to \mathbb{C} \\
(u,v) &\mapsto uz+v.
\end{align*}
Identifying $\mathbb{C}$ with $\mathbb{R}^2$, the map $T_z$ is invertible because its matrix in the standard real bases is
\begin{align*}
\begin{pmatrix}
x & 1 \\
y & 0
\end{pmatrix},
\end{align*}
whose determinant is $-y \neq 0$. Hence, by equivalence of norms on the finite-dimensional space $\mathbb{R}^2$, there exists a constant $c_z>0$ such that
\begin{align*}
|uz+v| \geq c_z (u^2+v^2)^{1/2}
\end{align*}
for all $(u,v)\in\mathbb{R}^2$.
Therefore, for every $(m,n)\in\mathbb{Z}^2\setminus\{(0,0)\}$,
\begin{align*}
\left|\frac{1}{(mz+n)^k}\right|
\leq c_z^{-k}(m^2+n^2)^{-k/2}.
\end{align*}
Since $k\geq4$, the comparison series
\begin{align*}
\sum_{(m,n)\in\mathbb{Z}^2\setminus\{(0,0)\}}(m^2+n^2)^{-k/2}
\end{align*}
converges. Indeed, grouping lattice points in annuli
\begin{align*}
A_j := \{(m,n)\in\mathbb{Z}^2 : j \leq (m^2+n^2)^{1/2} < j+1\},
\end{align*}
there is a constant $C_0>0$ such that $\#A_j \leq C_0 j$ for all $j\geq1$, and hence
\begin{align*}
\sum_{(m,n)\neq(0,0)}(m^2+n^2)^{-k/2}
\leq C_0\sum_{j=1}^{\infty} j^{1-k},
\end{align*}
which converges because $k>2$. Thus $G_k(z)$ is absolutely convergent, and we may group the terms according to the value of $m$.
[guided]
Fix $z \in \mathbb{H}$, and write $z = x+iy$ with $x,y\in\mathbb{R}$ and $y>0$. Before separating the lattice sum into rows, we must justify that such regrouping is legitimate. Absolute convergence gives this justification.
Define the real-linear map
\begin{align*}
T_z: \mathbb{R}^2 &\to \mathbb{C} \\
(u,v) &\mapsto uz+v.
\end{align*}
Under the identification $\mathbb{C}\cong\mathbb{R}^2$, this map has matrix
\begin{align*}
\begin{pmatrix}
x & 1 \\
y & 0
\end{pmatrix}.
\end{align*}
Its determinant is $-y$, which is nonzero because $z\in\mathbb{H}$. Hence $T_z$ is an invertible real-linear map. Since all norms on the finite-dimensional space $\mathbb{R}^2$ are equivalent, there is a constant $c_z>0$ such that
\begin{align*}
|uz+v| = |T_z(u,v)| \geq c_z (u^2+v^2)^{1/2}
\end{align*}
for every $(u,v)\in\mathbb{R}^2$.
Applying this estimate to integer pairs $(m,n)\neq(0,0)$ gives
\begin{align*}
\left|\frac{1}{(mz+n)^k}\right|
\leq c_z^{-k}(m^2+n^2)^{-k/2}.
\end{align*}
It remains to know that the right-hand comparison series converges. For each $j\in\mathbb{N}$, define the lattice annulus
\begin{align*}
A_j := \{(m,n)\in\mathbb{Z}^2 : j \leq (m^2+n^2)^{1/2} < j+1\}.
\end{align*}
The number of lattice points in $A_j$ is bounded above by $C_0j$ for a constant $C_0>0$, because the annulus has radius comparable to $j$ and bounded width. Thus
\begin{align*}
\sum_{(m,n)\neq(0,0)}(m^2+n^2)^{-k/2}
\leq C_0\sum_{j=1}^{\infty} j\cdot j^{-k}
= C_0\sum_{j=1}^{\infty}j^{1-k}.
\end{align*}
This final series converges since $k>2$. Therefore the Eisenstein lattice series is absolutely convergent, and separating it into the row $m=0$, the rows $m>0$, and the rows $m<0$ is valid.
[/guided]
[/step]
[step:Evaluate the row with $m=0$]
The contribution of the row $m=0$ is
\begin{align*}
\sum_{n\in\mathbb{Z}\setminus\{0\}}\frac{1}{n^k}
= \sum_{n=1}^{\infty}\frac{1}{n^k}+\sum_{n=1}^{\infty}\frac{1}{(-n)^k}.
\end{align*}
Since $k$ is even, $(-n)^k=n^k$, and hence this row contributes
\begin{align*}
2\sum_{n=1}^{\infty}\frac{1}{n^k}=2\zeta(k).
\end{align*}
[/step]
[step:Evaluate the positive rows using the row sum Fourier expansion]
For $w\in\mathbb{H}$, use the classical row sum expansion
(citing a result not yet in the wiki: Row Sum Fourier Expansion)
\begin{align*}
\sum_{n\in\mathbb{Z}}\frac{1}{(w+n)^k}
=
\frac{(-2\pi i)^k}{(k-1)!}\sum_{r=1}^{\infty}r^{k-1}e^{2\pi i r w}.
\end{align*}
For each $m\in\mathbb{N}$, the point $mz$ lies in $\mathbb{H}$ because $\operatorname{Im}(mz)=m\operatorname{Im}(z)>0$. Therefore the row sum expansion applies with $w=mz$, giving
\begin{align*}
\sum_{n\in\mathbb{Z}}\frac{1}{(mz+n)^k}
=
\frac{(-2\pi i)^k}{(k-1)!}\sum_{r=1}^{\infty}r^{k-1}e^{2\pi i r m z}.
\end{align*}
Since $q=e^{2\pi iz}$, this becomes
\begin{align*}
\sum_{n\in\mathbb{Z}}\frac{1}{(mz+n)^k}
=
\frac{(-2\pi i)^k}{(k-1)!}\sum_{r=1}^{\infty}r^{k-1}q^{mr}.
\end{align*}
Summing over $m\geq1$, the total contribution of the positive rows is
\begin{align*}
\frac{(-2\pi i)^k}{(k-1)!}
\sum_{m=1}^{\infty}\sum_{r=1}^{\infty}r^{k-1}q^{mr}.
\end{align*}
[guided]
For the positive rows, we fix $m\in\mathbb{N}$ and regard the inner sum over $n\in\mathbb{Z}$ as a one-dimensional periodic row sum. The theorem we use is the classical row sum Fourier expansion
(citing a result not yet in the wiki: Row Sum Fourier Expansion):
\begin{align*}
\sum_{n\in\mathbb{Z}}\frac{1}{(w+n)^k}
=
\frac{(-2\pi i)^k}{(k-1)!}\sum_{r=1}^{\infty}r^{k-1}e^{2\pi i r w},
\end{align*}
valid for $w\in\mathbb{H}$ and integer $k\geq2$.
We must verify its hypothesis. Since $z\in\mathbb{H}$ and $m\in\mathbb{N}$, we have
\begin{align*}
\operatorname{Im}(mz)=m\operatorname{Im}(z)>0.
\end{align*}
Thus $mz\in\mathbb{H}$, so the row sum expansion applies with $w=mz$. We obtain
\begin{align*}
\sum_{n\in\mathbb{Z}}\frac{1}{(mz+n)^k}
=
\frac{(-2\pi i)^k}{(k-1)!}\sum_{r=1}^{\infty}r^{k-1}e^{2\pi i r m z}.
\end{align*}
Now define $q:=e^{2\pi iz}$. Then
\begin{align*}
e^{2\pi i r m z}=q^{mr},
\end{align*}
so the $m$-th positive row is
\begin{align*}
\sum_{n\in\mathbb{Z}}\frac{1}{(mz+n)^k}
=
\frac{(-2\pi i)^k}{(k-1)!}\sum_{r=1}^{\infty}r^{k-1}q^{mr}.
\end{align*}
Summing this identity over all positive integers $m$ gives the total positive-row contribution
\begin{align*}
\frac{(-2\pi i)^k}{(k-1)!}
\sum_{m=1}^{\infty}\sum_{r=1}^{\infty}r^{k-1}q^{mr}.
\end{align*}
[/guided]
[/step]
[step:Show that the negative rows contribute the same amount]
Let $a\in\mathbb{N}$. The row with $m=-a$ is
\begin{align*}
\sum_{n\in\mathbb{Z}}\frac{1}{(-az+n)^k}.
\end{align*}
Changing the summation index by $n=-\ell$, where $\ell\in\mathbb{Z}$, gives
\begin{align*}
\sum_{n\in\mathbb{Z}}\frac{1}{(-az+n)^k}
=
\sum_{\ell\in\mathbb{Z}}\frac{1}{(-az-\ell)^k}
=
(-1)^{-k}\sum_{\ell\in\mathbb{Z}}\frac{1}{(az+\ell)^k}.
\end{align*}
Since $k$ is even, $(-1)^{-k}=1$. Therefore
\begin{align*}
\sum_{n\in\mathbb{Z}}\frac{1}{(-az+n)^k}
=
\sum_{\ell\in\mathbb{Z}}\frac{1}{(az+\ell)^k}.
\end{align*}
Thus the negative rows have the same total contribution as the positive rows.
[/step]
[step:Group the double series by the product $N=mr$]
Combining the row computations gives
\begin{align*}
G_k(z)
=
2\zeta(k)
+
2\frac{(-2\pi i)^k}{(k-1)!}
\sum_{m=1}^{\infty}\sum_{r=1}^{\infty}r^{k-1}q^{mr}.
\end{align*}
Because $z\in\mathbb{H}$, we have $|q|=e^{-2\pi\operatorname{Im}(z)}<1$. Hence
\begin{align*}
\sum_{m=1}^{\infty}\sum_{r=1}^{\infty}r^{k-1}|q|^{mr}
=
\sum_{r=1}^{\infty}r^{k-1}\frac{|q|^r}{1-|q|^r}
\leq
\frac{1}{1-|q|}
\sum_{r=1}^{\infty}r^{k-1}|q|^r
<\infty.
\end{align*}
Thus the double series is absolutely convergent, and we may group terms according to $N=mr$. For each $N\in\mathbb{N}$, the pairs $(m,r)\in\mathbb{N}^2$ with $mr=N$ are exactly the pairs $(N/d,d)$ where $d$ runs over the positive divisors of $N$. Therefore
\begin{align*}
\sum_{m=1}^{\infty}\sum_{r=1}^{\infty}r^{k-1}q^{mr}
=
\sum_{N=1}^{\infty}\left(\sum_{d\mid N}d^{k-1}\right)q^N
=
\sum_{N=1}^{\infty}\sigma_{k-1}(N)q^N.
\end{align*}
Consequently,
\begin{align*}
G_k(z)
=
2\zeta(k)
+
2\frac{(-2\pi i)^k}{(k-1)!}
\sum_{N=1}^{\infty}\sigma_{k-1}(N)q^N.
\end{align*}
[guided]
At this point we have computed the three pieces of the lattice sum: the row $m=0$, the positive rows, and the negative rows. Putting them together gives
\begin{align*}
G_k(z)
=
2\zeta(k)
+
2\frac{(-2\pi i)^k}{(k-1)!}
\sum_{m=1}^{\infty}\sum_{r=1}^{\infty}r^{k-1}q^{mr}.
\end{align*}
The factor $2$ appears because the negative rows equal the positive rows.
We now want to rewrite the double series as a single [Fourier series](/page/Fourier%20Series) in powers of $q$. Since $z\in\mathbb{H}$,
\begin{align*}
|q|=|e^{2\pi iz}|=e^{-2\pi\operatorname{Im}(z)}<1.
\end{align*}
This exponential decay justifies rearranging the double series. Indeed,
\begin{align*}
\sum_{m=1}^{\infty}\sum_{r=1}^{\infty}r^{k-1}|q|^{mr}
=
\sum_{r=1}^{\infty}r^{k-1}\sum_{m=1}^{\infty}|q|^{mr}
=
\sum_{r=1}^{\infty}r^{k-1}\frac{|q|^r}{1-|q|^r}.
\end{align*}
Since $1-|q|^r\geq 1-|q|$ for every $r\geq1$, we have
\begin{align*}
\sum_{r=1}^{\infty}r^{k-1}\frac{|q|^r}{1-|q|^r}
\leq
\frac{1}{1-|q|}
\sum_{r=1}^{\infty}r^{k-1}|q|^r.
\end{align*}
The final series converges because polynomial growth is dominated by exponential decay. Hence the double series is absolutely convergent.
We may now group terms by the exponent $N=mr$. For a fixed $N\in\mathbb{N}$, the condition $mr=N$ says that $r$ is a positive divisor of $N$ and $m=N/r$. Thus
\begin{align*}
\sum_{\substack{m,r\in\mathbb{N}\\ mr=N}}r^{k-1}
=
\sum_{d\mid N}d^{k-1}
=
\sigma_{k-1}(N).
\end{align*}
Therefore
\begin{align*}
\sum_{m=1}^{\infty}\sum_{r=1}^{\infty}r^{k-1}q^{mr}
=
\sum_{N=1}^{\infty}\sigma_{k-1}(N)q^N,
\end{align*}
and so
\begin{align*}
G_k(z)
=
2\zeta(k)
+
2\frac{(-2\pi i)^k}{(k-1)!}
\sum_{N=1}^{\infty}\sigma_{k-1}(N)q^N.
\end{align*}
[/guided]
[/step]
[step:Normalize and convert the scalar using Bernoulli numbers]
By definition,
\begin{align*}
E_k(z)=\frac{G_k(z)}{2\zeta(k)}.
\end{align*}
Using the formula for $G_k(z)$ obtained above,
\begin{align*}
E_k(z)
=
1+
\frac{(-2\pi i)^k}{(k-1)!\zeta(k)}
\sum_{N=1}^{\infty}\sigma_{k-1}(N)q^N.
\end{align*}
Since $k$ is even, write $k=2\ell$ with $\ell\in\mathbb{N}$ and $\ell\geq2$. Then
\begin{align*}
(-2\pi i)^k
=
(-2\pi i)^{2\ell}
=
(-1)^{\ell}(2\pi)^{2\ell}
=
(-1)^{k/2}(2\pi)^k.
\end{align*}
We use Euler's evaluation of the zeta value at an even positive integer
(citing a result not yet in the wiki: Euler's Formula for Even Zeta Values)
\begin{align*}
\zeta(k)=(-1)^{k/2+1}\frac{B_k(2\pi)^k}{2k!}.
\end{align*}
Therefore
\begin{align*}
\frac{(-2\pi i)^k}{(k-1)!\zeta(k)}
&=
\frac{(-1)^{k/2}(2\pi)^k}{(k-1)!}
\cdot
\frac{2k!}{(-1)^{k/2+1}B_k(2\pi)^k} \\
&=
-\frac{2k!}{(k-1)!B_k} \\
&=
-\frac{2k}{B_k}.
\end{align*}
Substituting this scalar into the expression for $E_k(z)$ yields
\begin{align*}
E_k(z)
=
1-\frac{2k}{B_k}\sum_{N=1}^{\infty}\sigma_{k-1}(N)q^N.
\end{align*}
This is the desired Fourier expansion.
[/step]