How to prove $T_{m}T_{n}=\sum_{d|\gcd(m,n)}d^{k-1}T_{mn/d^2}?$ $T_{n}$ is Hecke operator.

Let $f=\sum_{m=0} ^{\infty} a(m)q^m$ be a Fourier expansion of $f\in M_k(\Gamma(1))$, then we have: \begin{equation}(1) \ \ \ T_n f=\sum_{m=0} ^{\infty} \left(\sum_{d|(n,m)} d^{k-1}a\left(\frac{mn}{d^2}\right)\right)q^m . \end{equation} Let $b(l)$ denote the $l$-th Fourier coefficient of $T_n T_m f$. Then, by (1), \begin{equation}(2) \ \ \ b(l)=\sum_{\substack{{d|n}\\{d|l}}} d^{k-1} \sum_{\substack{{e|m}\\{e|{\frac{nl}{d^2}}}}} e^{k-1} a\left(\frac{mnl}{d^2 e^2}\right). \end{equation} On the other hand, the $l$-th Fourier coefficient $c(l)$ of $\sum_{u|(n,m)} u^{k-1} T_{\frac{mn}{u^2}}f$ is \begin{equation} (3) \ \ \ c(l)=\sum_{\substack{{u|n}\\{u|m}}} u^{k-1} \sum_{\substack{{v|l}\\{v|{\frac{mn}{u^2}}}}} v^{k-1} a\left(\frac{mnl}{u^2 v^2}\right). \end{equation} Consider a change of variable in (2) with $$ (4) \ \ \ u=\frac {e(d,\frac me)}{(e,\frac ld)}, \ \ v=\frac{(e,\frac ld)d}{(d,\frac me)}. $$

This is an invertible mapping between the sets of pairs $(d,e)$ in (2) and $(u,v)$ in (3).

Hence, $b(l)=c(l)$ as desired. Then it completes the proof of $$ T_nT_m=\sum_{u|(n,m)} u^{k-1} T_{\frac{mn}{u^2}}. $$

To see that the change of variable (4) is invertible, we have to show that $(u,v)$ in (4) satisfy the conditions in (3) under the assumptions of $(d,e)$ in (2). This can be shown by $u|e \frac me = m$, $u| \frac nd(d,\frac me) | n$, $v|\frac ld d = l$, $$ u^2 v = \frac{e^2 (d,\frac me)^2}{(e,\frac ld)^2} \frac{(e,\frac ld)d}{(d,\frac me)}=\frac{e^2(d,\frac me) d}{(e,\frac ld)}=ude \Bigg\vert n(d,\frac me) e \Bigg\vert n \frac me e = nm. $$

Also, the same procedure $$ (5) \ \ \ d=\frac{v(u,\frac lv)}{(v,\frac mu)}, \ \ e=\frac{(v,\frac mu) u}{(u,\frac lv)}, $$ takes the pairs $(u,v)$ in (3) to $(d,e)$ in (2).

Moreover, the change of variables (4) and (5) are in fact inverse of each other, since $$ \frac{(u,\frac lv)}{(v,\frac mu)} = \frac{\left( \frac{e(d,\frac me)}{(e,\frac ld)}, \frac{(d,\frac me)l}{(e,\frac ld)d}\right)}{\left(\frac{(e,\frac ld)d}{(d,\frac me)},\frac{(e,\frac ld) m}{(d,\frac me)e}\right)}=\frac{(d,\frac me)}{(e,\frac ld)}. $$

Note that this is essentially interchanging places of $l$ and $m$ in (2). Since places of $n$ and $l$ can be interchanged, the whole problem is equivalent to the result $T_mT_n = T_nT_m$.

The change of variable (4) is not easy to see. However, (4) becomes a lot easier if we impose $(m,n)=1$. So, an alternative method is to derive $T_{mn}=T_mT_n$ when $(m,n)=1$, then prove that $T_{p^s}$ is expressed as a polynomial of $T_p$.