[proofplan]
We first prove locally uniform absolute convergence on $\mathbb{H}$ by comparing $|mz+n|$ uniformly on compact subsets with the Euclidean length of the lattice vector $(m,n)$. This gives holomorphy on $\mathbb{H}$ by the locally [uniform convergence](/page/Uniform%20Convergence) theorem for holomorphic functions. The modular transformation law follows by rewriting $m\gamma z+n$ and reindexing the absolutely convergent lattice sum using the bijection of $\mathbb{Z}^2$ induced by $\gamma \in SL_2(\mathbb{Z})$. Finally, translation invariance gives a [holomorphic function](/page/Holomorphic%20Function) of $q=e^{2\pi i z}$ on the punctured unit disk, and a direct lattice estimate near $\operatorname{Im}(z)=\infty$ proves boundedness, so the removable singularity theorem gives holomorphy at the cusp.
[/proofplan]
[step:Prove locally uniform absolute convergence on compact subsets of $\mathbb{H}$]
Let $K \subset \mathbb{H}$ be compact. Define
\begin{align*}
\delta_K := \inf \left\{ |uz+v| : z \in K,\ (u,v) \in \mathbb{R}^2,\ u^2+v^2=1 \right\}.
\end{align*}
The set over which the infimum is taken is compact, and the function $(z,u,v) \mapsto |uz+v|$ is continuous. Also $uz+v \neq 0$ on this set: if $u=0$, then $|v|=1$; if $u \neq 0$, then $\operatorname{Im}(uz+v)=u\operatorname{Im}(z) \neq 0$. Hence $\delta_K>0$.
For every $(m,n) \in \mathbb{Z}^2 \setminus \{(0,0)\}$ and every $z \in K$, applying the definition of $\delta_K$ to
\begin{align*}
(u,v) := \frac{(m,n)}{\sqrt{m^2+n^2}}
\end{align*}
gives
\begin{align*}
|mz+n| \geq \delta_K \sqrt{m^2+n^2}.
\end{align*}
Therefore
\begin{align*}
\left|\frac{1}{(mz+n)^k}\right|
\leq
\delta_K^{-k}(m^2+n^2)^{-k/2}.
\end{align*}
Since $k \geq 4$, the lattice series
\begin{align*}
\sum_{(m,n) \in \mathbb{Z}^2 \setminus \{(0,0)\}} (m^2+n^2)^{-k/2}
\end{align*}
converges. Indeed, for $r \in \mathbb{N}$, the number of lattice points satisfying
\begin{align*}
r \leq \sqrt{m^2+n^2} < r+1
\end{align*}
is at most $(2r+3)^2-(2r-1)^2 = 16r+8$, and hence the lattice series is bounded above by
\begin{align*}
\sum_{r=1}^{\infty} (16r+8)r^{-k},
\end{align*}
which converges because $k>2$. The Weierstrass $M$-test gives absolute and uniform convergence on $K$. Since $K$ was arbitrary, the convergence is locally uniform on $\mathbb{H}$.
Each summand
\begin{align*}
z \mapsto (mz+n)^{-k}
\end{align*}
is holomorphic on $\mathbb{H}$, because $mz+n \neq 0$ for $z \in \mathbb{H}$ and $(m,n) \neq (0,0)$. By the locally uniform convergence theorem for holomorphic functions (citing a result not yet in the wiki: Weierstrass theorem for holomorphic functions), the locally uniform sum $G_k$ is holomorphic on $\mathbb{H}$.
[guided]
The first issue is convergence. We need convergence strong enough to preserve holomorphy, so we prove locally uniform absolute convergence. Fix a compact set $K \subset \mathbb{H}$. The quantity $|mz+n|$ must be bounded from below uniformly for $z \in K$ in terms of the size of the lattice vector $(m,n)$. To express this uniformly, define
\begin{align*}
\delta_K := \inf \left\{ |uz+v| : z \in K,\ (u,v) \in \mathbb{R}^2,\ u^2+v^2=1 \right\}.
\end{align*}
The domain of this infimum is compact, because it is the product of $K$ with the unit circle in $\mathbb{R}^2$. The map $(z,u,v) \mapsto |uz+v|$ is continuous. It remains to check that the infimum is positive. If $u=0$, then $v=\pm 1$, so $|uz+v|=1$. If $u \neq 0$, then
\begin{align*}
\operatorname{Im}(uz+v)=u\operatorname{Im}(z),
\end{align*}
which is nonzero because $z \in \mathbb{H}$. Thus $uz+v$ never vanishes on the compact set, and consequently $\delta_K>0$.
Now take any nonzero lattice vector $(m,n) \in \mathbb{Z}^2 \setminus \{(0,0)\}$. The normalized vector
\begin{align*}
(u,v) := \frac{(m,n)}{\sqrt{m^2+n^2}}
\end{align*}
lies on the unit circle in $\mathbb{R}^2$. Therefore, for every $z \in K$,
\begin{align*}
\frac{|mz+n|}{\sqrt{m^2+n^2}} = |uz+v| \geq \delta_K,
\end{align*}
or equivalently
\begin{align*}
|mz+n| \geq \delta_K \sqrt{m^2+n^2}.
\end{align*}
This gives the comparison
\begin{align*}
\left|\frac{1}{(mz+n)^k}\right|
\leq
\delta_K^{-k}(m^2+n^2)^{-k/2}.
\end{align*}
We now verify that the comparison series converges. For $r \in \mathbb{N}$, the number of lattice points in the annulus
\begin{align*}
r \leq \sqrt{m^2+n^2} < r+1
\end{align*}
is at most the number of lattice points in the square $[-r-1,r+1]^2$ minus the number in $[-r+1,r-1]^2$, which is bounded by
\begin{align*}
(2r+3)^2-(2r-1)^2 = 16r+8.
\end{align*}
Thus
\begin{align*}
\sum_{(m,n) \in \mathbb{Z}^2 \setminus \{(0,0)\}} (m^2+n^2)^{-k/2}
\leq
\sum_{r=1}^{\infty} (16r+8)r^{-k}.
\end{align*}
The right-hand side converges because $k>2$, and our hypothesis gives $k \geq 4$. Hence the Weierstrass $M$-test gives absolute and uniform convergence on $K$.
Each summand $z \mapsto (mz+n)^{-k}$ is holomorphic on $\mathbb{H}$: if $m=0$, then $n \neq 0$ and the summand is constant; if $m \neq 0$, then $mz+n=0$ would force $z=-n/m \in \mathbb{R}$, which is impossible for $z \in \mathbb{H}$. Since the series converges locally uniformly and consists of holomorphic functions, the locally uniform convergence theorem for holomorphic functions (citing a result not yet in the wiki: Weierstrass theorem for holomorphic functions) implies that $G_k$ is holomorphic on $\mathbb{H}$.
[/guided]
[/step]
[step:Reindex the lattice sum to prove the modular transformation law]
Let
\begin{align*}
\gamma =
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\in SL_2(\mathbb{Z})
\end{align*}
and let $z \in \mathbb{H}$. Since $cz+d \neq 0$, the point
\begin{align*}
\gamma z := \frac{az+b}{cz+d}
\end{align*}
is well-defined. For every $(m,n) \in \mathbb{Z}^2$, direct algebra gives
\begin{align*}
m\gamma z+n
&=
m\frac{az+b}{cz+d}+n \\
&=
\frac{(ma+nc)z+(mb+nd)}{cz+d}.
\end{align*}
Define the map
\begin{align*}
T_\gamma: \mathbb{Z}^2 &\to \mathbb{Z}^2 \\
(m,n) &\mapsto (ma+nc,\ mb+nd).
\end{align*}
This is right multiplication by the integer matrix $\gamma$. Since $\det(\gamma)=1$, the inverse matrix
\begin{align*}
\gamma^{-1} =
\begin{pmatrix}
d & -b \\
-c & a
\end{pmatrix}
\end{align*}
has integer entries, so $T_\gamma$ is a bijection of $\mathbb{Z}^2$ onto itself and maps $\mathbb{Z}^2 \setminus \{(0,0)\}$ onto itself.
Using absolute convergence to reindex the series,
\begin{align*}
G_k(\gamma z)
&=
\sum_{(m,n) \in \mathbb{Z}^2 \setminus \{(0,0)\}}
\frac{1}{(m\gamma z+n)^k} \\
&=
\sum_{(m,n) \in \mathbb{Z}^2 \setminus \{(0,0)\}}
\frac{(cz+d)^k}{\left((ma+nc)z+(mb+nd)\right)^k} \\
&=
(cz+d)^k
\sum_{(m',n') \in \mathbb{Z}^2 \setminus \{(0,0)\}}
\frac{1}{(m'z+n')^k} \\
&=
(cz+d)^kG_k(z).
\end{align*}
This proves the transformation law for every $\gamma \in SL_2(\mathbb{Z})$.
[guided]
Fix
\begin{align*}
\gamma =
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\in SL_2(\mathbb{Z})
\end{align*}
and $z \in \mathbb{H}$. First, $cz+d \neq 0$: if $c=0$, then $d=\pm 1$ because $ad-bc=1$; if $c \neq 0$ and $cz+d=0$, then $z=-d/c \in \mathbb{R}$, contradicting $z \in \mathbb{H}$. Thus
\begin{align*}
\gamma z := \frac{az+b}{cz+d}
\end{align*}
is defined.
We compute how one denominator changes under $\gamma$. For any $(m,n) \in \mathbb{Z}^2$,
\begin{align*}
m\gamma z+n
&=
m\frac{az+b}{cz+d}+n \\
&=
\frac{m(az+b)+n(cz+d)}{cz+d} \\
&=
\frac{(ma+nc)z+(mb+nd)}{cz+d}.
\end{align*}
The new pair of coefficients is
\begin{align*}
(ma+nc,\ mb+nd).
\end{align*}
This is not an arbitrary change of variables: it is the map
\begin{align*}
T_\gamma: \mathbb{Z}^2 &\to \mathbb{Z}^2 \\
(m,n) &\mapsto (ma+nc,\ mb+nd),
\end{align*}
namely right multiplication of the row vector $(m,n)$ by the integer matrix $\gamma$. Because $\gamma \in SL_2(\mathbb{Z})$, its determinant is $1$, and its inverse is
\begin{align*}
\gamma^{-1} =
\begin{pmatrix}
d & -b \\
-c & a
\end{pmatrix},
\end{align*}
which again has integer entries. Hence $T_\gamma$ is a bijection from $\mathbb{Z}^2$ to $\mathbb{Z}^2$, and it sends the nonzero lattice vectors bijectively to the nonzero lattice vectors.
Now absolute convergence matters: it allows us to reindex the infinite series without changing its value. We obtain
\begin{align*}
G_k(\gamma z)
&=
\sum_{(m,n) \in \mathbb{Z}^2 \setminus \{(0,0)\}}
\frac{1}{(m\gamma z+n)^k} \\
&=
\sum_{(m,n) \in \mathbb{Z}^2 \setminus \{(0,0)\}}
\frac{(cz+d)^k}{\left((ma+nc)z+(mb+nd)\right)^k}.
\end{align*}
Since $(cz+d)^k$ is independent of $(m,n)$, factor it out:
\begin{align*}
G_k(\gamma z)
&=
(cz+d)^k
\sum_{(m,n) \in \mathbb{Z}^2 \setminus \{(0,0)\}}
\frac{1}{\left((ma+nc)z+(mb+nd)\right)^k}.
\end{align*}
Finally set
\begin{align*}
(m',n') := (ma+nc,\ mb+nd).
\end{align*}
Because $T_\gamma$ is a bijection of $\mathbb{Z}^2 \setminus \{(0,0)\}$, the last sum is exactly
\begin{align*}
\sum_{(m',n') \in \mathbb{Z}^2 \setminus \{(0,0)\}}
\frac{1}{(m'z+n')^k}
=
G_k(z).
\end{align*}
Therefore
\begin{align*}
G_k(\gamma z)=(cz+d)^kG_k(z),
\end{align*}
which is the modular transformation law of weight $k$.
[/guided]
[/step]
[step:Show that $G_k$ is holomorphic at the cusp $\infty$]
Apply the transformation law to
\begin{align*}
T :=
\begin{pmatrix}
1 & 1 \\
0 & 1
\end{pmatrix}
\in SL_2(\mathbb{Z}).
\end{align*}
Since $Tz=z+1$ and $cz+d=1$, we have
\begin{align*}
G_k(z+1)=G_k(z)
\end{align*}
for every $z \in \mathbb{H}$.
Define the punctured unit disk
\begin{align*}
\mathbb{D}^* := \{q \in \mathbb{C} : 0<|q|<1\}.
\end{align*}
The periodicity of $G_k$ makes the following map well-defined:
\begin{align*}
g: \mathbb{D}^* &\to \mathbb{C} \\
e^{2\pi iz} &\mapsto G_k(z).
\end{align*}
It is holomorphic because locally one may choose a holomorphic branch of
\begin{align*}
z = \frac{1}{2\pi i}\log q.
\end{align*}
It remains to prove that $g$ is bounded near $q=0$, equivalently that $G_k$ is bounded on the vertical strip
\begin{align*}
S := \{x+iy \in \mathbb{H} : 0 \leq x \leq 1,\ y \geq 1\}.
\end{align*}
Let $z=x+iy \in S$. The terms with $m=0$ contribute
\begin{align*}
\sum_{n \in \mathbb{Z}\setminus\{0\}} \frac{1}{n^k}
=
2\sum_{n=1}^{\infty}\frac{1}{n^k},
\end{align*}
because $k$ is even.
For the terms with $m \neq 0$, define
\begin{align*}
C_k := 6+\frac{2}{k-1}.
\end{align*}
We claim that for every $a \in \mathbb{R}$ and every $A \geq 1$,
\begin{align*}
\sum_{n \in \mathbb{Z}} \frac{1}{\left((n+a)^2+A^2\right)^{k/2}}
\leq
C_k A^{1-k}.
\end{align*}
Indeed, for each integer $j \geq 0$, at most two integers $n$ satisfy
\begin{align*}
j \leq |n+a| < j+1.
\end{align*}
Thus
\begin{align*}
\sum_{n \in \mathbb{Z}} \frac{1}{\left((n+a)^2+A^2\right)^{k/2}}
&\leq
2\sum_{j=0}^{\infty} \frac{1}{(j^2+A^2)^{k/2}} \\
&\leq
2\sum_{0 \leq j \leq A} A^{-k}
+
2\sum_{j>A} j^{-k}.
\end{align*}
Since $A \geq 1$, the first sum is at most $4A^{1-k}$. Let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$. For the second sum, comparison with the decreasing function $t \mapsto t^{-k}$ gives
\begin{align*}
\sum_{j>A} j^{-k}
\leq
A^{-k}+\int_A^\infty t^{-k}\,d\mathcal{L}^1(t)
=
A^{-k}+\frac{A^{1-k}}{k-1}
\leq
\left(1+\frac{1}{k-1}\right)A^{1-k}.
\end{align*}
Combining the two estimates gives the claimed bound with $C_k=6+2/(k-1)$.
Now put $a=mx$ and $A=|m|y$. Since $m \neq 0$ and $y \geq 1$, we have $A \geq 1$, and hence
\begin{align*}
\sum_{n \in \mathbb{Z}} \left|\frac{1}{(mz+n)^k}\right|
&=
\sum_{n \in \mathbb{Z}}
\frac{1}{\left((mx+n)^2+(my)^2\right)^{k/2}} \\
&\leq
C_k |m|^{1-k} y^{1-k}.
\end{align*}
Therefore
\begin{align*}
\sum_{m \in \mathbb{Z}\setminus\{0\}}\sum_{n \in \mathbb{Z}}
\left|\frac{1}{(mz+n)^k}\right|
&\leq
C_k y^{1-k}
\sum_{m \in \mathbb{Z}\setminus\{0\}} |m|^{1-k} \\
&\leq
2C_k\sum_{m=1}^{\infty} m^{1-k},
\end{align*}
because $y^{1-k}\leq 1$ for $y \geq 1$ and $k \geq 4$. The right-hand side is finite. Thus $G_k$ is bounded on $S$.
Hence $g$ is bounded in a punctured neighbourhood of $0$. By the removable singularity theorem for holomorphic functions (citing a result not yet in the wiki: removable singularity theorem), $g$ extends holomorphically to $q=0$. This is precisely holomorphy of $G_k$ at the cusp $\infty$.
[guided]
To prove holomorphy at the cusp, we first need a function of the cusp parameter $q=e^{2\pi iz}$. The matrix
\begin{align*}
T :=
\begin{pmatrix}
1 & 1 \\
0 & 1
\end{pmatrix}
\end{align*}
lies in $SL_2(\mathbb{Z})$. Applying the modular transformation law already proved gives
\begin{align*}
G_k(Tz)=G_k(z+1)=G_k(z),
\end{align*}
because for $T$ we have $c=0$ and $d=1$. Thus $G_k$ is $1$-periodic.
Define
\begin{align*}
\mathbb{D}^* := \{q \in \mathbb{C} : 0<|q|<1\}.
\end{align*}
The map $z \mapsto q=e^{2\pi iz}$ sends $\mathbb{H}$ onto $\mathbb{D}^*$. Because $G_k(z+1)=G_k(z)$, the value of $G_k(z)$ depends only on $q$, not on the chosen logarithm of $q$. Hence the map
\begin{align*}
g: \mathbb{D}^* &\to \mathbb{C} \\
e^{2\pi iz} &\mapsto G_k(z)
\end{align*}
is well-defined. It is holomorphic on $\mathbb{D}^*$ because around every point of $\mathbb{D}^*$ there is a holomorphic branch of $\log q$, and locally
\begin{align*}
g(q)=G_k\left(\frac{1}{2\pi i}\log q\right).
\end{align*}
Now we prove that $g$ is bounded near $q=0$. Since
\begin{align*}
|q|=|e^{2\pi i(x+iy)}|=e^{-2\pi y},
\end{align*}
the limit $q \to 0$ corresponds to $y=\operatorname{Im}(z)\to\infty$. Because $G_k$ is $1$-periodic, it is enough to bound it on the strip
\begin{align*}
S := \{x+iy \in \mathbb{H} : 0 \leq x \leq 1,\ y \geq 1\}.
\end{align*}
Fix $z=x+iy \in S$. Separate the Eisenstein series into the terms with $m=0$ and the terms with $m \neq 0$. When $m=0$, the denominator is $n^k$, and since $k$ is even,
\begin{align*}
\sum_{n \in \mathbb{Z}\setminus\{0\}} \frac{1}{n^k}
=
2\sum_{n=1}^{\infty}\frac{1}{n^k},
\end{align*}
which is finite because $k>1$.
For $m \neq 0$, we need a uniform estimate for the sum over $n$. Define
\begin{align*}
C_k := 6+\frac{2}{k-1}.
\end{align*}
We claim that for every real number $a$ and every $A \geq 1$,
\begin{align*}
\sum_{n \in \mathbb{Z}} \frac{1}{\left((n+a)^2+A^2\right)^{k/2}}
\leq
C_k A^{1-k}.
\end{align*}
The idea is that the sum over $n$ behaves like the integral of $(t^2+A^2)^{-k/2}$, whose size is $A^{1-k}$. To make this discrete comparison explicit, group integers according to the size of $|n+a|$. For each $j \geq 0$, at most two integers $n$ satisfy
\begin{align*}
j \leq |n+a| < j+1.
\end{align*}
Therefore
\begin{align*}
\sum_{n \in \mathbb{Z}} \frac{1}{\left((n+a)^2+A^2\right)^{k/2}}
&\leq
2\sum_{j=0}^{\infty} \frac{1}{(j^2+A^2)^{k/2}} \\
&\leq
2\sum_{0 \leq j \leq A} A^{-k}
+
2\sum_{j>A} j^{-k}.
\end{align*}
The first part has at most $A+1$ integer values of $j$, and since $A \geq 1$,
\begin{align*}
2\sum_{0 \leq j \leq A} A^{-k}
\leq
4A^{1-k}.
\end{align*}
Let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$. For the tail, compare the decreasing function $t \mapsto t^{-k}$ with its integral:
\begin{align*}
\sum_{j>A} j^{-k}
\leq
A^{-k}+\int_A^\infty t^{-k}\,d\mathcal{L}^1(t)
=
A^{-k}+\frac{A^{1-k}}{k-1}
\leq
\left(1+\frac{1}{k-1}\right)A^{1-k}.
\end{align*}
Combining the two pieces gives the stated bound with $C_k=6+2/(k-1)$.
Apply this estimate with
\begin{align*}
a := mx,
\qquad
A := |m|y.
\end{align*}
Because $m \neq 0$ and $y \geq 1$, we have $A \geq 1$. Also
\begin{align*}
|mz+n|^2=(mx+n)^2+(my)^2.
\end{align*}
Thus
\begin{align*}
\sum_{n \in \mathbb{Z}} \left|\frac{1}{(mz+n)^k}\right|
&=
\sum_{n \in \mathbb{Z}}
\frac{1}{\left((mx+n)^2+(my)^2\right)^{k/2}} \\
&\leq
C_k |m|^{1-k} y^{1-k}.
\end{align*}
Now sum over all nonzero $m$:
\begin{align*}
\sum_{m \in \mathbb{Z}\setminus\{0\}}\sum_{n \in \mathbb{Z}}
\left|\frac{1}{(mz+n)^k}\right|
&\leq
C_k y^{1-k}
\sum_{m \in \mathbb{Z}\setminus\{0\}} |m|^{1-k}.
\end{align*}
Since $y \geq 1$ and $k \geq 4$, we have $y^{1-k} \leq 1$, and
\begin{align*}
\sum_{m \in \mathbb{Z}\setminus\{0\}} |m|^{1-k}
=
2\sum_{m=1}^{\infty} m^{1-k}
\end{align*}
converges because $k-1>1$. Therefore the entire Eisenstein series is bounded uniformly on $S$.
This boundedness says exactly that $g(q)$ is bounded for $q$ sufficiently close to $0$ with $q \neq 0$. Since $g$ is holomorphic on the punctured disk $\mathbb{D}^*$, the removable singularity theorem for holomorphic functions (citing a result not yet in the wiki: removable singularity theorem) gives a holomorphic extension of $g$ to $q=0$. This is the definition of holomorphy of $G_k$ at the cusp $\infty$.
[/guided]
[/step]
[step:Conclude that $G_k$ is a holomorphic modular form of weight $k$]
We have proved that $G_k$ is holomorphic on $\mathbb{H}$, satisfies
\begin{align*}
G_k\left(\frac{az+b}{cz+d}\right)=(cz+d)^kG_k(z)
\end{align*}
for every
\begin{align*}
\begin{pmatrix}
a & b \\
c & d
\end{pmatrix}
\in SL_2(\mathbb{Z}),
\end{align*}
and is holomorphic at the cusp $\infty$. Let $M_k(SL_2(\mathbb{Z}))$ denote the complex [vector space](/page/Vector%20Space) of holomorphic modular forms of weight $k$ for $SL_2(\mathbb{Z})$. The three properties just proved are precisely the defining conditions for membership in this space. Therefore $G_k \in M_k(SL_2(\mathbb{Z}))$.
[/step]