[proofplan]
First we identify the Hecke eigenvalue $\lambda_n$ with the $n$th Fourier coefficient $a_n$ by comparing the coefficient of $q^1$ in the $q$-expansion of $T_n f$. Then we apply the standard Hecke operator multiplication relations to the eigenform $f$. Since the Hecke eigenvalues have already been identified with the Fourier coefficients, the operator identities become the desired multiplicativity and prime-power recurrence for the coefficients.
[/proofplan]
[step:Compare the first Fourier coefficient to identify each Hecke eigenvalue]
For each $n \in \mathbb{N}$, let $\lambda_n \in \mathbb{C}$ denote the scalar satisfying $T_n f = \lambda_n f$. We use the standard $q$-expansion formula for the Hecke operator $T_n$ on $S_k(\mathrm{SL}_2(\mathbb{Z}))$:
\begin{align*}
T_n f
=
\sum_{\ell=1}^{\infty}
\left(
\sum_{d \mid \gcd(n,\ell)} d^{k-1} a_{n\ell/d^2}
\right)q^\ell.
\end{align*}
Taking $\ell=1$, the only divisor of $\gcd(n,1)$ is $d=1$, so the coefficient of $q^1$ in $T_n f$ is $a_n$. On the other hand,
\begin{align*}
\lambda_n f
=
\sum_{\ell=1}^{\infty} \lambda_n a_\ell q^\ell,
\end{align*}
whose coefficient of $q^1$ is $\lambda_n a_1=\lambda_n$, since $f$ is normalized and $a_1=1$. The equality $T_n f=\lambda_n f$ therefore gives $\lambda_n=a_n$. Hence
\begin{align*}
T_n f = a_n f
\end{align*}
for every $n \in \mathbb{N}$.
[guided]
Fix $n \in \mathbb{N}$. Because $f$ is a simultaneous Hecke eigenform, there is a scalar $\lambda_n \in \mathbb{C}$ such that $T_n f=\lambda_n f$. The goal is to prove that this scalar is not a new quantity: it is exactly the Fourier coefficient $a_n$.
We use the standard $q$-expansion formula for the Hecke operator $T_n$:
\begin{align*}
T_n f
=
\sum_{\ell=1}^{\infty}
\left(
\sum_{d \mid \gcd(n,\ell)} d^{k-1} a_{n\ell/d^2}
\right)q^\ell.
\end{align*}
To identify the eigenvalue, we inspect the coefficient of $q^1$. When $\ell=1$, we have $\gcd(n,1)=1$, so the inner divisor sum has only the divisor $d=1$. Therefore the coefficient of $q^1$ in $T_n f$ is
\begin{align*}
1^{k-1}a_{n\cdot 1/1^2}=a_n.
\end{align*}
Now compare this with the eigenvalue equation. Since
\begin{align*}
f=\sum_{\ell=1}^{\infty} a_\ell q^\ell,
\end{align*}
we have
\begin{align*}
\lambda_n f
=
\sum_{\ell=1}^{\infty}\lambda_n a_\ell q^\ell.
\end{align*}
Thus the coefficient of $q^1$ in $\lambda_n f$ is $\lambda_n a_1$. The normalization hypothesis gives $a_1=1$, so this coefficient is $\lambda_n$. Since the two modular forms $T_n f$ and $\lambda_n f$ are equal, their $q$-expansions have equal coefficients. Hence $\lambda_n=a_n$, and therefore $T_n f=a_n f$ for every $n \in \mathbb{N}$.
[/guided]
[/step]
[step:Apply the coprime Hecke product relation to obtain multiplicativity]
Let $m,n \in \mathbb{N}$ satisfy $\gcd(m,n)=1$. We use the standard Hecke operator multiplication formula (citing a result not yet in the wiki: Hecke operator multiplication formula), which in the coprime case gives
\begin{align*}
T_mT_n = T_{mn}.
\end{align*}
Applying both sides to $f$ gives
\begin{align*}
T_m(T_n f)=T_{mn}f.
\end{align*}
Using $T_j f=a_j f$ for $j=m,n,mn$, we obtain
\begin{align*}
a_n T_m f = a_{mn}f,
\end{align*}
and therefore
\begin{align*}
a_n a_m f = a_{mn}f.
\end{align*}
Since $f$ is a nonzero eigenform, equality of scalar multiples of $f$ implies
\begin{align*}
a_{mn}=a_m a_n.
\end{align*}
[guided]
Let $m,n \in \mathbb{N}$ with $\gcd(m,n)=1$. The relevant structural fact about Hecke operators is the standard multiplication formula (citing a result not yet in the wiki: Hecke operator multiplication formula). In the special case where $m$ and $n$ are coprime, the formula reduces to
\begin{align*}
T_mT_n = T_{mn}.
\end{align*}
This identity is an identity of linear operators on $S_k(\mathrm{SL}_2(\mathbb{Z}))$.
Apply the operator identity to the eigenform $f$. We get
\begin{align*}
T_m(T_n f)=T_{mn}f.
\end{align*}
From the previous step, the eigenvalue of $T_j$ acting on $f$ is $a_j$ for every $j \in \mathbb{N}$. Therefore
\begin{align*}
T_n f=a_n f,
\qquad
T_m f=a_m f,
\qquad
T_{mn}f=a_{mn}f.
\end{align*}
Substituting these identities gives
\begin{align*}
T_m(T_n f)
=
T_m(a_n f)
=
a_n T_m f
=
a_n a_m f,
\end{align*}
while the right-hand side is
\begin{align*}
T_{mn}f=a_{mn}f.
\end{align*}
Hence
\begin{align*}
a_n a_m f=a_{mn}f.
\end{align*}
Because $f$ is a nonzero eigenform, the equality of these scalar multiples forces the scalars to be equal:
\begin{align*}
a_{mn}=a_m a_n.
\end{align*}
[/guided]
[/step]
[step:Apply the prime-power Hecke relation to obtain the recursion]
Let $p$ be a prime number and let $r \geq 1$ be an integer. We use the standard prime-power Hecke relation (citing a result not yet in the wiki: Hecke operator multiplication formula):
\begin{align*}
T_pT_{p^r}=T_{p^{r+1}}+p^{k-1}T_{p^{r-1}}.
\end{align*}
Applying both sides to $f$ gives
\begin{align*}
T_p(T_{p^r}f)
=
T_{p^{r+1}}f+p^{k-1}T_{p^{r-1}}f.
\end{align*}
Using $T_j f=a_j f$ for $j=p,p^r,p^{r+1},p^{r-1}$, we obtain
\begin{align*}
a_{p^r}a_p f
=
a_{p^{r+1}}f+p^{k-1}a_{p^{r-1}}f.
\end{align*}
Since $f$ is nonzero, the scalar coefficients agree:
\begin{align*}
a_p a_{p^r}
=
a_{p^{r+1}}+p^{k-1}a_{p^{r-1}}.
\end{align*}
Rearranging gives
\begin{align*}
a_{p^{r+1}}
=
a_p a_{p^r}-p^{k-1}a_{p^{r-1}}.
\end{align*}
This is the stated prime-power recurrence.
[/step]
