Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Bound state spectra for the Coulomb problem

We now focus on solving the Coulomb problem:

−ℏ22μ∇⃗2ψ(x⃗)−e2∣x⃗∣ψ=(E−Ecom)ψ- \frac{\hbar^2}{2\mu}{\vec\nabla}^2 \psi({\vec x}) - \frac{e^2}{|{\vec x}|} \psi = (E - E_{com}) \psi

We will redefine E−Ecom→EE - E_{com} \to E just to make reading (and writing) easier.

Note that the reduced mass μ\mu is

μ=mempme+mp=me1+memp∼me(1−memp+O(me2mp2))\mu = \frac{m_e m_p}{m_e + m_p} = \frac{m_e}{1 + \frac{m_e}{m_p}} \sim m_e \left(1 - \frac{m_e}{m_p} + \cO(\frac{m_e^2}{m_p^2})\right)

since memp∼5×10−4\frac{m_e}{m_p} \sim 5\times 10^{-4}, the reduced mass in this case is very close to the electron mass.

Reduction to the radial equation

This Hamiltonian is invariant under rotations. It therefore makes sense to write the Schroedinger equation in spherical coordinates. This transformation is standard and it is worth it for you to work it out yourself. The end result is:

−ℏ22μ1r∂2∂r2(rψ)−ℏ22μr2L⃗2ψ−e2rψ=Eψ- \frac{\hbar^2}{2\mu} \frac{1}{r} \frac{\del^2}{\del r^2} (r\psi) - \frac{\hbar^2}{2\mu r^2} {\vec L}^2\psi - \frac{e^2}{r} \psi = E\psi

Now all of the angular derivatives are contained in L⃗2{\vec L}^2. Thus, we can us separation of variables, and write ψ=f(r)χ(θ,ϕ)\psi = f(r) \chi(\theta,\phi). Inserting this equation, and dividing by ψ\psi, we find

−ℏ22μf1r∂2∂r2(rf)−12μr2χL⃗2χ−e2r=E- \frac{\hbar^2}{2\mu f} \frac{1}{r} \frac{\del^2}{\del r^2} (rf) - \frac{1}{2\mu r^2 \chi} {\vec L}^2\chi - \frac{e^2}{r} = E

Further multiplying by 2μr22\mu r^2, there is only one term involving χ\chi or any other function of the angular coordinates, which is the term proportional to L⃗2χ{\vec L}^2\chi. Thus we can deduce that

1χL⃗2χ=C\frac{1}{\chi}{\vec L}^2 \chi = C

for some constant CC. Multiplying this by χ\chi, we have an eigenvalue equation for the operator L⃗2{\vec L}^2. But we know the solution -- χ\chi must be a spherical harmonic Yℓ,m(θ,ϕ)Y_{\ell,m}(\theta,\phi), and C=ℏ2ℓ(ℓ+1)C = \hbar^2\ell(\ell + 1). Insering this into (4), we find:

−ℏ22μ1r∂2∂r2(rf(r))=ℏ2ℓ(ℓ+1)2μr2f−e2rf=Ef- \frac{\hbar^2}{2\mu} \frac{1}{r} \frac{\del^2}{\del r^2}(r f(r)) = \frac{\hbar^2 \ell(\ell + 1)}{2\mu r^2} f - \frac{e^2}{r} f = E f

Finally, if we multiply the entire equation by rr and define u=rf(r)u = r f(r), we have and “effective” Schroedinger equation:

−ℏ22μ∂r2u(r)+[ℏ2ℓ(ℓ+1)2μr2−e2r]u=Eu- \frac{\hbar^2}{2\mu} \del_r^2 u(r) + \left[\frac{\hbar^2 \ell(\ell + 1)}{2\mu r^2} - \frac{e^2}{r}\right] u = E u

This is a one-dimensional Schroedinger equation with an “effective potential”

Veff=ℏ2ℓ(ℓ+1)2μr2−e2rV_{eff} = \frac{\hbar^2 \ell(\ell + 1)}{2\mu r^2} - \frac{e^2}{r}

