[proofplan] We prove the identity by comparing Fourier coefficients. The normalized Hecke operator has an explicit action on $q$-expansions, and this coefficient formula reduces the theorem to a finite arithmetic identity. We first verify the prime-power relation, then use the factorisation of the coefficient formula over distinct primes to assemble the general formula. [/proofplan]

[step:Write the normalized Hecke action on Fourier coefficients] Let $V$ be either $M_k(\operatorname{SL}_2(\mathbb{Z}))$ or $S_k(\operatorname{SL}_2(\mathbb{Z}))$. Let \begin{align*} \mathbb{H} &:= \{z\in\mathbb{C}:\operatorname{Im} z>0\},\\ f: \mathbb{H} &\to \mathbb{C} \end{align*} be an element of $V$, and write its Fourier expansion at infinity as \begin{align*} f(z) = \sum_{r=0}^{\infty} a_r q^r, \qquad q = e^{2\pi i z}. \end{align*} For $N \in \mathbb{N}$, we use the standard normalization of the Hecke operator $T_N: V \to V$, namely the defining Fourier-coefficient formula \begin{align*} (T_N f)(z) = \sum_{r=0}^{\infty} \left( \sum_{e \mid \gcd(N,r)} e^{k-1} a_{Nr/e^2} \right) q^r. \end{align*} For $r=0$, the inner sum is understood as the finite sum over all divisors $e \mid N$. Thus it is enough to prove the corresponding identity for the coefficient operators \begin{align*} L_N: \mathbb{C}^{\mathbb{N}_0} &\to \mathbb{C}^{\mathbb{N}_0}, \\ (a_r)_{r \geq 0} &\mapsto (b_r)_{r \geq 0}, \end{align*} where \begin{align*} b_r = \sum_{e \mid \gcd(N,r)} e^{k-1} a_{Nr/e^2}. \end{align*} Indeed, if the Fourier coefficients of two modular forms agree for every $r \geq 0$, then their difference has identically zero Fourier expansion near the cusp and therefore vanishes as a [holomorphic function](/page/Holomorphic%20Function) on $\mathbb{H}$. [/step]

