Solution (source code)

= Solution

Fix an even weight $k$ and write $\Gamma=SL_2(\mathbb Z)$. Extend the <slash operator for modular forms> to positive-determinant rational <matrices> by the <determinant-normalized slash operator>
$$
(f|_k\alpha)(\tau)=(\det\alpha)^{k/2}(c\tau+d)^{-k}f(\alpha\tau).
$$
The <automorphy factor> identity makes this a right action. Let $\mathcal D_n$ be the integral <matrices> of determinant $n>0$. Their left $\Gamma$-cosets have unique <determinant-n matrix representatives for Hecke operators>
$$
\alpha_{a,b,d}=\begin{pmatrix}a&b\\0&d\end{pmatrix},\qquad ad=n,\quad a,d>0,\quad0\leq b<d.
$$
An integral row operation obtained from <Bézout's identity> makes the first column $(a,0)^T$, and a row shear reduces the other top entry modulo $d$, proving this description. The usual <Hecke operator> is
$$
\boxed{T_nf=n^{k/2-1}\sum_{\Gamma\backslash\mathcal D_n}f|_k\alpha=n^{k-1}\sum_{ad=n}d^{-k}\sum_{b=0}^{d-1}f\!\left(\frac{a\tau+b}{d}\right).}
$$
This normalization has $T_1=1$. Right multiplication by $\Gamma$ permutes the left cosets in $\mathcal D_n$, so $T_nf$ has the same modular transformation law as $f$. Each summand is holomorphic in $\mathbb H$. Holomorphy at the cusp and preservation of <cusp forms> follow directly from the coefficient formula below; the full <modular group> has only one cusp orbit. Thus the operators act on both $M_k$ and $S_k$.

If $f=\sum_{r\geq0}a_rq^r$, summing over $b$ kills terms unless $d\mid r$. Writing $r=d\ell$ gives output exponent $a\ell$, so a term of exponent $m$ requires $a\mid m$ and has original index $r=dm/a=nm/a^2$. Its multiplier is $n^{k-1}d^{1-k}=a^{k-1}$. The <Fourier coefficients of a composite-index Hecke operator> are consequently
$$
\boxed{a_m(T_nf)=\sum_{a\mid\gcd(n,m)}a^{k-1}a_{nm/a^2}(f).}
$$
For $m=0$, this means $a_0(T_nf)=\sigma_{k-1}(n)a_0(f)$. In particular $a_1(T_nf)=a_n(f)$. For a <prime number> $p$,
$$
T_pf=p^{k-1}f(p\tau)+\frac1p\sum_{b=0}^{p-1}f\!\left(\frac{\tau+b}{p}\right),\qquad a_m(T_pf)=a_{pm}+p^{k-1}a_{m/p},
$$
with a nonexistent fractional-index coefficient interpreted as zero.

The <Hecke multiplication relations> follow from these coefficients. For coprime indices, the two divisor choices combine uniquely and give $T_mT_n=T_{mn}$. For prime powers put $A=p^{k-1}$, and, on formal <Fourier series>, let $U_p(\sum a_mq^m)=\sum a_{pm}q^m$ and $V_p(\sum a_mq^m)=\sum a_mq^{pm}$. These individual maps need not preserve level one, but $U_pV_p=1$, and the coefficient formula gives
$$
T_{p^r}=\sum_{j=0}^rA^jV_p^jU_p^{r-j}.
$$
Multiplying on the left by $T_p=U_p+AV_p$ and cancelling $U_pV_p$ yields
$$
T_pT_{p^r}=T_{p^{r+1}}+p^{k-1}T_{p^{r-1}}\qquad(r\geq1).
$$
Induction gives $T_{p^u}T_{p^v}=\sum_{j=0}^{\min(u,v)}p^{j(k-1)}T_{p^{u+v-2j}}$; combining distinct prime factors proves
$$
\boxed{T_mT_n=\sum_{d\mid\gcd(m,n)}d^{k-1}T_{mn/d^2}.}
$$
The <modular Hecke algebra> is therefore commutative and generated by its prime-index operators. The identity $U_pV_p=1$, rather than the generally false reversed identity, is what proves the prime-power recurrence.