This is sketched below. There is a Coulomb potential and for ℓ≠0\ell \neq 0, a “centrifugal barrier” keeping the electron from the proton. This barrier is absent for ℓ=0\ell = 0, and one might worry that the problem is ill-defined because the potential is singular. Classically, this is only a problem if we have truly point particles and fire them right at each other. Quantum-mechanically, we will see that the problem is well-defined for the correct choice of boundary conditions. In the end, rr is an operator, and ⟨e2r⟩\vev{\frac{e^2}{r}} in any well-defined quantum state is finite.

Effective potential for the Coulomb problem

It should be clear from this figure that we expect bound states for E<0E < 0, since the wavefunction will be in the classically forbidden region for r→∞r \to \infty and will die off exponentially. for E>0E > 0 we expect a continuum of propagating states, desribing scattering of the electron and proton.

Boundary conditions for u(r)u(r)

Before continuing, we need to examine the boundary conditions for u(r)u(r) at r=0r = 0, r=∞r = \infty. In the latter case, since Veff→0V_{eff} \to 0, the approximate solution is

u(r)∼e±γr ;  γ=2μ∣E∣ℏ2u(r) \sim e^{\pm \gamma r}\ ; \ \ \gamma = \sqrt{\frac{2\mu |E|}{\hbar^2}}

Clearly we need to choose the eponentially decaying solution (or the wavefunction is not normalizable).

For r→0r \to 0, the story is a little more subtle. Let us first consider ℓ>0\ell > 0. As r→0r \to 0, the centrifugal barrier term in the Schroedinger equation dominates EuE u or e2u/re^2 u/r, because 1/r21/r^2 blows up more quackly that 1/r1/r. Thus, as r→0r \to 0, the Schrodinger equation becomes

−ℏ22μ∂r2u+ℏ2ℓ2μr2u=0- \frac{\hbar^2}{2\mu} \del_r^2 u + \frac{\hbar^2 \ell}{2\mu r^2} u = 0

multiplying through by 2μ/ℏ22\mu/\hbar^2 we can see that both terms scale the same as r→λrr \to \lambda r; that is, they both scale as λ−2\lambda^{-2}. In this case, a solution of the form u∝rαu \propto r^{\alpha} makes sense. since both terms will reduce the power of α\alpha and can be balanced against each other. The result is

−α(α−1)+ℓ(ℓ+1)=0- \alpha(\alpha - 1) + \ell(\ell + 1) = 0

which has solutions α=−ℓ,ℓ+1\alpha = - \ell, \ell +1. The former leads to a wavefunction of the form

ψ=urYℓ,m∼r−ℓ−1\psi = \frac{u}{r} Y_{\ell,m} \sim r^{-\ell - 1}

When we compute the norm,

∫drr2dϕsin⁡θdθ∣ψ∣2∼∫dr1r2ℓ+1+…\int dr r^2 d\phi \sin\theta d\theta |\psi|^2 \sim \int dr \frac{1}{r^{2\ell + 1}} + \ldots

which clearly diverges; so we demand that u→rℓ+1u\to r^{\ell + 1} as r→0r \to 0.

When ℓ=0\ell = 0, the Coulomb term dominates, but scales differently in rr. We assume that uu can be expanded in a power series, u∼rα(1+a1t+a2r2+…)u \sim r^{\alpha} (1 + a_1 t + a_2 r^2 + \ldots). Inserting this into the equation

−ℏ22μ∂r2u+eru=Eu- \frac{\hbar^2}{2\mu} \del_r^2 u + \frac{e}{r} u = E u

the most singular term comes from the second derivative acting on rαr^{\alpha}, and is proportional to α(α−1)rα−2\alpha(\alpha - 1)r^{\alpha - 2}. This has the solutions u∝1u \propto 1 and u∝ru \propto r. In fact both are annihilated by the second derivative, but when we expand the solution to the full equation in a power series, we will see that it comes out OK.

