[proofplan]
We use the double-coset definition of the Hecke operator. The modular transformation law follows because right multiplication by an element of $\Gamma$ permutes the left cosets in $\Gamma \backslash \Delta_n$, and the weight-$k$ slash action is invariant under left multiplication by $\Gamma$ when applied to a modular form. Holomorphy on the upper half-plane is immediate from the fact that matrices in $\Delta_n$ preserve $\mathbb{H}$. Finally, holomorphy at the cusp is checked from the standard triangular representatives, which give an explicit Fourier coefficient formula with no negative powers of $q$.
[/proofplan]
[step:Declare the slash action and the Hecke operator]
Let $\mathbb{C}$ denote the field of complex numbers, let $\mathbb{R}\subset\mathbb{C}$ denote the real line, let $\mathbb{Z}$ denote the ring of integers, and let
\begin{align*}
\mathbb{H}:=\{z\in\mathbb{C}:\operatorname{Im}z>0\}
\end{align*}
denote the complex upper half-plane. Let $\Gamma:=SL_2(\mathbb{Z})$, and define
\begin{align*}
\Delta_n := \{\alpha \in M_2(\mathbb{Z}) : \det \alpha = n\}.
\end{align*}
For a matrix
\begin{align*}
\alpha =
\begin{pmatrix}
a & b\\
c & d
\end{pmatrix}
\in \Delta_n
\end{align*}
and a [holomorphic function](/page/Holomorphic%20Function) $f:\mathbb{H}\to \mathbb{C}$, define the weight-$k$ slash action
\begin{align*}
(f|_k \alpha): \mathbb{H} &\to \mathbb{C}\\
z &\mapsto \det(\alpha)^{k-1}(cz+d)^{-k} f\left(\frac{az+b}{cz+d}\right).
\end{align*}
The denominator $cz+d$ is nonzero for $z \in \mathbb{H}$. Indeed, if $cz+d=0$ and $c\ne 0$, then
\begin{align*}
z=-\frac{d}{c}\in\mathbb{R},
\end{align*}
contradicting $z\in\mathbb{H}$. If $c=0$, then $d\ne 0$ because $\det \alpha=n>0$.
Choose a finite set $R_n \subset \Delta_n$ of representatives for the left cosets $\Gamma\backslash \Delta_n$. Define
\begin{align*}
(T_n f): \mathbb{H} &\to \mathbb{C}\\
z &\mapsto \sum_{\alpha \in R_n} (f|_k \alpha)(z).
\end{align*}
Here $M_k(\Gamma)$ denotes the space of holomorphic functions $f:\mathbb{H}\to\mathbb{C}$ satisfying $f|_k\gamma=f$ for every $\gamma\in\Gamma$ and holomorphicity at the cusp, meaning that $f$ has a Fourier expansion in nonnegative powers of $q=e^{2\pi i z}$ for $\operatorname{Im}z$ sufficiently large. This definition is independent of the chosen representatives: if $\alpha'=\gamma\alpha$ with $\gamma\in\Gamma$, then
\begin{align*}
f|_k\alpha' = f|_k(\gamma\alpha) = (f|_k\gamma)|_k\alpha = f|_k\alpha,
\end{align*}
because $f\in M_k(\Gamma)$ satisfies $f|_k\gamma=f$ for every $\gamma\in\Gamma$.
[guided]
We first fix the precise normalization of the operator and the ambient notation. Let $\mathbb{C}$ denote the field of complex numbers, let $\mathbb{R}\subset\mathbb{C}$ denote the real line, let $\mathbb{Z}$ denote the ring of integers, and let
\begin{align*}
\mathbb{H}:=\{z\in\mathbb{C}:\operatorname{Im}z>0\}
\end{align*}
denote the complex upper half-plane. Let $\Gamma:=SL_2(\mathbb{Z})$. The relevant matrices are
\begin{align*}
\Delta_n := \{\alpha \in M_2(\mathbb{Z}) : \det \alpha = n\}.
\end{align*}
Each such matrix has positive determinant, so it acts on the upper half-plane $\mathbb{H}$ by the usual fractional linear transformation. If
\begin{align*}
\alpha =
\begin{pmatrix}
a & b\\
c & d
\end{pmatrix},
\end{align*}
we define
\begin{align*}
(f|_k \alpha): \mathbb{H} &\to \mathbb{C}\\
z &\mapsto \det(\alpha)^{k-1}(cz+d)^{-k} f\left(\frac{az+b}{cz+d}\right).
\end{align*}
The expression is well-defined on $\mathbb{H}$. Indeed, if $cz+d=0$ and $c\ne 0$, then
\begin{align*}
z=-\frac{d}{c}\in\mathbb{R},
\end{align*}
contradicting $z\in\mathbb{H}$. If $c=0$, then $d\ne 0$ because $\det\alpha=ad=n>0$.
Now choose a finite set $R_n\subset\Delta_n$ representing the left cosets in $\Gamma\backslash\Delta_n$, and define
\begin{align*}
(T_n f): \mathbb{H} &\to \mathbb{C}\\
z &\mapsto \sum_{\alpha\in R_n}(f|_k\alpha)(z).
\end{align*}
Here $M_k(\Gamma)$ denotes the space of holomorphic functions $f:\mathbb{H}\to\mathbb{C}$ satisfying $f|_k\gamma=f$ for every $\gamma\in\Gamma$ and holomorphicity at the cusp, meaning that $f$ has a Fourier expansion in nonnegative powers of $q=e^{2\pi i z}$ for $\operatorname{Im}z$ sufficiently large. Why does this not depend on the choice of representative? If $\alpha'$ is another representative of the same left coset, then $\alpha'=\gamma\alpha$ for some $\gamma\in\Gamma$. The slash action satisfies associativity:
\begin{align*}
f|_k(\gamma\alpha)=(f|_k\gamma)|_k\alpha.
\end{align*}
Since $f$ is a modular form of weight $k$ for $\Gamma$, it satisfies $f|_k\gamma=f$ for every $\gamma\in\Gamma$. Therefore
\begin{align*}
f|_k\alpha'=f|_k(\gamma\alpha)=(f|_k\gamma)|_k\alpha=f|_k\alpha.
\end{align*}
Thus the coset sum is intrinsic.
[/guided]
[/step]
[step:Show the Hecke sum has the modular transformation law]
Let $\gamma\in\Gamma$. Since right multiplication by $\gamma$ maps $\Delta_n$ bijectively to itself, it induces a permutation of the finite set $\Gamma\backslash\Delta_n$. Hence for each $\alpha\in R_n$ there are unique $\beta(\alpha)\in R_n$ and some $\delta_\alpha\in\Gamma$ such that
\begin{align*}
\alpha\gamma=\delta_\alpha \beta(\alpha).
\end{align*}
Using associativity of the slash action and the $\Gamma$-invariance of $f$, we obtain
\begin{align*}
(T_n f)|_k\gamma
&=
\sum_{\alpha\in R_n} f|_k(\alpha\gamma)\\
&=
\sum_{\alpha\in R_n} f|_k(\delta_\alpha\beta(\alpha))\\
&=
\sum_{\alpha\in R_n} (f|_k\delta_\alpha)|_k\beta(\alpha)\\
&=
\sum_{\alpha\in R_n} f|_k\beta(\alpha)\\
&=
\sum_{\beta\in R_n} f|_k\beta\\
&=
T_n f.
\end{align*}
Thus $T_n f$ transforms with weight $k$ under every $\gamma\in\Gamma$.
[guided]
We now prove the central invariance statement. Fix $\gamma\in\Gamma$. Since $\det\gamma=1$, right multiplication by $\gamma$ preserves determinants:
\begin{align*}
\det(\alpha\gamma)=\det(\alpha)\det(\gamma)=n.
\end{align*}
Thus $\alpha\mapsto \alpha\gamma$ is a bijection $\Delta_n\to\Delta_n$, with inverse $\alpha\mapsto \alpha\gamma^{-1}$. Consequently it permutes the left cosets in $\Gamma\backslash\Delta_n$.
Because $R_n$ contains exactly one representative from each left coset, for every $\alpha\in R_n$ there are a unique representative $\beta(\alpha)\in R_n$ and an element $\delta_\alpha\in\Gamma$ such that
\begin{align*}
\alpha\gamma=\delta_\alpha\beta(\alpha).
\end{align*}
The map $\alpha\mapsto\beta(\alpha)$ is a permutation of $R_n$.
Now compute using the slash action. First,
\begin{align*}
(T_n f)|_k\gamma
=
\left(\sum_{\alpha\in R_n} f|_k\alpha\right)|_k\gamma
=
\sum_{\alpha\in R_n} f|_k(\alpha\gamma),
\end{align*}
where the last equality uses linearity of the slash action and associativity. Substitute the coset relation $\alpha\gamma=\delta_\alpha\beta(\alpha)$:
\begin{align*}
(T_n f)|_k\gamma
=
\sum_{\alpha\in R_n} f|_k(\delta_\alpha\beta(\alpha))
=
\sum_{\alpha\in R_n} (f|_k\delta_\alpha)|_k\beta(\alpha).
\end{align*}
Since $\delta_\alpha\in\Gamma$ and $f\in M_k(\Gamma)$, the modular transformation law gives $f|_k\delta_\alpha=f$. Therefore
\begin{align*}
(T_n f)|_k\gamma
=
\sum_{\alpha\in R_n} f|_k\beta(\alpha).
\end{align*}
Because $\alpha\mapsto\beta(\alpha)$ is a permutation of $R_n$, the final sum is just the original Hecke sum:
\begin{align*}
(T_n f)|_k\gamma
=
\sum_{\beta\in R_n} f|_k\beta
=
T_n f.
\end{align*}
This proves that $T_n f$ satisfies the same weight-$k$ modular transformation law as $f$.
[/guided]
[/step]
[step:Verify holomorphy on the upper half-plane]
For each $\alpha\in R_n$, the map
\begin{align*}
\varphi_\alpha:\mathbb{H}&\to\mathbb{H}\\
z&\mapsto \frac{az+b}{cz+d}
\end{align*}
is holomorphic, and the factor $z\mapsto \det(\alpha)^{k-1}(cz+d)^{-k}$ is holomorphic on $\mathbb{H}$. Since $f:\mathbb{H}\to\mathbb{C}$ is holomorphic, each summand $f|_k\alpha$ is holomorphic on $\mathbb{H}$. The sum defining $T_n f$ is finite, so $T_n f$ is holomorphic on $\mathbb{H}$.
[/step]
[step:Use triangular representatives to compute the Fourier expansion at the cusp]
We first justify the triangular representatives. Let
\begin{align*}
\alpha=
\begin{pmatrix}
A & B\\
C & D
\end{pmatrix}
\in\Delta_n.
\end{align*}
Choose $a:=\gcd(A,C)>0$ and integers $r,s\in\mathbb{Z}$ with $rA+sC=a$. Then
\begin{align*}
\gamma_0:=
\begin{pmatrix}
r & s\\
-C/a & A/a
\end{pmatrix}
\in SL_2(\mathbb{Z}),
\end{align*}
because $\det\gamma_0=(rA+sC)/a=1$, and
\begin{align*}
\gamma_0\alpha=
\begin{pmatrix}
a & b_0\\
0 & d
\end{pmatrix}
\end{align*}
with $d=n/a>0$. Multiplying on the left by
\begin{align*}
\begin{pmatrix}
1 & t\\
0 & 1
\end{pmatrix}
\in\Gamma
\end{align*}
replaces $b_0$ by $b_0+td$, so a unique choice of $t\in\mathbb{Z}$ gives $0\le b<d$. Conversely, if two such triangular matrices lie in the same left $\Gamma$-coset, then comparing the bottom-left, determinant, and top-right entries forces the same $a$, the same $d$, and the same residue $b$ modulo $d$; the condition $0\le b<d$ gives equality. Hence the matrices
\begin{align*}
\alpha_{a,b,d}:=
\begin{pmatrix}
a & b\\
0 & d
\end{pmatrix},
\qquad
ad=n,\quad d>0,\quad 0\le b<d,
\end{align*}
are exactly one representative for each class in $\Gamma\backslash\Delta_n$. With this choice,
\begin{align*}
(T_n f)(z)
=
\sum_{\substack{ad=n\\ d>0}}
\sum_{b=0}^{d-1}
n^{k-1}d^{-k} f\left(\frac{az+b}{d}\right).
\end{align*}
Since $f$ is holomorphic at the cusp, there are coefficients $(c_m)_{m\ge 0}\subset\mathbb{C}$ such that, for $q=e^{2\pi i z}$ and $\operatorname{Im}z$ sufficiently large,
\begin{align*}
f(z)=\sum_{m=0}^{\infty} c_m q^m.
\end{align*}
For each fixed pair $(a,d)$ with $ad=n$ and $d>0$, the substitution $w=(az+b)/d$ satisfies $\operatorname{Im}w=a\operatorname{Im}z/d$. Thus, after increasing the lower bound on $\operatorname{Im}z$ if necessary, all finitely many points $(az+b)/d$ lie in the half-plane where the [Fourier series](/page/Fourier%20Series) of $f$ converges locally uniformly. Since the sums over $(a,d)$ and $b$ are finite, substitution and rearrangement with these finite sums are justified by local [uniform convergence](/page/Uniform%20Convergence) of the Fourier series on that half-plane. Substituting this expansion gives
\begin{align*}
(T_n f)(z)
&=
\sum_{\substack{ad=n\\ d>0}}
n^{k-1}d^{-k}
\sum_{b=0}^{d-1}
\sum_{m=0}^{\infty}
c_m
\exp\left(2\pi i m\frac{az+b}{d}\right)\\
&=
\sum_{\substack{ad=n\\ d>0}}
n^{k-1}d^{-k}
\sum_{m=0}^{\infty}
c_m q^{ma/d}
\sum_{b=0}^{d-1}\exp\left(\frac{2\pi i mb}{d}\right).
\end{align*}
The finite geometric sum satisfies
\begin{align*}
\sum_{b=0}^{d-1}\exp\left(\frac{2\pi i mb}{d}\right)
=
\begin{cases}
d, & d\mid m,\\
0, & d\nmid m.
\end{cases}
\end{align*}
Writing $m=d\ell$ in the surviving terms, with $\ell\in\mathbb{Z}_{\ge 0}$, we obtain
\begin{align*}
(T_n f)(z)
&=
\sum_{\substack{ad=n\\ d>0}}
n^{k-1}d^{1-k}
\sum_{\ell=0}^{\infty}
c_{d\ell}q^{a\ell}.
\end{align*}
Therefore the coefficient of $q^r$ is
\begin{align*}
\sum_{a\mid \gcd(n,r)} a^{k-1}c_{nr/a^2},
\end{align*}
with the usual interpretation that, for $r=0$, every $a\mid n$ contributes the constant coefficient $a^{k-1}c_0$. In particular, no negative powers of $q$ occur.
[guided]
The remaining condition for being a modular form is holomorphy at the cusp. We check it by using representatives adapted to Fourier expansions, and first verify why those representatives are valid. Let
\begin{align*}
\alpha=
\begin{pmatrix}
A & B\\
C & D
\end{pmatrix}
\in\Delta_n.
\end{align*}
Set $a:=\gcd(A,C)>0$. Choose $r,s\in\mathbb{Z}$ such that $rA+sC=a$, and define
\begin{align*}
\gamma_0:=
\begin{pmatrix}
r & s\\
-C/a & A/a
\end{pmatrix}.
\end{align*}
Then $\gamma_0\in SL_2(\mathbb{Z})$ because $\det\gamma_0=(rA+sC)/a=1$, and left multiplication gives
\begin{align*}
\gamma_0\alpha=
\begin{pmatrix}
a & b_0\\
0 & d
\end{pmatrix}
\end{align*}
with $d=n/a>0$. Multiplying on the left by $\begin{pmatrix}1&t\\0&1\end{pmatrix}\in\Gamma$ changes $b_0$ to $b_0+td$, so we may choose the unique residue $b$ satisfying $0\le b<d$. Uniqueness follows because if $\gamma\begin{pmatrix}a&b\\0&d\end{pmatrix}=\begin{pmatrix}a'&b'\\0&d'\end{pmatrix}$ with $\gamma=\begin{pmatrix}p&q\\r&s\end{pmatrix}\in\Gamma$, then the bottom-left entry gives $ra=0$, hence $r=0$; from $ps=1$ and $d,d'>0$ we get $p=s=1$, so $a'=a$, $d'=d$, and $b'=b+qd$. Since both $b$ and $b'$ lie in $\{0,\dots,d-1\}$, $b=b'$. Therefore a representative set for $\Gamma\backslash\Delta_n$ is
\begin{align*}
\alpha_{a,b,d}:=
\begin{pmatrix}
a & b\\
0 & d
\end{pmatrix},
\qquad
ad=n,\quad d>0,\quad 0\le b<d.
\end{align*}
For such a matrix, the slash action becomes
\begin{align*}
(f|_k\alpha_{a,b,d})(z)
=
n^{k-1}d^{-k}f\left(\frac{az+b}{d}\right),
\end{align*}
because the lower-left entry is $0$ and the lower-right entry is $d$. Hence
\begin{align*}
(T_n f)(z)
=
\sum_{\substack{ad=n\\ d>0}}
\sum_{b=0}^{d-1}
n^{k-1}d^{-k} f\left(\frac{az+b}{d}\right).
\end{align*}
Since $f$ is holomorphic at the cusp, it has a Fourier expansion with no negative powers:
\begin{align*}
f(z)=\sum_{m=0}^{\infty} c_m q^m,
\qquad q=e^{2\pi i z},
\end{align*}
for all $z\in\mathbb{H}$ with $\operatorname{Im}z$ sufficiently large. For each fixed divisor pair $(a,d)$ with $ad=n$ and $d>0$, the imaginary part of $w=(az+b)/d$ is $a\operatorname{Im}z/d$. Because there are only finitely many such pairs and finitely many values of $b$ for each $d$, we can require $\operatorname{Im}z$ to be large enough that every point $(az+b)/d$ lies in the half-plane where the Fourier expansion of $f$ converges locally uniformly. This local uniform convergence justifies substituting the expansion into each summand and then moving the finite sums over $a,d,b$ past the infinite Fourier series. Substituting $w=(az+b)/d$ gives
\begin{align*}
f\left(\frac{az+b}{d}\right)
=
\sum_{m=0}^{\infty}
c_m
\exp\left(2\pi i m\frac{az+b}{d}\right)
=
\sum_{m=0}^{\infty}
c_m q^{ma/d}\exp\left(\frac{2\pi i mb}{d}\right).
\end{align*}
Therefore
\begin{align*}
(T_n f)(z)
&=
\sum_{\substack{ad=n\\ d>0}}
n^{k-1}d^{-k}
\sum_{b=0}^{d-1}
\sum_{m=0}^{\infty}
c_m q^{ma/d}\exp\left(\frac{2\pi i mb}{d}\right)\\
&=
\sum_{\substack{ad=n\\ d>0}}
n^{k-1}d^{-k}
\sum_{m=0}^{\infty}
c_m q^{ma/d}
\sum_{b=0}^{d-1}\exp\left(\frac{2\pi i mb}{d}\right).
\end{align*}
The inner sum over $b$ is a finite geometric sum over the $d$th roots of unity. It equals $d$ exactly when $d$ divides $m$, and it equals $0$ otherwise:
\begin{align*}
\sum_{b=0}^{d-1}\exp\left(\frac{2\pi i mb}{d}\right)
=
\begin{cases}
d, & d\mid m,\\
0, & d\nmid m.
\end{cases}
\end{align*}
Thus only the terms $m=d\ell$, with $\ell\in\mathbb{Z}_{\ge 0}$, remain. Substituting $m=d\ell$ gives
\begin{align*}
(T_n f)(z)
&=
\sum_{\substack{ad=n\\ d>0}}
n^{k-1}d^{1-k}
\sum_{\ell=0}^{\infty}
c_{d\ell}q^{a\ell}.
\end{align*}
This is a Fourier series in nonnegative integral powers of $q$, because $a\ell\ge 0$. Equivalently, if $r\ge 0$, the coefficient of $q^r$ is obtained by requiring $r=a\ell$, so $a\mid r$ and $d=n/a$. The coefficient is then
\begin{align*}
\sum_{a\mid \gcd(n,r)} a^{k-1}c_{nr/a^2}.
\end{align*}
For $r=0$, the same formula means that every divisor $a$ of $n$ contributes to the constant term. The important conclusion is that no negative powers of $q$ appear, so $T_n f$ is holomorphic at the cusp.
[/guided]
[/step]
[step:Conclude that the Hecke operator preserves $M_k(\Gamma)$]
We have shown that $T_n f$ is holomorphic on $\mathbb{H}$, satisfies $(T_n f)|_k\gamma=T_n f$ for every $\gamma\in\Gamma$, and has a Fourier expansion at the cusp with no negative powers of $q$. Hence $T_n f\in M_k(\Gamma)$. Since $f\in M_k(\Gamma)$ was arbitrary, $T_n(M_k(\Gamma))\subseteq M_k(\Gamma)$.
[/step]