[proofplan]
We reduce the degeneracy computation to the separated-variable Coulomb spectral classification included in the theorem statement for the standard self-adjoint spinless Coulomb Hamiltonian on $L^2(\mathbb{R}^3)$. That classification supplies the completeness of separated eigenfunctions, the allowed angular labels $l=0,\dots,n-1$, the magnetic labels $m=-l,\dots,l$, and the one-dimensionality of the regular square-integrable radial solution space for each pair $(n,l)$. Once those spectral inputs are in place, the proof is a finite-dimensional count: each fixed $l$ contributes $2l+1$ states, and summing over $l$ gives $n^2$.
[/proofplan]
[step:Decompose the eigenspace into angular momentum sectors]
Let $\mathbb{N}=\{1,2,3,\dots\}$ denote the positive integers, and let $\mathcal{L}^1$ denote one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) on $\mathbb{R}$. Fix $n \in \mathbb{N}$. Let $H$ denote the standard self-adjoint spinless Coulomb Hamiltonian on $L^2(\mathbb{R}^3)$, and let $I:L^2(\mathbb{R}^3) \to L^2(\mathbb{R}^3)$ denote the identity operator. Define the eigenspace at energy $E_n$ by $\mathcal{E}_n := \ker(H - E_n I) \subset L^2(\mathbb{R}^3)$.
We now use the separated-variable Coulomb spectral classification stated in the theorem hypothesis. Its hypotheses are met because $H$ is the standard self-adjoint spinless Coulomb Hamiltonian on $L^2(\mathbb{R}^3)$ and $E_n$ is one of its bound-state energies, exactly as required in that classification. The classification states that every element of $\mathcal{E}_n$ is a finite linear combination of separated eigenfunctions with angular momentum label $l \in \{0,\dots,n-1\}$ and magnetic label $m \in \{-l,\dots,l\}$, and conversely that every such allowed pair gives an eigenfunction in $\mathcal{E}_n$ after multiplication by the corresponding regular square-integrable radial factor.
For each $l \in \{0,\dots,n-1\}$, let $\mathcal{Y}_l \subset L^2(S^2,\mathcal{H}^2)$ denote the spherical harmonic subspace of angular momentum $l$, where $S^2 \subset \mathbb{R}^3$ is the unit sphere and $\mathcal{H}^2$ is surface measure on $S^2$. Let $\mu_r$ denote the measure on $(0,\infty)$ defined by $d\mu_r(r)=r^2\,d\mathcal{L}^1(r)$. Let $\mathcal{R}_{n,l} \subset L^2((0,\infty),\mu_r)$ denote the [vector space](/page/Vector%20Space) of regular square-integrable functions $R:(0,\infty)\to\mathbb{C}$ solving the radial Coulomb eigenvalue equation
\begin{align*}
-a\left(R''(r)+\frac{2}{r}R'(r)-\frac{l(l+1)}{r^2}R(r)\right)-\frac{b}{r}R(r)=E_nR(r)
\end{align*}
for $r \in (0,\infty)$, where $a,b>0$ are the constants in $H=-a\Delta-b/|x|$. Let $U_{\mathrm{sph}}$ denote the spherical-coordinate Hilbert-space identification from $L^2(\mathbb{R}^3)$ to $L^2((0,\infty),\mu_r) \widehat{\otimes} L^2(S^2,\mathcal{H}^2)$ obtained from $x=r\omega$ with $r \in (0,\infty)$ and $\omega \in S^2$. The classification gives the vector-space direct sum decomposition, induced through $U_{\mathrm{sph}}$,
\begin{align*}
\mathcal{E}_n = \bigoplus_{l=0}^{n-1} \mathcal{R}_{n,l} \otimes \mathcal{Y}_l.
\end{align*}
The sum is direct because distinct spherical harmonic sectors are orthogonal in $L^2(S^2,\mathcal{H}^2)$, and within each fixed pair $(n,l)$ the classification supplies exactly one independent regular square-integrable radial solution.
[guided]
Let $\mathbb{N}=\{1,2,3,\dots\}$ denote the positive integers, and let $\mathcal{L}^1$ denote one-dimensional Lebesgue measure on $\mathbb{R}$. Fix $n \in \mathbb{N}$. Let $H$ denote the standard self-adjoint spinless Coulomb Hamiltonian on $L^2(\mathbb{R}^3)$, and let $I:L^2(\mathbb{R}^3) \to L^2(\mathbb{R}^3)$ denote the identity operator. The eigenspace whose dimension we want to compute is $\mathcal{E}_n := \ker(H - E_n I) \subset L^2(\mathbb{R}^3)$.
The structural input is the separated-variable Coulomb spectral classification stated in the theorem hypothesis. We use it as a stated hypothesis, not as a result proved inside this counting argument. Its hypotheses are satisfied here because $H$ is the standard self-adjoint spinless Coulomb Hamiltonian on $L^2(\mathbb{R}^3)$ and $E_n$ is the bound-state energy level under consideration. The classification supplies four pieces of information: the separated eigenfunctions are complete in $\mathcal{E}_n$, the allowed angular momentum labels are exactly $l \in \{0,\dots,n-1\}$, the corresponding magnetic labels are exactly $m \in \{-l,\dots,l\}$, and for each fixed pair $(n,l)$ the regular square-integrable radial solution space is one-dimensional.
For each $l \in \{0,\dots,n-1\}$, let $\mathcal{Y}_l \subset L^2(S^2,\mathcal{H}^2)$ be the spherical harmonic subspace of angular momentum $l$, where $S^2 \subset \mathbb{R}^3$ is the unit sphere and $\mathcal{H}^2$ is surface measure on $S^2$. Let $\mu_r$ denote the measure on $(0,\infty)$ defined by $d\mu_r(r)=r^2\,d\mathcal{L}^1(r)$. Let $\mathcal{R}_{n,l} \subset L^2((0,\infty),\mu_r)$ be the space of regular square-integrable functions $R:(0,\infty)\to\mathbb{C}$ solving the radial Coulomb eigenvalue equation
\begin{align*}
-a\left(R''(r)+\frac{2}{r}R'(r)-\frac{l(l+1)}{r^2}R(r)\right)-\frac{b}{r}R(r)=E_nR(r)
\end{align*}
for $r \in (0,\infty)$, where $a,b>0$ are the constants in $H=-a\Delta-b/|x|$. Let $U_{\mathrm{sph}}$ denote the spherical-coordinate Hilbert-space identification from $L^2(\mathbb{R}^3)$ to $L^2((0,\infty),\mu_r) \widehat{\otimes} L^2(S^2,\mathcal{H}^2)$ obtained from $x=r\omega$ with $r \in (0,\infty)$ and $\omega \in S^2$. With this notation, the classification identifies the eigenspace as the direct sum
\begin{align*}
\mathcal{E}_n = \bigoplus_{l=0}^{n-1} \mathcal{R}_{n,l} \otimes \mathcal{Y}_l.
\end{align*}
This tensor-product notation records the separated-variable construction through $U_{\mathrm{sph}}$: a radial solution is paired with a spherical harmonic. The directness is part of the same classification and is also reflected by angular orthogonality in $L^2(S^2,\mathcal{H}^2)$: different angular momentum sectors do not overlap. Therefore the dimension of $\mathcal{E}_n$ is computed by adding the dimensions of the finite-dimensional summands $\mathcal{R}_{n,l} \otimes \mathcal{Y}_l$.
[/guided]
[/step]
[step:Count the contribution from each fixed angular momentum]
Fix $l \in \{0,\dots,n-1\}$. By the radial conclusion of the separated-variable Coulomb bound-state classification assumed in the theorem statement, the space $\mathcal{R}_{n,l}$ of regular square-integrable radial solutions in $L^2((0,\infty),\mu_r)$ is one-dimensional:
\begin{align*}
\dim \mathcal{R}_{n,l} = 1.
\end{align*}
The spherical harmonic sector of degree $l$ has magnetic labels $m=-l,-l+1,\dots,l$, hence
\begin{align*}
\dim \mathcal{Y}_l = |\{-l,-l+1,\dots,l\}| = 2l+1.
\end{align*}
Therefore the dimension of the $l$-th summand is
\begin{align*}
\dim(\mathcal{R}_{n,l} \otimes \mathcal{Y}_l) = \dim \mathcal{R}_{n,l} \cdot \dim \mathcal{Y}_l = 1 \cdot (2l+1) = 2l+1.
\end{align*}
[/step]
[step:Sum the sector dimensions and evaluate the arithmetic series]
Since $\mathcal{E}_n$ is the direct sum of the sectors $\mathcal{R}_{n,l} \otimes \mathcal{Y}_l$ for $l=0,\dots,n-1$, finite-dimensional additivity of dimension gives
\begin{align*}
\dim \mathcal{E}_n = \sum_{l=0}^{n-1} \dim(\mathcal{R}_{n,l} \otimes \mathcal{Y}_l).
\end{align*}
Using the computation from the previous step,
\begin{align*}
\dim \mathcal{E}_n = \sum_{l=0}^{n-1}(2l+1).
\end{align*}
We evaluate the finite sum:
\begin{align*}
\sum_{l=0}^{n-1}(2l+1) = 2\sum_{l=0}^{n-1}l + \sum_{l=0}^{n-1}1.
\end{align*}
The two sums are
\begin{align*}
\sum_{l=0}^{n-1}l = \frac{n(n-1)}{2}
\end{align*}
and
\begin{align*}
\sum_{l=0}^{n-1}1 = n.
\end{align*}
Substituting these values gives
\begin{align*}
\sum_{l=0}^{n-1}(2l+1) = 2 \cdot \frac{n(n-1)}{2} + n = n(n-1)+n = n^2.
\end{align*}
Hence
\begin{align*}
\dim \mathcal{E}_n = n^2.
\end{align*}
This is precisely the spinless hydrogenic degeneracy of the Coulomb energy level $E_n$.
[guided]
The direct sum decomposition reduces the degeneracy computation to ordinary finite-dimensional linear algebra. Since
\begin{align*}
\mathcal{E}_n = \bigoplus_{l=0}^{n-1} \mathcal{R}_{n,l} \otimes \mathcal{Y}_l,
\end{align*}
finite-dimensional additivity of dimension gives
\begin{align*}
\dim \mathcal{E}_n = \sum_{l=0}^{n-1} \dim(\mathcal{R}_{n,l} \otimes \mathcal{Y}_l).
\end{align*}
For each allowed angular momentum label $l$, the radial solution space has dimension $1$ and the spherical harmonic space has dimension $2l+1$. Therefore
\begin{align*}
\dim(\mathcal{R}_{n,l} \otimes \mathcal{Y}_l) = 1 \cdot (2l+1) = 2l+1.
\end{align*}
Substituting this sector count gives
\begin{align*}
\dim \mathcal{E}_n = \sum_{l=0}^{n-1}(2l+1).
\end{align*}
Now evaluate the arithmetic series by splitting it into the sum of even parts and the sum of ones:
\begin{align*}
\sum_{l=0}^{n-1}(2l+1) = 2\sum_{l=0}^{n-1}l + \sum_{l=0}^{n-1}1.
\end{align*}
The finite sums are
\begin{align*}
\sum_{l=0}^{n-1}l = \frac{n(n-1)}{2}
\end{align*}
and
\begin{align*}
\sum_{l=0}^{n-1}1 = n.
\end{align*}
Therefore
\begin{align*}
\sum_{l=0}^{n-1}(2l+1) = 2 \cdot \frac{n(n-1)}{2} + n = n(n-1)+n = n^2.
\end{align*}
Hence $\dim \mathcal{E}_n = n^2$, which is the claimed spinless degeneracy of the energy level $E_n$.
[/guided]
[/step]