[proofplan]
We use the double-coset form of the full-level Hecke operator. A single determinant-$n$ slash operator is not self-adjoint on the quotient by $\Gamma$, so we first pass to the finite cover on which that single term is invariant and then change variables. This identifies its adjoint term with the adjugate matrix. Finally, the adjugate map is the involution of the determinant-$n$ Hecke correspondence that swaps the two projections. Summing over the finite correspondence gives the same Hecke operator in the second variable.
[/proofplan]
[step:Write the Hecke operator as a determinant-$n$ correspondence]
Let $M_2(\mathbb{Z})$ denote the ring of $2\times 2$ integer matrices, and let
\begin{align*}
\Delta_n
:=
\{\alpha\in M_2(\mathbb{Z}):\det\alpha=n\}.
\end{align*}
For
\begin{align*}
\alpha=\begin{pmatrix}a&b\\ c&d\end{pmatrix}\in GL_2^+(\mathbb{R}),
\end{align*}
define the weight-$k$ slash operator by
\begin{align*}
(f|_k\alpha)(z)
:=
(\det\alpha)^{k/2}(cz+d)^{-k}f(\alpha z),
\qquad
\alpha z:=\frac{az+b}{cz+d}.
\end{align*}
For the full modular group, the matrices
\begin{align*}
R_n
:=
\left\{
\begin{pmatrix}a&b\\0&d\end{pmatrix}
:
a,d\in\mathbb{N},\ ad=n,\ 0\leq b<d
\right\}
\end{align*}
form a set of representatives for the left cosets $\Gamma\backslash\Delta_n$. For $\alpha=\begin{pmatrix}a&b\\0&d\end{pmatrix}\in R_n$, one has
\begin{align*}
(f|_k\alpha)(z)
=
n^{k/2}d^{-k}f\left(\frac{az+b}{d}\right).
\end{align*}
Therefore the stated normalization is equivalently
\begin{align*}
T_n f
=
n^{k/2-1}
\sum_{\alpha\in R_n} f|_k\alpha .
\end{align*}
[guided]
The operator is best handled as a finite correspondence. We define
\begin{align*}
\Delta_n
:=
\{\alpha\in M_2(\mathbb{Z}):\det\alpha=n\}
\end{align*}
and use the slash operator
\begin{align*}
(f|_k\alpha)(z)
:=
(\det\alpha)^{k/2}(cz+d)^{-k}f(\alpha z),
\qquad
\alpha z:=\frac{az+b}{cz+d}.
\end{align*}
For $\Gamma=SL_2(\mathbb{Z})$, a standard representative set for $\Gamma\backslash\Delta_n$ is
\begin{align*}
R_n
:=
\left\{
\begin{pmatrix}a&b\\0&d\end{pmatrix}
:
a,d\in\mathbb{N},\ ad=n,\ 0\leq b<d
\right\}.
\end{align*}
For $\alpha=\begin{pmatrix}a&b\\0&d\end{pmatrix}$ in this set,
\begin{align*}
(f|_k\alpha)(z)
=
n^{k/2}d^{-k}f\left(\frac{az+b}{d}\right),
\end{align*}
so the formula in the statement becomes
\begin{align*}
T_n f
=
n^{k/2-1}
\sum_{\alpha\in R_n} f|_k\alpha .
\end{align*}
[/guided]
[/step]
[step:Compute the adjoint of one sheet on a finite cover]
Fix $\alpha\in\Delta_n$, and define its adjugate-normalized partner
\begin{align*}
\alpha^\vee:=n\alpha^{-1}.
\end{align*}
Let
\begin{align*}
\Gamma_\alpha:=\Gamma\cap\alpha^{-1}\Gamma\alpha.
\end{align*}
This subgroup has finite index in $\Gamma$. If $\mathcal{F}$ is a measurable fundamental domain for $\Gamma$ and
\begin{align*}
\Gamma=\bigsqcup_{\delta\in C_\alpha}\Gamma_\alpha\delta,
\end{align*}
then
\begin{align*}
\mathcal{F}_\alpha:=\bigcup_{\delta\in C_\alpha}\delta\mathcal{F}
\end{align*}
is a measurable fundamental domain for $\Gamma_\alpha$, up to boundary sets of $\mathcal{L}^2$-measure zero.
For $\gamma\in\Gamma_\alpha$, the relation $\alpha\gamma\alpha^{-1}\in\Gamma$ gives
\begin{align*}
(f|_k\alpha)|_k\gamma
=
f|_k(\alpha\gamma)
=
f|_k\bigl((\alpha\gamma\alpha^{-1})\alpha\bigr)
=
f|_k\alpha,
\end{align*}
and $g|_k\gamma=g$. Hence the Petersson integrand
\begin{align*}
(f|_k\alpha)(z)\overline{g(z)}\,(\operatorname{Im}z)^k\,d\mu_{\mathbb{H}}(z)
\end{align*}
is $\Gamma_\alpha$-invariant. Thus
\begin{align*}
\int_{\mathcal{F}}
(f|_k\alpha)(z)\overline{g(z)}\,(\operatorname{Im}z)^k\,d\mu_{\mathbb{H}}(z)
=
\frac{1}{[\Gamma:\Gamma_\alpha]}
\int_{\mathcal{F}_\alpha}
(f|_k\alpha)(z)\overline{g(z)}\,(\operatorname{Im}z)^k\,d\mu_{\mathbb{H}}(z).
\end{align*}
Now set $w=\alpha z$. If $\alpha=\begin{pmatrix}a&b\\ c&d\end{pmatrix}$, then
\begin{align*}
\operatorname{Im}(\alpha z)
=
\frac{n\,\operatorname{Im}z}{|cz+d|^2},
\end{align*}
and the hyperbolic measure $d\mu_{\mathbb{H}}$ is invariant under $GL_2^+(\mathbb{R})$. Also
\begin{align*}
\alpha^\vee=\begin{pmatrix}d&-b\\-c&a\end{pmatrix},
\qquad
-cw+a=\frac{n}{cz+d}.
\end{align*}
Therefore
\begin{align*}
\overline{(g|_k\alpha^\vee)(w)}
=
n^{-k/2}(c\overline{z}+d)^k\overline{g(z)}.
\end{align*}
Combining this with
\begin{align*}
(f|_k\alpha)(z)
=
n^{k/2}(cz+d)^{-k}f(w)
\end{align*}
and
\begin{align*}
(\operatorname{Im}w)^k
=
n^k|cz+d|^{-2k}(\operatorname{Im}z)^k
\end{align*}
gives
\begin{align*}
(f|_k\alpha)(z)\overline{g(z)}\,(\operatorname{Im}z)^k\,d\mu_{\mathbb{H}}(z)
=
f(w)\overline{(g|_k\alpha^\vee)(w)}\,(\operatorname{Im}w)^k\,d\mu_{\mathbb{H}}(w).
\end{align*}
The image $\alpha\mathcal{F}_\alpha$ is a fundamental domain for
\begin{align*}
\alpha\Gamma_\alpha\alpha^{-1}
=
\Gamma\cap\alpha\Gamma\alpha^{-1}.
\end{align*}
The integrand on the right is invariant under this subgroup: if $\eta=\alpha\gamma\alpha^{-1}$ with $\gamma\in\Gamma_\alpha$, then
\begin{align*}
\alpha^\vee\eta(\alpha^\vee)^{-1}
=
\gamma\in\Gamma,
\end{align*}
so $(g|_k\alpha^\vee)|_k\eta=g|_k\alpha^\vee$, and $f|_k\eta=f$. Averaging back down to the $\Gamma$-quotient yields
\begin{align*}
\int_{\mathcal{F}}
(f|_k\alpha)(z)\overline{g(z)}\,(\operatorname{Im}z)^k\,d\mu_{\mathbb{H}}(z)
=
\int_{\mathcal{F}}
f(w)\overline{(g|_k\alpha^\vee)(w)}\,(\operatorname{Im}w)^k\,d\mu_{\mathbb{H}}(w).
\end{align*}
[guided]
The single slash term is not invariant under all of $\Gamma$, so the computation must be made on the finite cover where it is invariant. Define
\begin{align*}
\alpha^\vee:=n\alpha^{-1},
\qquad
\Gamma_\alpha:=\Gamma\cap\alpha^{-1}\Gamma\alpha.
\end{align*}
Choose $C_\alpha\subset\Gamma$ with
\begin{align*}
\Gamma=\bigsqcup_{\delta\in C_\alpha}\Gamma_\alpha\delta,
\end{align*}
and put
\begin{align*}
\mathcal{F}_\alpha:=\bigcup_{\delta\in C_\alpha}\delta\mathcal{F}.
\end{align*}
For $\gamma\in\Gamma_\alpha$,
\begin{align*}
(f|_k\alpha)|_k\gamma
=
f|_k(\alpha\gamma)
=
f|_k\bigl((\alpha\gamma\alpha^{-1})\alpha\bigr)
=
f|_k\alpha,
\end{align*}
and $g|_k\gamma=g$. Hence
\begin{align*}
\int_{\mathcal{F}}
(f|_k\alpha)(z)\overline{g(z)}\,(\operatorname{Im}z)^k\,d\mu_{\mathbb{H}}(z)
=
\frac{1}{[\Gamma:\Gamma_\alpha]}
\int_{\mathcal{F}_\alpha}
(f|_k\alpha)(z)\overline{g(z)}\,(\operatorname{Im}z)^k\,d\mu_{\mathbb{H}}(z).
\end{align*}
Under $w=\alpha z$ one has
\begin{align*}
\operatorname{Im}(\alpha z)
=
\frac{n\,\operatorname{Im}z}{|cz+d|^2},
\qquad
d\mu_{\mathbb{H}}(z)=d\mu_{\mathbb{H}}(w),
\end{align*}
and, for $\alpha^\vee=\begin{pmatrix}d&-b\\-c&a\end{pmatrix}$,
\begin{align*}
-cw+a=\frac{n}{cz+d}.
\end{align*}
Thus
\begin{align*}
\overline{(g|_k\alpha^\vee)(w)}
=
n^{-k/2}(c\overline{z}+d)^k\overline{g(z)}.
\end{align*}
Together with
\begin{align*}
(f|_k\alpha)(z)
=
n^{k/2}(cz+d)^{-k}f(w),
\qquad
(\operatorname{Im}w)^k
=
n^k|cz+d|^{-2k}(\operatorname{Im}z)^k,
\end{align*}
this gives the transformed integrand
\begin{align*}
(f|_k\alpha)(z)\overline{g(z)}\,(\operatorname{Im}z)^k\,d\mu_{\mathbb{H}}(z)
=
f(w)\overline{(g|_k\alpha^\vee)(w)}\,(\operatorname{Im}w)^k\,d\mu_{\mathbb{H}}(w).
\end{align*}
The image $\alpha\mathcal{F}_\alpha$ is a fundamental domain for
\begin{align*}
\alpha\Gamma_\alpha\alpha^{-1}
=
\Gamma\cap\alpha\Gamma\alpha^{-1}.
\end{align*}
If $\eta=\alpha\gamma\alpha^{-1}$, then
\begin{align*}
\alpha^\vee\eta(\alpha^\vee)^{-1}
=
\gamma,
\end{align*}
so $g|_k\alpha^\vee$ is invariant under this image subgroup in the same way that $f|_k\alpha$ was invariant under $\Gamma_\alpha$. Averaging the image fundamental domain back down to $\mathcal{F}$ gives
\begin{align*}
\int_{\mathcal{F}}
(f|_k\alpha)(z)\overline{g(z)}\,(\operatorname{Im}z)^k\,d\mu_{\mathbb{H}}(z)
=
\int_{\mathcal{F}}
f(w)\overline{(g|_k\alpha^\vee)(w)}\,(\operatorname{Im}w)^k\,d\mu_{\mathbb{H}}(w).
\end{align*}
[/guided]
[/step]
[step:Sum over the symmetric determinant-$n$ correspondence]
Summing the one-sheet identity over $\alpha\in R_n$ and using the double-coset expression for $T_n$ gives
\begin{align*}
(T_n f,g)_{\mathrm{Pet}}
=
n^{k/2-1}
\sum_{\alpha\in R_n}
\int_{\mathcal{F}}
f(w)\overline{(g|_k\alpha^\vee)(w)}\,(\operatorname{Im}w)^k\,d\mu_{\mathbb{H}}(w).
\end{align*}
It remains only to identify the finite sum in the second variable with $T_n g$.
Define the determinant-$n$ Hecke correspondence as the quotient
\begin{align*}
\mathcal{C}_n
:=
(\mathbb{H}\times\Delta_n)/(\Gamma\times\Gamma),
\end{align*}
where
\begin{align*}
(\gamma_1,\gamma_2)\cdot(z,\alpha)
:=
(\gamma_2^{-1}z,\gamma_1\alpha\gamma_2).
\end{align*}
The two maps
\begin{align*}
p_1([z,\alpha])&:=[z],\\
p_2([z,\alpha])&:=[\alpha z]
\end{align*}
from $\mathcal{C}_n$ to $\Gamma\backslash\mathbb{H}$ are well defined. In the sheet decomposition determined by the left-coset representatives $R_n$, the operator $T_n$ is the pull-push operator $n^{k/2-1}p_{1*}p_2^*$, which is precisely the sum
\begin{align*}
n^{k/2-1}\sum_{\alpha\in R_n}g|_k\alpha.
\end{align*}
The map
\begin{align*}
\tau:\mathcal{C}_n&\to\mathcal{C}_n\\
[z,\alpha]&\mapsto[\alpha z,\alpha^\vee]
\end{align*}
is well defined. Indeed, for $\gamma_1,\gamma_2\in\Gamma$,
\begin{align*}
(\gamma_1\alpha\gamma_2)^\vee
=
\gamma_2^{-1}\alpha^\vee\gamma_1^{-1},
\end{align*}
which is exactly the [equivalence relation](/page/Equivalence%20Relation) in $\mathcal{C}_n$. Also
\begin{align*}
(\alpha^\vee)^\vee=\alpha,
\end{align*}
so $\tau$ is an involution. Finally,
\begin{align*}
p_1\circ\tau=p_2,
\qquad
p_2\circ\tau=p_1.
\end{align*}
Thus the same correspondence is obtained after swapping the two projections. Consequently the adjoint finite sum
\begin{align*}
n^{k/2-1}\sum_{\alpha\in R_n}g|_k\alpha^\vee
\end{align*}
is the same pull-push operator as
\begin{align*}
n^{k/2-1}\sum_{\alpha\in R_n}g|_k\alpha
=
T_n g.
\end{align*}
This is an equality of the correspondence operator, not a termwise replacement of right cosets by left cosets.
Substituting this identification into the previous displayed formula gives
\begin{align*}
(T_n f,g)_{\mathrm{Pet}}
&=
\int_{\mathcal{F}}
f(w)\overline{(T_n g)(w)}\,(\operatorname{Im}w)^k\,d\mu_{\mathbb{H}}(w)\\
&=
(f,T_n g)_{\mathrm{Pet}}.
\end{align*}
The cusp condition gives exponential decay at the cusp, so all integrals above are absolutely convergent and all finite sums may be interchanged with integration.
[guided]
After the one-sheet calculation, summing over $R_n$ gives
\begin{align*}
(T_n f,g)_{\mathrm{Pet}}
=
n^{k/2-1}
\sum_{\alpha\in R_n}
\int_{\mathcal{F}}
f(w)\overline{(g|_k\alpha^\vee)(w)}\,(\operatorname{Im}w)^k\,d\mu_{\mathbb{H}}(w).
\end{align*}
The point is not that each right-equivalent matrix can be replaced term by term. Instead, the whole finite sum is the trace over a symmetric correspondence.
The correspondence is
\begin{align*}
\mathcal{C}_n
:=
(\mathbb{H}\times\Delta_n)/(\Gamma\times\Gamma),
\end{align*}
with action
\begin{align*}
(\gamma_1,\gamma_2)\cdot(z,\alpha)
:=
(\gamma_2^{-1}z,\gamma_1\alpha\gamma_2).
\end{align*}
It has projections
\begin{align*}
p_1([z,\alpha])&:=[z],\\
p_2([z,\alpha])&:=[\alpha z].
\end{align*}
Using the left-coset representatives $R_n$, the pull-push operator $n^{k/2-1}p_{1*}p_2^*$ is the usual Hecke sum
\begin{align*}
n^{k/2-1}\sum_{\alpha\in R_n}g|_k\alpha.
\end{align*}
Now define
\begin{align*}
\tau:\mathcal{C}_n&\to\mathcal{C}_n\\
[z,\alpha]&\mapsto[\alpha z,\alpha^\vee].
\end{align*}
The identity
\begin{align*}
(\gamma_1\alpha\gamma_2)^\vee
=
\gamma_2^{-1}\alpha^\vee\gamma_1^{-1}
\end{align*}
shows that $\tau$ is well defined on the quotient, and
\begin{align*}
(\alpha^\vee)^\vee=\alpha
\end{align*}
shows that $\tau$ is an involution. Moreover
\begin{align*}
p_1\circ\tau=p_2,
\qquad
p_2\circ\tau=p_1.
\end{align*}
Thus the finite correspondence is symmetric. Therefore the finite sum with $\alpha^\vee$ in the second slot is the same Hecke pull-push operator as the original sum:
\begin{align*}
n^{k/2-1}\sum_{\alpha\in R_n}g|_k\alpha^\vee
=
T_n g.
\end{align*}
Substituting this into the summed adjoint identity gives
\begin{align*}
(T_n f,g)_{\mathrm{Pet}}
&=
\int_{\mathcal{F}}
f(w)\overline{(T_n g)(w)}\,(\operatorname{Im}w)^k\,d\mu_{\mathbb{H}}(w)\\
&=
(f,T_n g)_{\mathrm{Pet}}.
\end{align*}
This proves self-adjointness.
[/guided]
[/step]