Now u∼1⇒ψ∼1ru \sim 1 \Rightarrow \psi \sim \frac{1}{r}. While this is putatively ormalizable at the origin, it is in fact problematic; the wvafunction has a sharp kin at the origin. One can in fact show that ∇⃗21r∝δ(r){\vec \nabla}^2 \frac{1}{r} \propto \delta(r). This delta function singularity cannot balance against the other terms in the Schroedinger equation. Thus, we demand that u∼ru \sim r.

Power series solution

Before continuing, we will follow our treatment of the simple harmonic oscillator and try to make the Schrodinger equation one in terms of dimensionless variables. First, we note that for an energy EE, as r→∞r \to \infty the equation becomes

∂r2u=2μ∣E∣ℏ2u\del_r^2 u = \frac{2\mu |E|}{\hbar^2} u

with the solution

u∼e−2μ∣E∣ℏ2ru \sim e^{- \sqrt{\frac{2\mu |E|}{\hbar^2}} r}

This dies off with a characteristic length scale L=ℏ22μ∣E∣L = \sqrt{\frac{\hbar^2}{2\mu |E|}}. Defining r=Lρr = L \rho, the radial Schroedinger equation becomes:

−∂ρ2u+ℓ(ℓ+1)u−2μℏ2∣E∣e2ρu=u- \del_{\rho}^2 u + \frac{\ell(\ell + 1)} u - \sqrt{\frac{2\mu}{\hbar^2|E|}}\frac{e^2}{\rho} u = u

We further define λ=2μℏ2∣E∣\lambda = \sqrt{\frac{2\mu}{\hbar^2 |E|}}. Next, we strip off the large-ρ\rho behavior by writing u=e−ρG(ρ)u = e^{-\rho} G(\rho). The resulting equation is

ddρ2G−2ddρG+[e2λρ−ℓ(ℓ+1)ρ2]G=0\frac{d}{d\rho^2}G - 2 \frac{d}{d\rho} G + \left[\frac{e^2\lambda}{\rho} - \frac{\ell(\ell + 1)}{\rho^2}\right]G = 0

These terms scale as either σ−2\sigma^{-2} or σ−1\sigma^{-1} as ρ→σρ\rho \to \sigma \rho. If we expand GG in a power series, these terms will relate successive terms, leaving the possibility that the coefficients are at least determined by a recursion relation. From the discussion above, we try

G(ρ)=ρ∑n=0∞cnρnG(\rho) = \rho \sum_{n = 0}^{\infty} c_n \rho^n

Inserting this into Eq. (18), and collecting all terms of the same order in ρ\rho, we find:

∑m=1∞[cm(m+ℓ+1)(m+ℓ)−ℓ(ℓ+1)cm−(2(m+ℓ)−e2λ)cm−1]ρm−1\sum_{m = 1}^{\infty} \left[c_m(m + \ell + 1)(m + \ell) - \ell(\ell + 1) c_m - (2(m + \ell) - e^2\lambda) c_{m-1}\right]\rho^{m - 1}

This is solved if

cmcm−1=2(m+ℓ)−e2λ(m+ℓ+1)(m+ℓ)−ℓ(ℓ+1)\frac{c_m}{c_{m-1}} = \frac{2(m + \ell) - e^2 \lambda}{(m + \ell + 1)(m +\ell) - \ell(\ell + 1)}

Thus each cmc_m is determined by the previous cm−1c_{m-1}. This leaves c0c_0 undetermined, but it will be fixed by normalization of the wavefunction.

Unfortunately, this power series will in general descibe a diverging wavefunction. So see this note that for large mm, cm/cm−1∼2/mc_m/c_{m-1} \sim 2/m, so that cm∼2m/m!c_m \sim 2^m/m!. Summing

ρ∑mcmρm∼ρ∑m(2ρ)mm!∼ρe2ρ\rho \sum_m c_m \rho^m \sim \rho \sum_m \frac{(2 \rho)^m}{m!} \sim \rho e^{2\rho}

This overwhelms the e−ρe^{-\rho} behavior and renders the wavefunction non-normalizable.

