[proofplan]
We expand each summand in the triangular formula for $T_n f$ using the Fourier expansion of $f$. The finite sum over $b$ is a roots-of-unity filter: it keeps exactly those Fourier terms whose index is divisible by $d$. After writing the surviving index as $r=dt$, the contribution from the pair $(a,d)$ has only powers $q^{at}$. Extracting the coefficient of $q^m$ forces $a\mid m$, and reindexing the pairs $ad=n$ gives the stated divisor sum.
[/proofplan]
[step:Insert the Fourier expansion into the triangular formula]
For each factorization $ad=n$ with $a,d\in\mathbb{N}$ and each $b\in\{0,\dots,d-1\}$, the Fourier expansion of $f$ at the cusp gives
\begin{align*}
f\left(\frac{az+b}{d}\right)
&=
\sum_{r=0}^{\infty}
a_r(f)
\exp\left(2\pi i r\frac{az+b}{d}\right) \\
&=
\sum_{r=0}^{\infty}
a_r(f)
\exp\left(\frac{2\pi i rb}{d}\right)
q^{ar/d}.
\end{align*}
Substituting this into the definition of $T_n f$ gives
\begin{align*}
(T_n f)(z)
&=
\sum_{\substack{a,d\in \mathbb{N}\\ ad=n}}
\frac{a^{k-1}}{d}
\sum_{b=0}^{d-1}
\sum_{r=0}^{\infty}
a_r(f)
\exp\left(\frac{2\pi i rb}{d}\right)
q^{ar/d}.
\end{align*}
Since the outer sums over $a,d,b$ are finite and the [Fourier series](/page/Fourier%20Series) is convergent for $\operatorname{Im}(z)>0$, this substitution is legitimate term by term.
[guided]
Fix a factorization $ad=n$ with $a,d\in\mathbb{N}$ and fix $b\in\{0,\dots,d-1\}$. The variable in the Fourier expansion of $f$ is
\begin{align*}
q\left(\frac{az+b}{d}\right)
&=
\exp\left(2\pi i\frac{az+b}{d}\right) \\
&=
\exp\left(\frac{2\pi i b}{d}\right)q^{a/d}.
\end{align*}
Therefore the $r$-th Fourier term of $f$ becomes
\begin{align*}
a_r(f)
\left(q\left(\frac{az+b}{d}\right)\right)^r
=
a_r(f)
\exp\left(\frac{2\pi i rb}{d}\right)
q^{ar/d}.
\end{align*}
Substituting this expression into the triangular formula for $T_n$ yields
\begin{align*}
(T_n f)(z)
&=
\sum_{\substack{a,d\in \mathbb{N}\\ ad=n}}
\frac{a^{k-1}}{d}
\sum_{b=0}^{d-1}
\sum_{r=0}^{\infty}
a_r(f)
\exp\left(\frac{2\pi i rb}{d}\right)
q^{ar/d}.
\end{align*}
The only infinite sum here is the Fourier expansion of $f$, which converges for $\operatorname{Im}(z)>0$. The remaining sums are finite, so no limiting interchange beyond substituting the convergent Fourier series into finitely many summands is involved.
[/guided]
[/step]
[step:Use the roots-of-unity sum to keep only indices divisible by $d$]
For $d\in\mathbb{N}$ and $r\in\mathbb{N}\cup\{0\}$, define
\begin{align*}
S_d(r):=\sum_{b=0}^{d-1}\exp\left(\frac{2\pi i rb}{d}\right).
\end{align*}
If $d\mid r$, then every summand equals $1$, so $S_d(r)=d$. If $d\nmid r$, set $\omega:=\exp(2\pi i r/d)$. Then $\omega\neq 1$ and $\omega^d=1$, hence the finite geometric sum gives
\begin{align*}
S_d(r)
=
\sum_{b=0}^{d-1}\omega^b
=
\frac{1-\omega^d}{1-\omega}
=
0.
\end{align*}
Thus
\begin{align*}
S_d(r)
=
\begin{cases}
d, & d\mid r,\\
0, & d\nmid r.
\end{cases}
\end{align*}
It follows that, for each fixed pair $(a,d)$ with $ad=n$, only indices $r=dt$ with $t\in\mathbb{N}\cup\{0\}$ survive. Hence the contribution of this pair to $T_n f$ is
\begin{align*}
\frac{a^{k-1}}{d}
\sum_{t=0}^{\infty}
a_{dt}(f)dq^{at}
=
a^{k-1}
\sum_{t=0}^{\infty}
a_{dt}(f)q^{at}.
\end{align*}
[guided]
The sum over $b$ is the mechanism that removes all fractional powers of $q$. Define
\begin{align*}
S_d(r):=\sum_{b=0}^{d-1}\exp\left(\frac{2\pi i rb}{d}\right).
\end{align*}
If $d\mid r$, then $r/d$ is an integer, so
\begin{align*}
\exp\left(\frac{2\pi i rb}{d}\right)=1
\end{align*}
for every $b\in\{0,\dots,d-1\}$. Therefore $S_d(r)=d$.
If $d\nmid r$, set $\omega:=\exp(2\pi i r/d)$. Then $\omega^d=\exp(2\pi i r)=1$, while $\omega\neq 1$ because $r/d$ is not an integer. The finite geometric-series identity gives
\begin{align*}
S_d(r)
=
\sum_{b=0}^{d-1}\omega^b
=
\frac{1-\omega^d}{1-\omega}
=
0.
\end{align*}
Thus the $b$-sum kills precisely the terms with $d\nmid r$ and multiplies the terms with $d\mid r$ by $d$.
Now write each surviving index as $r=dt$, where $t\in\mathbb{N}\cup\{0\}$. The corresponding power of $q$ is
\begin{align*}
q^{ar/d}=q^{a(dt)/d}=q^{at},
\end{align*}
and the prefactor $\frac{a^{k-1}}{d}$ cancels the factor $d$ from the roots-of-unity sum. Therefore the total contribution from the fixed factorization $ad=n$ is
\begin{align*}
a^{k-1}
\sum_{t=0}^{\infty}
a_{dt}(f)q^{at}.
\end{align*}
[/guided]
[/step]
[step:Extract the coefficient of $q^m$]
Combining the previous step over all factorizations $ad=n$ gives
\begin{align*}
(T_n f)(z)
=
\sum_{\substack{a,d\in \mathbb{N}\\ ad=n}}
a^{k-1}
\sum_{t=0}^{\infty}
a_{dt}(f)q^{at}.
\end{align*}
Fix $m\in\mathbb{N}\cup\{0\}$. The term indexed by $(a,d,t)$ contributes to the coefficient of $q^m$ exactly when $at=m$. This is equivalent to $a\mid m$ and $t=m/a$, with the convention that every $a\in\mathbb{N}$ divides $0$. Therefore
\begin{align*}
a_m(T_n f)
&=
\sum_{\substack{a,d\in\mathbb{N}\\ ad=n\\ a\mid m}}
a^{k-1}a_{d(m/a)}(f).
\end{align*}
Since $ad=n$, the index of the Fourier coefficient satisfies
\begin{align*}
d\frac{m}{a}
=
\frac{n}{a}\frac{m}{a}
=
\frac{mn}{a^2}.
\end{align*}
Thus
\begin{align*}
a_m(T_n f)
&=
\sum_{\substack{a\mid n\\ a\mid m}}
a^{k-1}a_{mn/a^2}(f) \\
&=
\sum_{a\mid \gcd(m,n)}
a^{k-1}a_{mn/a^2}(f).
\end{align*}
Renaming the divisor $a$ as $c$ gives
\begin{align*}
a_m(T_n f)
=
\sum_{c\mid \gcd(m,n)} c^{k-1}a_{mn/c^2}(f),
\end{align*}
which is the desired formula.
[guided]
After the roots-of-unity filter, the Hecke transform has the expansion
\begin{align*}
(T_n f)(z)
=
\sum_{\substack{a,d\in \mathbb{N}\\ ad=n}}
a^{k-1}
\sum_{t=0}^{\infty}
a_{dt}(f)q^{at}.
\end{align*}
We now identify which terms contribute to a fixed coefficient $a_m(T_n f)$. A term in the inner sum has exponent $at$, so it contributes to the coefficient of $q^m$ exactly when
\begin{align*}
at=m.
\end{align*}
For fixed $a$, this equation has a solution $t\in\mathbb{N}\cup\{0\}$ exactly when $a\mid m$, in which case $t=m/a$. When $m=0$, this condition is satisfied for every positive integer $a$, and $t=0$.
Therefore the coefficient of $q^m$ is
\begin{align*}
a_m(T_n f)
&=
\sum_{\substack{a,d\in\mathbb{N}\\ ad=n\\ a\mid m}}
a^{k-1}a_{d(m/a)}(f).
\end{align*}
The condition $ad=n$ means $d=n/a$, so the Fourier index becomes
\begin{align*}
d\frac{m}{a}
=
\frac{n}{a}\frac{m}{a}
=
\frac{mn}{a^2}.
\end{align*}
The simultaneous conditions $a\mid n$ and $a\mid m$ are exactly the condition $a\mid \gcd(m,n)$, with $\gcd(0,n)=n$. Hence
\begin{align*}
a_m(T_n f)
&=
\sum_{a\mid \gcd(m,n)}
a^{k-1}a_{mn/a^2}(f).
\end{align*}
Finally, the letter used for the divisor is immaterial; writing $c$ instead of $a$ gives
\begin{align*}
a_m(T_n f)
=
\sum_{c\mid \gcd(m,n)} c^{k-1}a_{mn/c^2}(f).
\end{align*}
This is precisely the stated coefficient formula.
[/guided]
[/step]
