[proofplan]
We use the explicit Fourier-coefficient formula for the Hecke operator $T_n$ on weight-$k$ cusp forms. First, comparing the coefficient of $q^1$ in $T_n f=\lambda_n f$ identifies the Hecke eigenvalue $\lambda_n$ with the normalized Fourier coefficient $a_n$. Then comparing the coefficient of $q^m$ in $T_n f=a_n f$ gives the general Hecke coefficient relation. The coprime multiplicativity and the prime-power recursion are the two specializations of that relation when $\gcd(m,n)=1$ and when $(m,n)=(p^r,p)$.
[/proofplan]
[step:Write the Fourier coefficient formula for the Hecke action]
Let
\begin{align*}
g: \mathbb{H} &\to \mathbb{C}
\end{align*}
be a cusp form in $S_k(SL_2(\mathbb{Z}))$ with Fourier expansion
\begin{align*}
g(z)=\sum_{\ell=1}^{\infty} b_\ell q^\ell.
\end{align*}
By the definition of the normalized Hecke operator $T_n$ on weight-$k$ cusp forms for $SL_2(\mathbb{Z})$, the Fourier expansion of $T_n g$ is
\begin{align*}
(T_n g)(z)
=
\sum_{\ell=1}^{\infty}
\left(
\sum_{d \mid \gcd(n,\ell)} d^{k-1} b_{n\ell/d^2}
\right)q^\ell.
\end{align*}
We apply this formula to $g=f$, so the coefficient of $q^\ell$ in $T_n f$ is
\begin{align*}
\sum_{d \mid \gcd(n,\ell)} d^{k-1} a_{n\ell/d^2}.
\end{align*}
[/step]
[step:Compare the first Fourier coefficient to identify the eigenvalues]
Fix $n \in \mathbb{N}$. Since $f$ is a simultaneous Hecke eigenform, there exists $\lambda_n \in \mathbb{C}$ such that
\begin{align*}
T_n f=\lambda_n f.
\end{align*}
Taking the coefficient of $q^1$ on both sides gives
\begin{align*}
\sum_{d \mid \gcd(n,1)} d^{k-1} a_{n/d^2}
=
\lambda_n a_1.
\end{align*}
Because $\gcd(n,1)=1$, the only divisor $d$ is $1$. Since $a_1=1$, this becomes
\begin{align*}
a_n=\lambda_n.
\end{align*}
Thus $T_n f=a_n f$ for every $n \in \mathbb{N}$.
[/step]
[step:Extract the coefficient identity relating $a_m$, $a_n$, and $a_{mn/d^2}$]
Fix $m,n \in \mathbb{N}$. From the identity $T_n f=a_n f$, compare the coefficient of $q^m$. The coefficient of $q^m$ in $T_n f$ is
\begin{align*}
\sum_{d \mid \gcd(m,n)} d^{k-1} a_{mn/d^2},
\end{align*}
while the coefficient of $q^m$ in $a_n f$ is $a_n a_m$. Therefore
\begin{align*}
a_m a_n
=
\sum_{d \mid \gcd(m,n)} d^{k-1} a_{mn/d^2}.
\end{align*}
[guided]
The point of comparing the coefficient of $q^m$ is that the Hecke operator $T_n$ mixes the Fourier coefficients of $f$ in a controlled divisor sum. Since we already proved that the eigenvalue of $T_n$ is exactly $a_n$, the equality $T_n f=a_n f$ can be read coefficient by coefficient.
Fix $m,n \in \mathbb{N}$. The coefficient formula for $T_n f$ gives the coefficient of $q^m$ as
\begin{align*}
\sum_{d \mid \gcd(m,n)} d^{k-1} a_{mn/d^2}.
\end{align*}
On the other hand, multiplying the [Fourier series](/page/Fourier%20Series) of $f$ by the scalar $a_n$ gives
\begin{align*}
a_n f(z)
=
\sum_{\ell=1}^{\infty} a_n a_\ell q^\ell,
\end{align*}
so the coefficient of $q^m$ in $a_n f$ is $a_n a_m$. Equality of Fourier expansions therefore gives
\begin{align*}
a_m a_n
=
\sum_{d \mid \gcd(m,n)} d^{k-1} a_{mn/d^2}.
\end{align*}
This identity is the common source of both the coprime multiplicativity and the prime-power recursion.
[/guided]
[/step]
[step:Specialize the coefficient identity to coprime indices]
Assume that $m,n \in \mathbb{N}$ satisfy $\gcd(m,n)=1$. Then the only positive divisor $d$ of $\gcd(m,n)$ is $d=1$. The coefficient identity becomes
\begin{align*}
a_m a_n
=
1^{k-1}a_{mn}
=
a_{mn}.
\end{align*}
Hence
\begin{align*}
a_{mn}=a_m a_n.
\end{align*}
[/step]
[step:Specialize the coefficient identity to powers of a prime]
Let $p$ be a prime number, and let $r \ge 1$ be an integer. Apply the coefficient identity with $m=p^r$ and $n=p$. Since
\begin{align*}
\gcd(p^r,p)=p,
\end{align*}
the positive divisors of $\gcd(p^r,p)$ are $1$ and $p$. Thus
\begin{align*}
a_{p^r}a_p
&=
1^{k-1}a_{p^{r+1}}
+
p^{k-1}a_{p^{r+1}/p^2} \\
&=
a_{p^{r+1}}+p^{k-1}a_{p^{r-1}}.
\end{align*}
Rearranging this equality gives
\begin{align*}
a_{p^{r+1}}=a_p a_{p^r}-p^{k-1}a_{p^{r-1}}.
\end{align*}
This is the asserted prime-power recursion.
[/step]
