Define Hecke operator is $$T_{n}: M_{k}(\Gamma(1))\rightarrow M_{k}(\Gamma(1))$$ $$ T_{n}: f \longmapsto n^{k-1}\sum_{a,d>0, ad=n}\sum_{b=0}^{d-1} d^{-k}f(\frac{a\tau+b}{d}).$$ How to prove $T_{m}T_{n}=\sum_{d|\gcd(m,n)}d^{k-1}T_{mn/d^2}?$
I have proved that $$(T_{n}f)(\tau)=\sum_{a,d>0, ad=n}\Big(\sum_{d|\gcd(m,n)}d^{k-1} a\big(\frac{mn}{d^2}\big)\Big)q^m$$
The proof is rather technical. Here we show the multiplicative property of the Hecke operators by following mostly verbatim Theorem 6.12 and 6.13 of T. M. Apostol's Modular Functions and Dirichlet Series in Number Theory.
The proof is organised in two parts. First it is shown for positive integers $m,n$ relatively prime and then for prime powers.
If $f\in M_k$ we have \begin{align*} (T_nf)(\tau)=\frac{1}{n}\sum_{{a\ge 1,ad=n}\atop{0\leq b<d}}a^kf(A\tau) \end{align*} where $A=\begin{pmatrix}a&b\\0&d\end{pmatrix}$. Applying $T_m$ to each member we have \begin{align*} \{T_m(T_n(f))\}(\tau)=\frac{1}{m}\sum_{{\alpha\geq 1,\alpha\delta=m}\atop{0\leq \beta<\delta}} \alpha^k\frac{1}{n}\sum_{{a\geq 1,ad=m}\atop{0\leq b<d}}a^kf(BA\tau) \end{align*} where $B=\begin{pmatrix}\alpha&\beta\\0&\delta\end{pmatrix}$. This can be written as \begin{align*} \{T_m(T_n(f))\}(\tau)=\frac{1}{mn}\sum_{{\alpha\geq 1,\alpha\delta=m}\atop{0\leq \beta<\delta}} \sum_{{a\geq 1,ad=m}\atop{0\leq b<d}}(\alpha a)^kf(C\tau)\tag{2} \end{align*} where \begin{align*} C=BA=\begin{pmatrix}\alpha&\beta\\0&\delta\end{pmatrix}\begin{pmatrix}a&b\\0&d\end{pmatrix} =\begin{pmatrix}\alpha a&\alpha b+\beta d\\0&d\delta\end{pmatrix} \end{align*} As $d$ and $\delta$ run through the positive divisors of $n$ and $m$, respectively, the product $d\delta$ runs through the positive divisors of $mn$ since $(m,n)=1$. The linear combination $\alpha b+\beta d$ runs through a complete residue system mod $d\delta$ as $b$ and $\beta$ run through complete residue systems mod $d$ and $\delta$, respectively. Therefore the matrix $C$ runs through a complete set of non-equivalent elements of $\Gamma(mn)$ and we see that (2) implies (1).
If $(m,n)=1$ the formula (3) reduces to (1). Therefore it suffices to treat the case when $m$ and $n$ are powers of the same prime $p$. First we consider the case $m=p$ and $n=p^r$, where $r \geq 1$. In this case we are to prove that \begin{align*} T(p)T(p^r)=T(p^{r+1})+p^{k-1}T(p^{r-1}) \end{align*} Using the representation \begin{align*} (T(n)f)(\tau)=n^{k-1}\sum_{{a\geq 1,ad=n}\atop{0\leq b<d}}d^{-k}f\left(\frac{a\tau+b}{d}\right) =\frac{1}{n}\sum_{{a\geq 1,ad=n}\atop{0\leq b<d}}a^k f\left(\frac{a\tau+b}{d}\right)\tag{4} \end{align*} we have \begin{align*} \{T(p^r)f\}(\tau)=p^{-r}\sum_{{0\leq t\leq r}\atop{0\leq b_t<p^t}}p^{(r-t)k}f\left(\frac{p^{r-t}\tau+b_t}{p^t}\right)\tag{5} \end{align*} From the representation \begin{align*} (T(n)f)(\tau)=n^{k-1}\sum_{d|n}d^{-k}\sum_{b=0}^{d-1}f\left(\frac{n\tau+bd}{d^2}\right) \end{align*} we have in the special case when $n$ is prime, say $n=p$ \begin{align*} (T(p)f)(\tau)=p^{k-1}f(p\tau)+\frac{1}{p}\sum_{b=0}^{p-1}f\left(\frac{\tau+b}{p}\right) \end{align*} and so when we apply $T(p)$ to each member of (5) we find \begin{align*} \{T(p)T(p^r)f\}(\tau)&=p^{k-1-r}\sum_{{0\leq t\leq r}\atop{0\leq b_t<p^t}}p^{(r-t)k}f\left(\frac{p^{r+1-t}\tau+pb_t}{p^t}\right)\\ &\qquad +p^{-1-r}\sum_{{0\leq t\leq r}\atop{0\leq b_t<p^t}}p^{(r-t)k}\sum_{b=0}^{p-1}f\left(\frac{p^{r-t}\tau+b_t+bp^t}{p^{t+1}}\right) \end{align*} In the second sum the linear combination $b_t+bp^t$ runs through a complete residue system mod $p^{t+1}$. Since $r-t=(r+1)-(t+1)$ the second sum, together with the term $t=0$ from the first sum, is equal to $\{T(p^{r+1}f\}(\tau)$. In the remaining terms we cancel a factor $p$ in the argument of $f$, then transfer the facctor $p^k$ to each summand to obtain \begin{align*} \{T(p)T(p^r)f\}(\tau)=\{T(p^{r+1})f\}(\tau)+p^{-1-r}\sum_{{1\leq t\leq r}\atop{0\leq b_t<p^t}}p^{(r+1-t)k}f\left(\frac{p^{r-t}\tau+b_t}{p^{t-1}}\right). \end{align*} Dividing each $b_t$ b $p^{t-1}$ we can write \begin{align*} b_t=q_tp^{t-1}+r_t \end{align*} where $0\leq r_t<p^{t-1}$ and $q_t$ runs through a complete residue system mod $p$. Since $f$ is periodic with period $1$ we have \begin{align*} f\left(\frac{p^{r-t}\tau+b_t}{p^{t-1}}\right)=f\left(\frac{p^{r-t}\tau+r_t}{p^{t-1}}\right), \end{align*} so as $q_t$ runs through a complete residue system mod $p$ each term is repeated $p$ times. Replacing the index $t$ by $t-1$ we see the last sum is $p^{k-1}$ times the sum defining $\{T(p^{r-1})f\}(\tau)$. This proves (4).
The last step is to consider general powers of the same prime, say $m=p^s$ and $n=p^r$. Assuming without loss of generality that $r\leq s$ it can be shown by induction on $r$ that \begin{align*} T(p^r)T(p^s)=\sum_{t=0}^rp^{t(k-1)}T(p^{r+s-2t})=\sum_{d|(p^r,p^s)}d^{k-1}T\left(\frac{p^{r+s}}{d^2}\right) \end{align*} for all $r$ and all $s\geq r$ and the claim (3) follows.
According to Introduction to Arithmetical Functions by P. J. McCarthy we have:
Prominent members of this function class are the divisor function $\sigma_k$ fulfilling the so-called Busche-Ramanujan identity \begin{align*} \sigma_k(m)\sigma_k(n)=\sum_{d|(m,n)}d^k\sigma_k\left(\frac{mn}{d^2}\right) \end{align*} and the Ramanujan Tau function with \begin{align*} \tau(m)\tau(n)=\sum_{d|(m,n)}d^{11}\tau\left(\frac{mn}{d^2}\right) \end{align*}