[proofplan]
We write the Hecke operator using the weight-$k$ slash action and the finite double-coset set of integral matrices of determinant $n$. The key point is that right multiplication by an element of $SL_2(\mathbb{Z})$ permutes the left cosets in this determinant-$n$ set, so the Hecke sum remains modular of weight $k$. Holomorphy follows because the operator is a finite sum of holomorphic pullbacks. Finally, the explicit Fourier expansion formula shows that the constant term of $T_n f$ is a multiple of the constant term of $f$, hence cusp forms are preserved.
[/proofplan]
[step:Represent the Hecke operator by determinant-$n$ double cosets]
Let $\mathbb{H} = \{z \in \mathbb{C} : \operatorname{Im}(z) > 0\}$ denote the complex upper half-plane, and let
\begin{align*}
\Delta_n := \left\{ A \in M_2(\mathbb{Z}) : \det A = n \right\}.
\end{align*}
For a matrix $A = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \in GL_2^+(\mathbb{R})$, define the weight-$k$ slash action on a [holomorphic function](/page/Holomorphic%20Function) $f: \mathbb{H} \to \mathbb{C}$ by
\begin{align*}
(f|_k A)(z) := \det(A)^{k/2}(cz+d)^{-k} f\left(\frac{az+b}{cz+d}\right).
\end{align*}
The standard level-one Hecke operator $T_n$ is given by
\begin{align*}
T_n f := n^{k/2-1} \sum_{A \in \Gamma \backslash \Delta_n} f|_k A,
\end{align*}
where one representative is chosen from each left $\Gamma$-coset in $\Delta_n$.
This definition is independent of the chosen representatives: if $A'=\gamma A$ with $\gamma \in \Gamma$, then
\begin{align*}
f|_k A' = f|_k(\gamma A) = (f|_k \gamma)|_k A = f|_k A,
\end{align*}
because $f \in M_k(\Gamma)$ satisfies $f|_k\gamma=f$ for every $\gamma \in \Gamma$.
[/step]
[step:Show that the Hecke sum is invariant under $SL_2(\mathbb{Z})$]
Let $\gamma \in \Gamma$. Since $\det \gamma=1$, right multiplication by $\gamma$ defines a bijection
\begin{align*}
R_\gamma: \Gamma \backslash \Delta_n &\to \Gamma \backslash \Delta_n,\\
\Gamma A &\mapsto \Gamma A\gamma.
\end{align*}
It is well-defined because if $\Gamma A=\Gamma A'$, then $A'=\delta A$ for some $\delta \in \Gamma$, hence $A'\gamma=\delta A\gamma$ and $\Gamma A'\gamma=\Gamma A\gamma$. It is bijective with inverse $R_{\gamma^{-1}}$.
Now let $f \in M_k(\Gamma)$. Using associativity of the slash action and the preceding bijection,
\begin{align*}
(T_n f)|_k\gamma
&= n^{k/2-1}\sum_{A \in \Gamma\backslash\Delta_n} (f|_k A)|_k\gamma\\
&= n^{k/2-1}\sum_{A \in \Gamma\backslash\Delta_n} f|_k(A\gamma)\\
&= n^{k/2-1}\sum_{B \in \Gamma\backslash\Delta_n} f|_k B\\
&= T_n f.
\end{align*}
Thus $T_n f$ transforms with weight $k$ under every element of $\Gamma$.
[/step]
[step:Verify holomorphy on the upper half-plane and at the cusp]
For each $A=\begin{pmatrix}a&b\\c&d\end{pmatrix}\in\Delta_n$, the map
\begin{align*}
\phi_A: \mathbb{H} &\to \mathbb{H},\\
z &\mapsto \frac{az+b}{cz+d}
\end{align*}
is holomorphic, and $cz+d\neq 0$ for $z\in\mathbb{H}$. Since $f$ is holomorphic on $\mathbb{H}$, each summand $f|_kA$ is holomorphic on $\mathbb{H}$. The quotient $\Gamma\backslash\Delta_n$ is finite, so $T_n f$ is holomorphic on $\mathbb{H}$.
Choose the standard representatives
\begin{align*}
A_{a,b,d}:=\begin{pmatrix}a&b\\0&d\end{pmatrix},
\qquad ad=n,\qquad 0\leq b<d.
\end{align*}
Using these representatives, the Hecke operator is
\begin{align*}
(T_n f)(z)
= n^{k-1}\sum_{\substack{ad=n\\d>0}} d^{-k}\sum_{b=0}^{d-1} f\left(\frac{az+b}{d}\right).
\end{align*}
Write the Fourier expansion of $f$ at infinity as
\begin{align*}
f(z)=\sum_{m=0}^{\infty} a_m q^m,
\qquad q:=e^{2\pi i z}.
\end{align*}
Then
\begin{align*}
f\left(\frac{az+b}{d}\right)
=
\sum_{m=0}^{\infty} a_m
\exp\left(2\pi i m\frac{az+b}{d}\right).
\end{align*}
Averaging over $b$ gives
\begin{align*}
\sum_{b=0}^{d-1}\exp\left(2\pi i m\frac{b}{d}\right)
=
\begin{cases}
d, & d\mid m,\\
0, & d\nmid m.
\end{cases}
\end{align*}
Therefore $T_n f$ has a Fourier expansion in nonnegative powers of $q$, so it is holomorphic at the cusp infinity. Hence $T_n f\in M_k(\Gamma)$.
[/step]
[step:Compute the constant term and preserve cuspidality]
The constant term of $T_n f$ comes only from the $m=0$ term in the Fourier expansion of $f$. From the explicit formula above,
\begin{align*}
a_0(T_n f)
&=
n^{k-1}\sum_{\substack{ad=n\\d>0}} d^{-k}\sum_{b=0}^{d-1} a_0\\
&=
n^{k-1}\left(\sum_{\substack{ad=n\\d>0}} d^{1-k}\right)a_0.
\end{align*}
Thus if $f\in S_k(\Gamma)$, then its constant term satisfies $a_0=0$, and consequently
\begin{align*}
a_0(T_n f)=0.
\end{align*}
For level one, $\Gamma=SL_2(\mathbb{Z})$ has a single cusp represented by infinity, so vanishing of the constant term at infinity is exactly the cusp condition. Hence $T_n f\in S_k(\Gamma)$.
Combining the modularity, holomorphy, and cusp calculations proves that $T_n$ maps $M_k(\Gamma)$ to itself and maps $S_k(\Gamma)$ to itself.
[/step]