[step:Prove the prime-power relation on coefficients]Fix a prime number $p$. Let $A,B \in \mathbb{N}_0$, and set \begin{align*} K := k-1. \end{align*} We prove that, for every sequence $(a_r)_{r \geq 0}$, \begin{align*} L_{p^A}L_{p^B} = \sum_{j=0}^{\min(A,B)} p^{jK} L_{p^{A+B-2j}}. \end{align*} Let $r \in \mathbb{N}$. Write \begin{align*} r = p^R r_0, \end{align*} where $R \in \mathbb{N}_0$ and $r_0 \in \mathbb{N}$ satisfies $p \nmid r_0$. For $C \in \mathbb{N}_0$, the coefficient formula gives \begin{align*} (L_{p^C}a)_r = \sum_{i=0}^{\min(C,R)} p^{iK} a_{p^{C+R-2i}r_0}. \end{align*} Applying this first with $C=B$ and then with $C=A$, we obtain \begin{align*} (L_{p^A}L_{p^B}a)_r = \sum_{i=0}^{\min(A,R)} \sum_{j=0}^{\min(B,A+R-2i)} p^{(i+j)K} a_{p^{A+B+R-2(i+j)}r_0}. \end{align*} For a fixed integer $\ell \geq 0$, the coefficient of \begin{align*} p^{\ell K}a_{p^{A+B+R-2\ell}r_0} \end{align*} on the left is the number of pairs $(i,j) \in \mathbb{N}_0^2$ satisfying \begin{align*} i+j &= \ell,\\ 0 \leq i &\leq \min(A,R),\\ 0 \leq j &\leq \min(B,A+R-2i). \end{align*} Equivalently, \begin{align*} \max(0,\ell-B) \leq i \leq \min(\ell,A,R,A+R-\ell). \end{align*} On the other hand, the coefficient of the same term on \begin{align*} \sum_{h=0}^{\min(A,B)} p^{hK} L_{p^{A+B-2h}}a \end{align*} is the number of pairs $(h,t) \in \mathbb{N}_0^2$ satisfying \begin{align*} h+t &= \ell,\\ 0 \leq h &\leq \min(A,B),\\ 0 \leq t &\leq \min(A+B-2h,R). \end{align*} Equivalently, \begin{align*} \max(0,\ell-R) \leq h \leq \min(\ell,A,B,A+B-\ell). \end{align*} The two displayed intervals have the same number of integers. Indeed, writing the number of integers in $[u,v]\cap\mathbb{Z}$ as $\max(0,v-u+1)$, this is the elementary min-max identity \begin{align*} \max\bigl(0,\min(\ell,A,R,A+R-\ell)-\max(0,\ell-B)+1\bigr) \\ = \max\bigl(0,\min(\ell,A,B,A+B-\ell)-\max(0,\ell-R)+1\bigr). \end{align*} To verify the identity, split into the four cases determined by $\ell\leq B$ or $\ell>B$, and by $\ell\leq R$ or $\ell>R$. In each case subtract the positive lower endpoint, if there is one, from the corresponding upper endpoint. The remaining upper bound is the same on both sides: \begin{align*} \min(\ell,A,\min(B,R),A+\min(B,R)-\ell), \end{align*} with the convention that a negative displayed upper bound gives the empty interval. Hence the two intervals have the same cardinality. Therefore the two expressions have the same coefficient of every sequence entry $a_s$ at every positive index $r \in \mathbb{N}$. It remains to check the constant coefficient $r=0$. Define the finite geometric sum \begin{align*} \sigma_C := \sum_{i=0}^{C} p^{iK} \end{align*} for each $C \in \mathbb{N}_0$. Since every divisor of $p^C$ divides $\gcd(p^C,0)$, the coefficient formula gives \begin{align*} (L_{p^C}a)_0 = \sigma_C a_0. \end{align*} Thus the constant coefficient of $L_{p^A}L_{p^B}a$ is $\sigma_A\sigma_B a_0$. The constant coefficient of \begin{align*} \sum_{h=0}^{\min(A,B)} p^{hK} L_{p^{A+B-2h}}a \end{align*} is \begin{align*} \left(\sum_{h=0}^{\min(A,B)} p^{hK}\sigma_{A+B-2h}\right)a_0. \end{align*} Let $x := p^K$. If $A \leq B$, the finite geometric-series identity gives \begin{align*} \sum_{h=0}^{A} x^h\sum_{t=0}^{A+B-2h}x^t &= \frac{1}{1-x}\left(\sum_{h=0}^{A}x^h-\sum_{h=0}^{A}x^{A+B-h+1}\right)\\ &= \frac{1}{1-x}\left(\frac{1-x^{A+1}}{1-x}-x^{B+1}\frac{1-x^{A+1}}{1-x}\right)\\ &= \left(\sum_{i=0}^{A}x^i\right)\left(\sum_{j=0}^{B}x^j\right). \end{align*} The case $B \leq A$ is the same calculation with $A$ and $B$ interchanged. Therefore \begin{align*} \sigma_A\sigma_B = \sum_{h=0}^{\min(A,B)} p^{hK}\sigma_{A+B-2h}, \end{align*} so the constant coefficients also agree. Hence \begin{align*} L_{p^A}L_{p^B} = \sum_{j=0}^{\min(A,B)} p^{jK} L_{p^{A+B-2j}}. \end{align*}[/step]

