How to prove $T_{m}T_{n}=\sum_{d|\gcd(m,n)}d^{k-1}T_{mn/d^2}?$ $T_{n}$ is Hecke operator. 
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$$
 A: 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$. 
A: 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.

First part (Theorem 6.12): If $m,n$ are positive integers with $(m,n)=1$ the following is valid
  \begin{align*}
T_nT_m=T_{mn}\tag{1}
\end{align*}

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).

Second part (Theorem 6.13): Any two Hecke operators $T(n)$ an $T(m)$ defined on $M_k$ commute with each other. Moreover, we have the composition formula
  \begin{align*}
T(m)T(n)=\sum_{d|(m,n)}d^{k-1}T(mn/d^2)\tag{3}
\end{align*}

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.

Note: It might be interesting to note that multiplicative  arithmetic functions $f$ which share the property corresponding to (3) are called specially multiplicative functions.

According to Introduction to Arithmetical Functions by P. J. McCarthy we have:

If $f$ is a multiplicative function then the following statements are equivalent:
  
  
*
  
*$f$ is a convolution of two completely multiplicative functions.
  
*There is a multiplicative function $F$ such that for all $m$ and $n$,
  \begin{align*}
f(mn)=\sum_{d|(m,n)}f\left(\frac{m}{d}\right)f\left(\frac{n}{d}\right)F(d)
\end{align*}
  
*There is a completely multiplicative function $B$ such that for all $m$ and $n$
  \begin{align*}
f(m)f(n)=\sum_{d|(m,n)}f\left(\frac{mn}{d^2}\right)B(d)
\end{align*}
  
*For all primes $p$ and all $r \geq 1$,
  \begin{align*}
f(p^{r+1})=f(p)f(p^r)+f(p^{r-1})\left(f(p^2)-f(p)^2\right)
\end{align*}

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*}