On $S_k$ use the positive definite <Petersson inner product>
$$
\langle f,g\rangle=\int_{\Gamma\backslash\mathbb H}f(\tau)\overline{g(\tau)}y^k\frac{dx\,dy}{y^2}.
$$
Exponential cusp decay makes this integral finite. The <Hecke operators are self-adjoint for the Petersson inner product>. Here is the correspondence argument. For a positive-determinant <matrix> $\alpha$, set $H=\Gamma\cap\alpha^{-1}\Gamma\alpha$. The finite sum over the double coset $\Gamma\alpha\Gamma$ unfolds its pairing into
$$
\int_{H\backslash\mathbb H}(f|_k\alpha)(\tau)\overline{g(\tau)}y^k\,d\mu(\tau).
$$
Changing variables $z=\alpha\tau$ preserves hyperbolic area, and $\operatorname{Im}(\alpha\tau)=\det(\alpha)y/|c\tau+d|^2$ shows that the determinant-normalized slash factors transfer $\alpha$ from $f$ to $\alpha^{-1}$ on $g$. The new quotient is by $\alpha H\alpha^{-1}=\Gamma\cap\alpha\Gamma\alpha^{-1}$; folding it back gives the operator for the reversed double coset. For $\alpha\in\mathcal D_n$, its adjugate $n\alpha^{-1}$ again lies in $\mathcal D_n$. Positive scalar <matrices> act trivially in the determinant normalization, so $\alpha^{-1}$ and its adjugate give the same action. The adjugate involution preserves the entire determinant-$n$ correspondence, and the external normalization $n^{k/2-1}$ is real. Thus
$$
\langle T_nf,g\rangle=\langle f,T_ng\rangle.
$$

Finite-dimensional spectral theory now gives a simultaneous orthonormal <eigenbasis> on $S_k$, since the operators commute and are <self-adjoint>. For a common <Hecke eigenform> $f$ with <eigenvalues> $\lambda_n$, the first-coefficient formula gives $a_n(f)=\lambda_na_1(f)$. If $a_1(f)=0$, every coefficient vanishes and $f=0$, a contradiction. Rescale each eigenform to have $a_1=1$. The resulting basis is orthogonal, and
$$
\boxed{T_nf=a_n(f)f.}
$$
All these <eigenvalues> and coefficients are real. A fixed simultaneous <eigenspace> is one-dimensional, because its normalized <Fourier coefficients> are uniquely determined by the <eigenvalues>.

Applying the multiplication relations to an eigenform gives
$$
a_ma_n=\sum_{d\mid\gcd(m,n)}d^{k-1}a_{mn/d^2},\qquad a_{mn}=a_ma_n\ \text{if }\gcd(m,n)=1,
$$
$$
a_{p^{r+1}}=a_pa_{p^r}-p^{k-1}a_{p^{r-1}}.
$$
The recurrence yields $\sum_{r\geq0}a_{p^r}X^r=(1-a_pX+p^{k-1}X^2)^{-1}$. Multiplicativity then gives the <Euler product of a Hecke eigenform>
$$
\boxed{L(f,s)=\prod_p\left(1-a_pp^{-s}+p^{k-1-2s}\right)^{-1},}
$$
initially in a right half-plane, for example $\operatorname{Re}s>k/2+1$ using the <Fourier coefficient bound for a cusp form>. Modular inversion and the <Mellin transform of a cusp-form L-function> continue $\Lambda(f,s)=(2\pi)^{-s}\Gamma(s)L(f,s)$ to an <entire function> and give $\Lambda(f,s)=i^k\Lambda(f,k-s)$.

There is also an integral structure. If $d=\dim S_k$, choose $F_j=\Delta^jR_{k-12j}$, $1\leq j\leq d$, with $R_w$ a monomial in $E_4,E_6$ of weight $w$ and constant term one. These forms have integral <Fourier coefficients> and start with $q^j$; their first $d$ coefficients form a unitriangular <matrix>. They are a basis, and every <cusp form> with integral coefficients is an integral combination of them, by successive elimination of its first $d$ coefficients. The coefficient formula shows that $T_n$ preserves this lattice. Its characteristic polynomial is consequently monic with integer coefficients, proving that the <Hecke eigenvalues are algebraic integers>. Together with self-adjointness, this gives real algebraic integers without requiring any unproved bound for their sizes.

For even $k\geq4$, the full space decomposes as $M_k=\mathbb CE_k\oplus S_k$. The normalized <Eisenstein series> is a simultaneous eigenform with
$$
T_nE_k=\sigma_{k-1}(n)E_k.
$$
Its coefficients are $-2k\,\sigma_{k-1}(m)/B_k$, and the elementary divisor identity $\sigma_{k-1}(m)\sigma_{k-1}(n)=\sum_{d\mid\gcd(m,n)}d^{k-1}\sigma_{k-1}(mn/d^2)$ verifies the eigenvalue formula, including the constant coefficient. Thus the cusp eigenbasis plus $E_k$ diagonalizes the operators on $M_k$. Weight two and odd weights have no forms; in weight zero the constant function has eigenvalue $\sigma_{-1}(n)$. As a concrete cuspidal example, $S_{12}=\mathbb C\Delta$, so $T_n\Delta=\tau(n)\Delta$ for its <Fourier coefficients>, and the recurrence gives $\tau(p^2)=\tau(p)^2-p^{11}$.