[guided]The point of the prime-power case for positive coefficient indices is that divisibility by $p$ is measured by one integer: the $p$-adic valuation. Fix a prime $p$, write $K=k-1$, and let $r \in \mathbb{N}$. Write $r=p^R r_0$, where $R \in \mathbb{N}_0$ and $r_0 \in \mathbb{N}$ satisfies $p \nmid r_0$. Applying $L_{p^C}$ to a coefficient indexed by $r$ gives \begin{align*} (L_{p^C}a)_r = \sum_{i=0}^{\min(C,R)} p^{iK} a_{p^{C+R-2i}r_0}. \end{align*} Here $i$ records the exponent in the divisor $p^i$ of both $p^C$ and $r$. Now apply $L_{p^B}$ first and then $L_{p^A}$. The first application changes the exponent of $p$ in the index, so the second divisor parameter is constrained by the new exponent. This gives \begin{align*} (L_{p^A}L_{p^B}a)_r = \sum_{i=0}^{\min(A,R)} \sum_{j=0}^{\min(B,A+R-2i)} p^{(i+j)K} a_{p^{A+B+R-2(i+j)}r_0}. \end{align*} Thus the contribution with total divisor exponent $\ell=i+j$ is counted by admissible pairs $(i,j)$. The right-hand side first chooses an overlap exponent $h$, contributing the scalar $p^{hK}$, and then applies $L_{p^{A+B-2h}}$. Its contribution with total exponent $\ell=h+t$ is counted by admissible pairs $(h,t)$ satisfying \begin{align*} 0 \leq h \leq \min(A,B), \qquad 0 \leq t \leq \min(A+B-2h,R), \qquad h+t=\ell. \end{align*} Equivalently, $h$ ranges through \begin{align*} \max(0,\ell-R) \leq h \leq \min(\ell,A,B,A+B-\ell), \end{align*} whereas the left-hand admissible values of $i$ range through \begin{align*} \max(0,\ell-B) \leq i \leq \min(\ell,A,R,A+R-\ell). \end{align*} The min-max identity proved in the exact argument shows that these two intervals have the same number of integers. Therefore every term \begin{align*} p^{\ell K}a_{p^{A+B+R-2\ell}r_0} \end{align*} appears with the same multiplicity on both sides for every $r \in \mathbb{N}$. The constant coefficient must be handled separately because $0$ has no finite $p$-adic valuation. Define \begin{align*} \sigma_C := \sum_{i=0}^{C} p^{iK} \end{align*} for $C \in \mathbb{N}_0$. Since every divisor of $p^C$ divides $\gcd(p^C,0)$, we have \begin{align*} (L_{p^C}a)_0 = \sigma_Ca_0. \end{align*} Hence the left-hand constant coefficient is $\sigma_A\sigma_Ba_0$. The right-hand constant coefficient is \begin{align*} \left(\sum_{h=0}^{\min(A,B)}p^{hK}\sigma_{A+B-2h}\right)a_0. \end{align*} Let $x:=p^K$. If $A\leq B$, then \begin{align*} \sum_{h=0}^{A}x^h\sum_{t=0}^{A+B-2h}x^t &= \frac{1}{1-x}\left(\sum_{h=0}^{A}x^h-\sum_{h=0}^{A}x^{A+B-h+1}\right)\\ &= \left(\sum_{i=0}^{A}x^i\right)\left(\sum_{j=0}^{B}x^j\right). \end{align*} If $B\leq A$, the same finite geometric-series computation with $A$ and $B$ interchanged gives the same identity. Therefore the constant coefficients agree. Combining the positive-index and constant-index computations proves the prime-power relation.[/guided]

[step:Factor the coefficient operators over distinct primes]Let $M,N \in \mathbb{N}$ with $\gcd(M,N)=1$. For $r \in \mathbb{N}_0$, the coefficient of $L_ML_Na$ at $r$ is \begin{align*} (L_ML_Na)_r = \sum_{e_1 \mid \gcd(M,r)} e_1^{k-1} \sum_{e_2 \mid \gcd(N,Mr/e_1^2)} e_2^{k-1} a_{N(Mr/e_1^2)/e_2^2}. \end{align*} For such $e_1$, every prime divisor of $e_2$ divides $N$ and is therefore coprime to $M$ and $e_1$. Hence \begin{align*} e_2 \mid \frac{Mr}{e_1^2} \quad \Longleftrightarrow \quad e_2 \mid r. \end{align*} Thus the admissible pairs are exactly the pairs \begin{align*} e_1 \mid \gcd(M,r), \qquad e_2 \mid \gcd(N,r). \end{align*} Since $\gcd(M,N)=1$, multiplication gives a bijection from these pairs to divisors $e=e_1e_2$ of $\gcd(MN,r)$. Under this bijection, \begin{align*} e^{k-1}=e_1^{k-1}e_2^{k-1}, \qquad \frac{MNr}{e^2} = \frac{N(Mr/e_1^2)}{e_2^2}. \end{align*} Substituting these identities into the coefficient formula gives \begin{align*} L_ML_N = L_{MN}. \end{align*} The same argument gives $L_NL_M=L_{MN}$, so the coefficient operators attached to coprime indices commute and multiply.[/step]