The out is if the power series truncates at some order. Then the wavfunction scales as ρme−ρ\rho^m e^{-\rho} which is normalizable as ρ→∞\rho \to \infty. This means that there is some m=1,2,3,4...m = 1, 2, 3, 4... such that e2λ=2(m+ℓ)e^2 \lambda = 2(m + \ell). In other words,

E=−μe42ℏ2n2E = - \frac{\mu e^4}{2\hbar^2 n^2}

where n=1,2,3,…=(m+ℓ)n = 1,2,3,\ldots = (m + \ell). This means that

ρ−μe4ℏ2nr\rho - \frac{\mu e^4}{\hbar^2 n} r

so that the charateristic scale of exponential falloff is Ln∼ℏ2nμe4L_n \sim \frac{\hbar^2 n}{\mu e^4}. We call L1=ℏ2μe4=a0L_1 = \frac{\hbar^2}{\mu e^4} = a_0 the *Bohr radius$, defining the characteristic side of the ground state.

At the end of the day, the wavefunction is a polynomial times an exponential. The polynomials turn out to be well-known functions. We can write the full solution as

ψ(r,θ,ϕ)∝Yℓ,m(θ,ϕ)e−r/(na0)(rna0)ℓLn−ℓ−12ℓ+1(2rna0)\psi(r,\theta,\phi) \propto Y_{\ell,m}(\theta,\phi) e^{-r/(n a_0)} \left(\frac{r}{n a_0}\right)^{\ell} L_{n-\ell-1}^{2\ell + 1}\left(\frac{2 r}{n a_0}\right)

where LnmL_n^m are the associated Laguerre polynomials whose properties can be looked up in a good book on special functions. Note that the characteristic scale is na0n a_0.

If we replace the proton with a nucleus with charge ZZ, the Coulomb potential gets replaced by e2/r→Ze2/re^2/r \to Z e^2/r. Then a0→a0/Z2a_0 \to a_0/Z^2, and En→Z2EnE_n \to Z^2 E_n.

The spectrum

As we stated, the bound states have energy (23). Looking back, we see that n=ℓ+mn = \ell + m where m≥1m \geq 1. Thus, for any value of mm, we have the possible values ℓ=n−1,n−2,…0\ell = n - 1, n -2, \ldots 0. For n>1n > 1 this is a highly degenerate spectrum because for each value of ℓ\ell there is a 2ℓ+12\ell + 1-dimensional subspace. This degeneracy is the result of an additional symmetry, that we will discuss below.

It is standard to refer to ℓ=0\ell = 0 states by SS; ℓ=1\ell = 1 states by PP, ℓ=2\ell = 2 states by DD, and so on.

Finally, note that the low-lying states are separated by energies of order a fraction of the magnitude of the ground state energy. As n→∞n \to \infty the spacing becomes finer and finer, as the energy approaches zero and the wavefunction can spread out more and more.

Explaining the degeneracy

Whenever there is a large degeneracy, the lore is that there is some meaningful symmetry group for which the degenerate subspaces form irreps. “Meaningful” is what makes this statement interesting. Given a degeneracy one could always design a set of operators commuting with the Hamiltonian for which this is true. But usually it comes from some interesting feature of the dynamics.

(More generally, what happens in the real world is that we have some approximate degeneracies, as seen below:

Nearly degenerate spectrum

In this case, when δϵ≪ΔE\delta\eps \ll \Delta E, we can often explain this by an “approximate symmetry”, one in which the Hamiltonian is of the form H0+H1H_0 + H_1 where H0H_0 is degenerate with level spacing ΔE\Delta E, and has a symmetry that forces this degeneracy; while H1H_1 breaks this symmetry but has O(δϵ)\cO(\delta\eps) matrix elements. We will see examples of this soon.)

As it happens, in classical mechanics, there is a a conserved vectorial quantity called the Runge-Lenz vector:

R⃗=1μp⃗×L⃗−e2r⃗r{\vec R} = \frac{1}{\mu} {\vec p}\times{\vec L} - \frac{e^2 {\vec r}}{r}

