= Solution
Write $\Gamma=SL_2(\mathbb Z)$ and $G=GL_2(\mathbb Q)^+$. The <double-coset Hecke algebra> consists of $\Gamma$-bi-invariant complex <functions> on $G$ supported on finitely many <double cosets>, with convolution
$$
(h_1*h_2)(g)=\sum_{\Gamma x\in\Gamma\backslash G}h_1(gx^{-1})h_2(x).
$$
Each <double coset> has finitely many orbits under left multiplication by $\Gamma$, by <rational conjugation of finite-index modular subgroups>. Thus the sum is finite and independent of representatives. The characteristic functions of the <double cosets> form its basis, and its right action on invariant <modular forms> is $f*h=\sum_{\Gamma x}h(x)f|_kx$, using the <determinant-normalized slash operator>.
For a positive integer $n$, let $\mathcal M_n$ be all integral two-by-two <matrices> of positive <determinant> $n$. It is $\Gamma$-bi-invariant. Define $T(n)$ to be its <indicator function>, equivalently the sum of the distinct <double cosets> it contains, each with coefficient one. This is the convention consistent with the requested formula. For composite $n$, it need not be the single <double coset> of $\operatorname{diag}(1,n)$: for instance $2I\in\mathcal M_4$ cannot lie in that <double coset>, since multiplying by unimodular integral <matrices> preserves the greatest common divisor of the entries.
Prove all the needed subgroup facts directly. For a <finite-index subgroup> $L\subseteq\mathbb Z^2$, take $a$ to be the least positive first coordinate appearing in $L$ and $d$ the least positive second coordinate on its intersection with the second axis. Euclidean division then shows that the first-coordinate projection is $a\mathbb Z$ and $L\cap(\{0\}\times\mathbb Z)=\{0\}\times d\mathbb Z$. Positivity follows, for example, because the finite <quotient group> kills a nonzero multiple of each coordinate vector. Choose $(a,b)\in L$ and reduce $b$ modulo $d$ to $0\leq b<d$. Every vector of $L$ has first coordinate a multiple of $a$, and subtracting that multiple of $(a,b)$ leaves a multiple of $(0,d)$. Therefore these two vectors form a <basis> of $L$. The parameters $a,d,b$ are unique. Reducing the first coordinate modulo $a$ and then the second modulo $d$ gives precisely $ad$ quotient representatives, so $[\mathbb Z^2:L]=ad$.
Apply this <row Hermite normal form in rank two> to the <row lattice> of an integral <matrix> $A$ with $\det A=n$. Its <row lattice> contains $n\mathbb Z^2$, since $\operatorname{adj}(A)A=nI$, so it has finite index. Its two rows and the displayed two rows are <bases> of the same <row lattice>. The two inverse change-of-basis <matrices> have integer entries, so their <determinants> are integers whose product is one. Thus the change-of-basis <matrix> has <determinant> $\pm1$; since both orientations are positive, its <determinant> is one. Thus each orbit under left multiplication by $\Gamma$ has exactly one of the <determinant-n matrix representatives for Hecke operators>
$$
\begin{pmatrix}a&b\\0&d\end{pmatrix},\qquad a,d>0,\quad ad=n,\quad0\leq b<d.
$$
There are $\sum_{d\mid n}d$ such representatives, proving finiteness as well as the formula. The subgroup argument is a proof of the relevant <Hermite normal form>, not an invocation of an unproved lattice classification.
Consequently the normalized <Hecke operator> is
$$
\boxed{T_nf=n^{k/2-1}\sum_{ad=n}\sum_{b=0}^{d-1}
f|_k\begin{pmatrix}a&b\\0&d\end{pmatrix}.}
$$
It preserves $M_k(\Gamma)$: right multiplication by $\Gamma$ permutes the left-multiplication orbits in $\mathcal M_n$, giving invariance, and <cusp holomorphy under rational slash operators> gives the holomorphy of each term at every cusp.
For a triangular representative the <determinant-normalized slash operator> is $f|_kA=n^{k/2}d^{-k}f((az+b)/d)$. Summing the <Fourier expansion of a modular form> over $b$ kills every index not divisible by $d$, by finite exponential orthogonality. Thus
$$
T_nf=n^{k-1}\sum_{ad=n}d^{1-k}\sum_{\ell\geq0}a_{d\ell}(f)q^{a\ell}.
$$
The <Fourier coefficients of a composite-index Hecke operator> are therefore
$$
\boxed{a_r(T_nf)=\sum_{a\mid\gcd(n,r)}a^{k-1}a_{nr/a^2}(f)\quad(r\geq1),
\qquad a_0(T_nf)=\sigma_{k-1}(n)a_0(f).}
$$
In particular $a_1(T_nf)=a_n(f)$. If $T_nf=\lambda f$, comparison of the $q$ coefficients gives
$$
\boxed{a_n(f)=\lambda a_1(f).}
$$
Finally suppose $a_0(f)\ne0$ and $f$ is a simultaneous eigenfunction. Constant coefficients force $\lambda_n=\sigma_{k-1}(n)$, so $a_n(f)=a_1(f)\sigma_{k-1}(n)$ for all $n\geq1$. Weight two cannot occur, by <vanishing of weight-two level-one modular forms>. For even $k\geq4$, use the normalized <Eisenstein series> with its <Fourier expansion of a normalized Eisenstein series>
$$
E_k=1+c_k\sum_{n\geq1}\sigma_{k-1}(n)q^n,
\qquad c_k=-\frac{2k}{B_k}\ne0,
$$
where $B_k$ is the <Bernoulli number>. Then $h=f-a_0(f)E_k$ is a <cusp form> with coefficients $(a_1(f)-a_0(f)c_k)\sigma_{k-1}(n)$. The <Fourier coefficient bound for a cusp form> bounds these by $O(n^{k/2})$, but at arbitrarily large <primes> $p$ their magnitude is $|a_1-a_0c_k|(1+p^{k-1})$. Since $k-1>k/2$, the coefficient factor must vanish. All coefficients of $h$ then vanish, including its constant coefficient, so $h=0$ by its cusp expansion and the <identity theorem>. This proves the <noncuspidal level-one Hecke eigenform> characterization:
$$
\boxed{f=a_0(f)E_k\qquad(k\geq4).}
$$
Back to article page