[proofplan]
We equip $S_k(SL_2(\mathbb{Z}))$ with the Petersson inner product and use the standard Hecke-theoretic facts that the [Hecke operators preserve cusp forms](/theorems/4245), are self-adjoint for this inner product, and satisfy the usual multiplication relations. Since the space is finite-dimensional, the commuting normal operators generated by the prime-index Hecke operators admit a common orthogonal eigenbasis by the finite-dimensional simultaneous spectral theorem. The Hecke relations then imply that each vector in this basis is an eigenvector for every $T_n$. Finally, the Fourier-coefficient formula for $T_n$ shows that a nonzero level-one Hecke eigenform has nonzero first Fourier coefficient, so it can be normalized by scaling.
[/proofplan]
[step:Put the cusp form space in finite-dimensional Hermitian form]
Let
\begin{align*}
V := S_k(SL_2(\mathbb{Z}))
\end{align*}
be the complex [vector space](/page/Vector%20Space) of cusp forms of weight $k$ and level one. This space is finite-dimensional by the standard finite-dimensionality theorem for spaces of modular forms (citing a result not yet in the wiki: finite-dimensionality of spaces of modular forms).
Let $\mathcal{F} \subset \mathbb{H}$ be a measurable fundamental domain for the action of $SL_2(\mathbb{Z})$ on the upper half-plane $\mathbb{H}$. Define the hyperbolic measure $\mu_{\mathrm{hyp}}$ on $\mathbb{H}$ by
\begin{align*}
d\mu_{\mathrm{hyp}}(z) := y^{-2}\,d\mathcal{L}^2(x,y),
\end{align*}
where $z = x + iy$ and $\mathcal{L}^2$ denotes two-dimensional Lebesgue measure on $\mathbb{R}^2$. The Petersson inner product on $V$ is the Hermitian form
\begin{align*}
(\cdot,\cdot)_{\mathrm{Pet}}: V \times V &\to \mathbb{C}, \\
(f,g) &\mapsto \int_{\mathcal{F}} f(z)\overline{g(z)}\,y^k\,d\mu_{\mathrm{hyp}}(z).
\end{align*}
For cusp forms this integral converges, and $(\cdot,\cdot)_{\mathrm{Pet}}$ is positive definite. Thus $V$ is a finite-dimensional Hermitian [inner product space](/page/Inner%20Product%20Space).
[guided]
The first point is to put the problem into finite-dimensional linear algebra. Define
\begin{align*}
V := S_k(SL_2(\mathbb{Z})).
\end{align*}
The standard finite-dimensionality theorem for spaces of modular forms says that $V$ is finite-dimensional over $\mathbb{C}$ (citing a result not yet in the wiki: finite-dimensionality of spaces of modular forms).
We now choose the inner product with respect to which Hecke operators have good adjointness properties. Let $\mathcal{F} \subset \mathbb{H}$ be a measurable fundamental domain for $SL_2(\mathbb{Z}) \curvearrowright \mathbb{H}$. Write $z = x + iy$, and let $\mathcal{L}^2$ denote two-dimensional Lebesgue measure on $\mathbb{R}^2$. Define the hyperbolic measure $\mu_{\mathrm{hyp}}$ by
\begin{align*}
d\mu_{\mathrm{hyp}}(z) := y^{-2}\,d\mathcal{L}^2(x,y).
\end{align*}
The Petersson inner product is the map
\begin{align*}
(\cdot,\cdot)_{\mathrm{Pet}}: V \times V &\to \mathbb{C}, \\
(f,g) &\mapsto \int_{\mathcal{F}} f(z)\overline{g(z)}\,y^k\,d\mu_{\mathrm{hyp}}(z).
\end{align*}
Because $f$ and $g$ are cusp forms, the exponential decay at the cusp gives convergence of the integral. The standard positivity argument for the Petersson product gives $(f,f)_{\mathrm{Pet}} > 0$ whenever $f \neq 0$. Hence $V$ is a finite-dimensional Hermitian inner product space, which is exactly the setting needed for spectral theory.
[/guided]
[/step]
[step:Diagonalize the commuting prime-index Hecke operators]
For each $n \in \mathbb{N}$, let
\begin{align*}
T_n: V \to V
\end{align*}
denote the $n$th Hecke operator. The standard adjointness theorem for level-one Hecke operators gives
\begin{align*}
(T_n f,g)_{\mathrm{Pet}} = (f,T_n g)_{\mathrm{Pet}}
\end{align*}
for all $f,g \in V$ and all $n \in \mathbb{N}$; hence every $T_n$ is self-adjoint and therefore normal.
For primes $p$ and $q$, the Hecke multiplication relation gives
\begin{align*}
T_pT_q = T_qT_p.
\end{align*}
Thus the family $\{T_p : p \text{ prime}\}$ is a commuting family of normal operators on the finite-dimensional Hermitian space $V$. By the finite-dimensional simultaneous spectral theorem for commuting normal operators (citing a result not yet in the wiki: simultaneous diagonalization of commuting normal operators), there exists an orthogonal basis
\begin{align*}
f_1,\dots,f_d \in V
\end{align*}
such that, for every prime $p$ and every $j \in \{1,\dots,d\}$, there is a scalar $\lambda_j(p) \in \mathbb{C}$ satisfying
\begin{align*}
T_p f_j = \lambda_j(p) f_j.
\end{align*}
If $d = \dim_{\mathbb{C}} V = 0$, this statement means that the empty family is a basis, and the theorem is vacuous.
[guided]
We now use the compatibility between the Petersson product and Hecke operators. For each $n \in \mathbb{N}$, the Hecke operator is a complex-[linear map](/page/Linear%20Map)
\begin{align*}
T_n: V \to V.
\end{align*}
The standard adjointness theorem for level-one Hecke operators states that
\begin{align*}
(T_n f,g)_{\mathrm{Pet}} = (f,T_n g)_{\mathrm{Pet}}
\end{align*}
for all $f,g \in V$ (citing a result not yet in the wiki: Petersson self-adjointness of level-one Hecke operators). Therefore each $T_n$ is self-adjoint. A self-adjoint operator on a Hermitian inner product space is normal, because $T_nT_n^* = T_n^2 = T_n^*T_n$.
The next input is commutativity. The Hecke multiplication relation says that for positive integers $m,n$,
\begin{align*}
T_mT_n = \sum_{r \mid \gcd(m,n)} r^{k-1} T_{mn/r^2}.
\end{align*}
In particular, if $p$ and $q$ are primes, then this formula is symmetric in $p$ and $q$, so
\begin{align*}
T_pT_q = T_qT_p.
\end{align*}
Thus the prime-index Hecke operators form a commuting family of normal operators on the finite-dimensional Hermitian space $V$.
We may now apply the finite-dimensional simultaneous spectral theorem for commuting normal operators (citing a result not yet in the wiki: simultaneous diagonalization of commuting normal operators). Its hypotheses are exactly satisfied: the space is finite-dimensional and Hermitian, each $T_p$ is normal, and the family $\{T_p : p \text{ prime}\}$ commutes pairwise. Hence there is an orthogonal basis
\begin{align*}
f_1,\dots,f_d \in V
\end{align*}
such that every $f_j$ is an eigenvector for every prime-index operator. That is, for each prime $p$ and each $j \in \{1,\dots,d\}$, there exists $\lambda_j(p) \in \mathbb{C}$ with
\begin{align*}
T_p f_j = \lambda_j(p) f_j.
\end{align*}
If $V = \{0\}$, then $d = 0$ and the empty basis already proves the theorem.
[/guided]
[/step]
[step:Use the Hecke relations to get eigenvectors for every $T_n$]
Fix $j \in \{1,\dots,d\}$. We prove that $f_j$ is an eigenvector for $T_n$ for every $n \in \mathbb{N}$.
For a prime $p$ and an integer $r \geq 1$, the Hecke recursion gives
\begin{align*}
T_{p^{r+1}} = T_pT_{p^r} - p^{k-1}T_{p^{r-1}},
\end{align*}
with $T_{p^0} = T_1 = \operatorname{id}_V$. By induction on $r$, each $T_{p^r}f_j$ is a scalar multiple of $f_j$.
Now let
\begin{align*}
n = \prod_{\ell=1}^{s} p_\ell^{r_\ell}
\end{align*}
be the prime factorization of $n$, where $p_1,\dots,p_s$ are distinct primes and $r_1,\dots,r_s \in \mathbb{N}$. The multiplicativity relation for coprime indices gives
\begin{align*}
T_n = T_{p_1^{r_1}}\cdots T_{p_s^{r_s}}.
\end{align*}
Since each factor sends $f_j$ to a scalar multiple of $f_j$, their product does too. Therefore there exists $\lambda_j(n) \in \mathbb{C}$ such that
\begin{align*}
T_n f_j = \lambda_j(n) f_j.
\end{align*}
Thus the basis $f_1,\dots,f_d$ consists of simultaneous eigenforms for all Hecke operators $T_n$.
[guided]
The spectral theorem gave eigenvectors for the operators $T_p$ with $p$ prime. We must now explain why this implies the same statement for every $T_n$.
Fix one basis vector $f_j$. First consider prime powers. For a prime $p$ and an integer $r \geq 1$, the Hecke recursion is
\begin{align*}
T_{p^{r+1}} = T_pT_{p^r} - p^{k-1}T_{p^{r-1}},
\end{align*}
with $T_{p^0} = T_1 = \operatorname{id}_V$. We already know that $T_p f_j = \lambda_j(p)f_j$. The case $r=0$ is also an eigenvalue statement, since
\begin{align*}
T_1f_j = f_j.
\end{align*}
Assume inductively that $T_{p^r}f_j = \lambda_j(p^r)f_j$ and $T_{p^{r-1}}f_j = \lambda_j(p^{r-1})f_j$ for some scalars $\lambda_j(p^r),\lambda_j(p^{r-1}) \in \mathbb{C}$. Applying the recursion to $f_j$ gives
\begin{align*}
T_{p^{r+1}}f_j
&= T_pT_{p^r}f_j - p^{k-1}T_{p^{r-1}}f_j \\
&= T_p\bigl(\lambda_j(p^r)f_j\bigr) - p^{k-1}\lambda_j(p^{r-1})f_j \\
&= \lambda_j(p^r)\lambda_j(p)f_j - p^{k-1}\lambda_j(p^{r-1})f_j \\
&= \bigl(\lambda_j(p^r)\lambda_j(p)-p^{k-1}\lambda_j(p^{r-1})\bigr)f_j.
\end{align*}
So $f_j$ is an eigenvector for $T_{p^{r+1}}$.
Now let $n \in \mathbb{N}$ have prime factorization
\begin{align*}
n = \prod_{\ell=1}^{s} p_\ell^{r_\ell},
\end{align*}
where the primes $p_1,\dots,p_s$ are distinct and $r_1,\dots,r_s \in \mathbb{N}$. The Hecke multiplicativity relation for coprime indices gives
\begin{align*}
T_n = T_{p_1^{r_1}}\cdots T_{p_s^{r_s}}.
\end{align*}
Each factor sends $f_j$ to a scalar multiple of $f_j$, so the product also sends $f_j$ to a scalar multiple of $f_j$. Hence there is a scalar $\lambda_j(n) \in \mathbb{C}$ such that
\begin{align*}
T_n f_j = \lambda_j(n) f_j.
\end{align*}
This proves that every vector in the basis is a simultaneous eigenform for all Hecke operators, not only the prime-index ones.
[/guided]
[/step]
[step:Show every nonzero simultaneous eigenform has nonzero first Fourier coefficient]
Let $f \in V$ be a nonzero simultaneous Hecke eigenform. Write its Fourier expansion at infinity as
\begin{align*}
f(z) = \sum_{m=1}^{\infty} a_m q^m, \qquad q := e^{2\pi i z},
\end{align*}
where $a_m \in \mathbb{C}$. Since $f$ is nonzero, choose the least integer $m_0 \in \mathbb{N}$ such that $a_{m_0} \neq 0$.
For the normalized level-one Hecke operator $T_{m_0}$, the coefficient of $q^1$ in $T_{m_0}f$ is $a_{m_0}$. Since $f$ is a simultaneous eigenform, there exists $\lambda(m_0) \in \mathbb{C}$ such that
\begin{align*}
T_{m_0}f = \lambda(m_0)f.
\end{align*}
Comparing the coefficient of $q^1$ gives
\begin{align*}
a_{m_0} = \lambda(m_0)a_1.
\end{align*}
Because $a_{m_0} \neq 0$, it follows that $a_1 \neq 0$.
[guided]
Let $f \in V$ be a nonzero simultaneous Hecke eigenform. Since $f$ is a cusp form, its Fourier expansion at the cusp infinity has no constant term, so we may write
\begin{align*}
f(z) = \sum_{m=1}^{\infty} a_m q^m, \qquad q := e^{2\pi i z},
\end{align*}
with coefficients $a_m \in \mathbb{C}$. Because $f \neq 0$, not all Fourier coefficients vanish. Let $m_0 \in \mathbb{N}$ be the least integer such that
\begin{align*}
a_{m_0} \neq 0.
\end{align*}
We use the Fourier-coefficient formula for the normalized level-one Hecke operator. If
\begin{align*}
T_n f(z) = \sum_{r=1}^{\infty} b_r q^r,
\end{align*}
then
\begin{align*}
b_r = \sum_{d \mid \gcd(n,r)} d^{k-1}a_{nr/d^2}.
\end{align*}
Apply this with $n=m_0$ and $r=1$. Since $\gcd(m_0,1)=1$, the only divisor is $d=1$, so the coefficient of $q^1$ in $T_{m_0}f$ is
\begin{align*}
b_1 = a_{m_0}.
\end{align*}
Because $f$ is a simultaneous eigenform, it is in particular an eigenform for $T_{m_0}$. Hence there exists $\lambda(m_0) \in \mathbb{C}$ such that
\begin{align*}
T_{m_0}f = \lambda(m_0)f.
\end{align*}
The coefficient of $q^1$ on the right-hand side is $\lambda(m_0)a_1$. Comparing the two $q^1$ coefficients gives
\begin{align*}
a_{m_0} = \lambda(m_0)a_1.
\end{align*}
The left-hand side is nonzero by the choice of $m_0$, so $a_1 \neq 0$.
[/guided]
[/step]
[step:Normalize the simultaneous eigenbasis]
For each basis element $f_j$, write
\begin{align*}
f_j(z) = \sum_{m=1}^{\infty} a_j(m)q^m.
\end{align*}
By the previous step, $a_j(1) \neq 0$. Define
\begin{align*}
g_j := a_j(1)^{-1}f_j.
\end{align*}
Then $g_j$ is still a simultaneous eigenform for every $T_n$, because scalar multiples of eigenvectors have the same eigenvalue equations. Its first Fourier coefficient is $1$. Since each $g_j$ is obtained from $f_j$ by multiplying by a nonzero scalar, the family $g_1,\dots,g_d$ is again a basis of $V$. This proves the theorem.
[/step]