[guided]When $M$ and $N$ are coprime, their prime factors do not interact, but we still have to check the divisibility condition after the first operator changes the coefficient index. For $r \in \mathbb{N}_0$, applying $L_N$ first and then $L_M$ gives \begin{align*} (L_ML_Na)_r = \sum_{e_1 \mid \gcd(M,r)} e_1^{k-1} \sum_{e_2 \mid \gcd(N,Mr/e_1^2)} e_2^{k-1} a_{N(Mr/e_1^2)/e_2^2}. \end{align*} The outer divisor $e_1$ is the divisor chosen for $L_M$, and the inner divisor $e_2$ is the divisor chosen for $L_N$ after the index has become $Mr/e_1^2$. Because $e_1 \mid M$ and $\gcd(M,N)=1$, every divisor $e_2$ of $N$ is coprime to both $M$ and $e_1$. Therefore the condition \begin{align*} e_2 \mid \frac{Mr}{e_1^2} \end{align*} is equivalent to $e_2 \mid r$: the factors of $M$ and $e_1$ cannot contribute any prime factor of $e_2$. Thus the iterated coefficient formula is indexed exactly by pairs \begin{align*} e_1 \mid \gcd(M,r), \qquad e_2 \mid \gcd(N,r). \end{align*} Multiplication sends such a pair to $e=e_1e_2$, and this is a bijection onto the divisors of $\gcd(MN,r)$ because $M$ and $N$ have disjoint prime divisors. The weight separates as \begin{align*} (e_1e_2)^{k-1}=e_1^{k-1}e_2^{k-1}, \end{align*} and the coefficient index is the same index appearing in $L_{MN}$: \begin{align*} \frac{MNr}{(e_1e_2)^2} = \frac{N(Mr/e_1^2)}{e_2^2}. \end{align*} Therefore each summand in the coefficient formula for $L_ML_N$ corresponds to exactly one summand in the coefficient formula for $L_{MN}$, with the same scalar weight and the same sequence entry. Hence \begin{align*} L_ML_N=L_{MN}. \end{align*}[/guided]

[step:Assemble the formula from prime powers] Write the prime factorizations \begin{align*} m=\prod_p p^{\alpha_p}, \qquad n=\prod_p p^{\beta_p}, \end{align*} where all but finitely many exponents $\alpha_p,\beta_p \in \mathbb{N}_0$ are zero. By the coprime factorisation just proved, \begin{align*} L_mL_n = \prod_p L_{p^{\alpha_p}}L_{p^{\beta_p}}. \end{align*} Using the prime-power relation for each prime $p$ gives \begin{align*} L_mL_n = \prod_p \left( \sum_{\gamma_p=0}^{\min(\alpha_p,\beta_p)} p^{\gamma_p(k-1)} L_{p^{\alpha_p+\beta_p-2\gamma_p}} \right). \end{align*} Expanding the finite product, each choice of exponents $\gamma_p$ determines a divisor \begin{align*} d=\prod_p p^{\gamma_p} \end{align*} of $\gcd(m,n)$, and every divisor of $\gcd(m,n)$ arises uniquely in this way. Moreover, \begin{align*} \prod_p p^{\gamma_p(k-1)}=d^{k-1}, \qquad \prod_p p^{\alpha_p+\beta_p-2\gamma_p} = \frac{mn}{d^2}. \end{align*} Therefore \begin{align*} L_mL_n = \sum_{d \mid \gcd(m,n)} d^{k-1}L_{mn/d^2}. \end{align*} [/step]

[step:Transfer the coefficient identity back to Hecke operators] Let $f \in V$. Applying the coefficient identity to the Fourier coefficient sequence of $f$ gives \begin{align*} T_mT_n f = \sum_{d \mid \gcd(m,n)} d^{k-1}T_{mn/d^2}f. \end{align*} Since this holds for every $f \in V$, the operators are equal in $\operatorname{End}_{\mathbb{C}}(V)$: \begin{align*} T_mT_n = \sum_{d \mid \gcd(m,n)} d^{k-1}T_{mn/d^2}. \end{align*} The argument used only the normalized Hecke action on Fourier expansions, which is valid on both $M_k(\operatorname{SL}_2(\mathbb{Z}))$ and its Hecke-stable subspace $S_k(\operatorname{SL}_2(\mathbb{Z}))$. This proves the theorem. [/step]