The classical Coulomb problem, being mathematically identical to the Kepler problem, supports elliptical orbits. The Runge-Lenz vector points from the force center (for us, r=0r = 0) to the point of closest approach. This is true anywhere along the orbit, thus it is conserved. Note that it is only conserved for inverse-square law forces; if we change to some 1/rα1/r^{\alpha} central force with α≠2\alpha \neq 2, we lose that conservation law.

Deriving the degeeracy for this takes several steps. First, a highly nontrivial calculation yields:

R⃗2=e4+2H(L⃗2+ℏ2)μ{\vec R}^2 = e^4 + \frac{2H({\vec L}^2 + \hbar^2)}{\mu}

Next, wf we promote R⃗{\vec R} to an operator acting on the Hilbert space of the Hydrogen atom, we find that

[Ri,H]=0[Ri,Lj]=iℏϵijkRk[Ri,Rj]=−2iℏHμϵijkLk\begin{align} [R^i, H] & = 0\\ [R^i, L_j] & = i\hbar \epsilon_{ijk} R^k\\ [R^i, R^j] & = - \frac{2i\hbar H}{\mu} \epsilon^{ijk}L_k \end{align}

We consider only E<0E < 0 states, and redefine

K⃗=−μ2HR⃗{\vec K} = \sqrt{\frac{ - \mu}{2H}} {\vec R}

with this, we can calculate:

H=−μe22(K⃗2+L⃗2+ℏ2)H = - \frac{\mu e^2}{2({\vec K}^2 + {\vec L}^2 + \hbar^2)}

Next, if we define

M⃗=L⃗+K⃗N⃗=L⃗−K⃗\begin{align} {\vec M} & = {\vec L} + {\vec K} \\ {\vec N} & = {\vec L} - {\vec K} \end{align}

Here MiM^i commutes with any component of N⃗{\vec N}, while

[Mi,Mj]=iℏϵijkMk[Ni,Nj]=iℏϵ6ijkNk\begin{align} [M^i, M^j] & = i\hbar \eps^{ijk} M^k\\ [N^i, N^j] & = i\hbar\eps6{ijk} N^k \end{align}

This is just two SU(2)SU(2) algebras. (Note that SO(4)SO(4), the space of rotations in 4d, is isomorphic to SU(2)×SU(2)SU(2)\times SU(2). However, their casimirs are not distinct. Since L⃗{\vec L} is always perpendicular to the ellipse, L⃗⋅r⃗=0{\vec L}\cdot{\vec r} = 0; and L⃗⋅(p⃗×L⃗)=0{\vec L}\cdot({\vec p}\times{\vec L}) = 0 automatically, we have K⃗dotL⃗=0{\vec K}dot{\vec L} = 0. This means that M⃗2=N⃗2{\vec M}^2 = {\vec N}^2 = K⃗2+L⃗2{\vec K}^2 + {\vec L}^2. Finally, we can show that

H=−μe42(4M⃗2+ℏ2)H = - \frac{\mu e^4}{2(4{\vec M}^2 + \hbar^2)}

Now because it is part of an SU(2)SU(2) algebra, M⃗2{\vec M}^2 has eigenvalues ℏ2m(m+1)\hbar^2 m(m+1), and the corresponding irreps with this eigenvalue have degeneracy 2m+12m + 1. In this case N⃗2=ℏ2(2m+1){\vec N}^2 = \hbar^2(2m+1) also, so the total degeneracy is (2m+1)2(2m + 1)^2. The energy eigenvalues are:

H=−μe42ℏ2(4m(m+1)+1)=−μe42ℏ2(2m+1)2H = - \frac{\mu e^4}{2\hbar^2(4m(m+1) + 1)} = - \frac{\mu e^4}{2\hbar^2(2m + 1)^2}

Now if we let m=0,1/2,1,3/2,…m = 0,1/2,1,3/2,\ldots, then 2m+1=1,2,3,4,…2m + 1 = 1, 2, 3, 4,\ldots> In other words we can identify n=(2m+1)n = (2m + 1). We have thus explained the n2n^2 degeneracy of the eigenstates of the Coulomb problem.