# Modern theory of the hydrogen atom

(Part of the Wolverhampton Lectures of Physics's Quantum Physics Course)

However accurate, Bohr's model was however clearly incomplete and handwaving. In particular, although it provides correctly the energy spectrum of Hydrogen, this was for this particular case only, for which it works exceedingly well, but even there, it finds quickly its limitation, for instance, it does not say anything about the details of the transitions between the level (the rate at which they occur and the underlying physical mechanism).

The full picture is given by solving Schrödinger's equation for the wavefunction of the electron:

$$H\psi=E\psi$$

where the Hamiltonian is the electron's kinetic energy in the Coulomb potential

$$H=-{\hbar^2\over 2m}\nabla^2-{e^2\over4\pi\epsilon_0}{1\over r}$$

Note that doing so, we do not describe the proton's dynamics, which only provides the potential. This assumes that the proton has an infinite mass as compared to the electron, which can be relaxed by working in the reference frame not of the proton itself but of the center of mass of the compound system, in which case, $m\rightarrow mM/(m+M)$ where $M$ is the mass of the proton.

Since $H$ is spherically symmetric, the problem is solved for the wavefunction $\psi(r,\theta,\phi)=R(r)Y(\theta,\phi)$, with $Y$ a spherical harmonic and the corresponding radial equation for $u(r)=rR(r)$:

$$-{\hbar^2\over2m}{d^2u\over dr^2}+\left[-{e^2\over4\pi\epsilon_0}{1\over r}+{\hbar^2\over 2m}{l(l+1)\over r^2}\right]u=Eu\,.$$

Introducing the usual $k$ number, this time for negative energies $E<$ as the problem will have both bound and scattering states:

$$\kappa={\sqrt{-2mE}\over\hbar}$$

and also introducing a rescaled radius:

$$\rho\equiv\kappa r\quad\text{where}\quad\rho_0\equiv{me^2\over2\pi\epsilon_0\hbar^2\kappa}$$

$${d^2u\over d\rho^2}=\left[1-{\rho_0\over\rho}+{l(l+1)\over\rho^2}\right]u\,.\tag{1}$$

To solve this equation, we consider first the asymptotics:

• At $\rho\rightarrow\infty$, the constant term dominates so that the equation reads ${d^2u\over d\rho^2}=u$ with general solution $u(\rho)=A\exp(-\rho)+B\exp(\rho)$, which means, as the wavefunction should be normalizable, that $B=0$ and

$$u(\rho)\rightarrow A\exp(-\rho)\,.$$

• On the other hand, as $\rho\rightarrow0$, the centrifugal term dominates and the equation tends to

$${d^2u\over d\rho^2}={l(l+1)\over\rho^2}u\,.$$

A solution in this case is $u(\rho)=C\rho^{l+1}+D\rho^{-l}$, which one can check easily since $u'(\rho)=C(l+1)\rho^{l}-Dl\rho^{-l-1}$ and, iterating, $u''(\rho)=Cl(l+1)\rho^{l-1}+Dl(l+1)\rho^{-l-2}={l(l+1)\over\rho^2}u$. Having two independent terms for a second-order linear differential equation, this is actually the most general solution. This time the term $\rho^{-l}$ blows up (as $\rho\rightarrow0$) so that in this limit, the solution reads:

$$u(\rho)\rightarrow C\rho^{l+1}\,.$$

We now introduce still another variable $v(\rho)$, to "peel off" the asymptotic behaviour:

$$u(\rho)=\rho^{l+1}e^{-\rho}v(\rho)\,.$$

If we now compute the successive derivatives, we find:

\begin{align} {du\over d\rho}&=\rho^le^{-\rho}\left[(l+1-\rho)v+\rho{du\over d\rho}\right]\,,\\ {d^2u\over d\rho^2}&=\rho^le^{-\rho}\left\{\left[-2l-2+\rho+{l(l+1)\over\rho}\right]v+2(l+1-\rho){dv\over d\rho}+\rho{d^2v\over d\rho^2}\right\}\,, \end{align}

so that, inserting this into the radial equation (1), we find the equation in terms of $v$:

$$\rho\frac{d^2v}{d\rho^2}+2(l+1-\rho)\frac{dv}{d\rho}+[\rho_0-2(l+1)]v=0\,.\tag{2}$$

The reason to bring it in this form is that we can now solve it with a series expansion for $v$:

$$v(\rho)=\sum_{j=0}^\infty a_j\rho^j$$

for which, also computing the respective derivatives:

\begin{align} \frac{dv}{d\rho}&=\sum_{j=0}^\infty ja_j\rho^{j-1}=\sum_{j=0}^\infty(j+1)a_{j+1}\rho^j\,,\\ {d^2v\over d\rho^2}&=\sum_{j=0}^\infty j(j+1)a_{j+1}\rho^{j-1}\,, \end{align}

which, once inserted into Eq. (2), yields

$$\sum_{j=0}^\infty\left\{j(j+1)a_{j+1}\rho^j+2(l+1)(j+1)a_{j+1}\rho^j-2ja_j\rho^j+[\rho_0-2(l+1)]a_j\rho^j\right\}=0$$

Equating coefficients with the same power we find a recursive relation:

$$a_{j+1}=\frac{2(j+l+1)-\rho_0}{(j+1)(j+2l+2)}a_j\,.\tag{3}$$

Let us consider again the asymptotic behaviour. For $j$ large enough, Eq. (3) becomes $a_{j+1}\approx(2/j)a_j$ with solution, neglecting the departures for small~$j$

$$a_j={2^j\over j!}a_0\tag{4}$$

in which case, the solution would be $v(\rho)=a_0\sum_{j=0}^\infty(2\rho)^j/j!=a_0e^{2\rho}$, but since $u(\rho)=A\rho^{l+1}e^{-\rho}v(\rho)$, we have again an exponential divergence as $\rho\rightarrow\infty$. At such, the solution is also not physical. Unless... the series is finite, in which case the asymptotic behaviour of Eq. (4) can be excluded. This therefore brings a new constrain:

$$2(j+l+1)-\rho_0=0\,.$$

This brings another integer, the principal quantum number $n$ (already introduced by Bohr):

$$n\equiv j_\mathrm{max}+l+1$$

for $j_\mathrm{max}$ the integer in the series expansion that cancels the numerator of the recursive relation, at which point:

$$\rho_0=2n\,.$$

In this way, we have linked $n$ to the energy, since $\rho_0$ was itself energy dependent. This yields:

$$E_n=-{m\over2\hbar^2}\left({e^2\over4\pi\epsilon_0}\right)^2{1\over n^2}$$

and therefore:

$$E_n={E_1\over n^2}$$

where $E_1$ is as defined in the old Bohr theory, but it has now be derived as part of the Schrödinger equation. This also provides a sound derivation for the Bohr radius $a_\mathrm{B}$ which is linked to $\rho$ as

$$\rho={r\over a_\mathrm{B}n}\,.$$

We can now work out recursively the wavefunctions for the various quantum numbers by plugging back $v(\rho)$ in the radial wavefunction

$$R_{nl}(r)={1\over r}\rho^{l+1}e^{-\rho}v(\rho)\,.$$

for the various quantum numbers $n=1, 2, \ldots$ with $j_\mathrm{max}=n-l-1$ for all the possible positive integer values of $l$. In terms of $l$, since $l=n-1-j_\mathrm{max}$, one can see that one term of the sum is always a solution, when $l=n-1$, in which case $j_\mathrm{max}=0$ and the series is a single term. Since $l\ge0$, $j_\mathrm{max}$ is at most $n-1$, providing $n$ terms for the $l=0$ wavefunction. We thus have $n$ values of $l$, which we can all work out in turns:

• The ground state comes from $n=1$, in which case $l=0$ (as otherwise $j$ would be negative) in which case, the recursive relation yields $a_1=0$. The series truncates at the first term and $v(\rho)=a_0$ and thus $R_{10}(r)={1\over r}\rho e^{-\rho}a_0$ or, expressing all terms as function of $r$

$${\displaystyle R_{10}(r)={a_0\over a_\mathrm{B}}\exp(-r/a_\mathrm{B})}\,.$$

The normalization yields

$$\int_0^\infty|R_{10}|^2r^2\,dr={a_0^2\over a_\mathrm{B}^2}\int_0^\infty e^{-2r/a_\mathrm{B}}r^2\,dr={a_0^2a_\mathrm{B}^2\over4}=1$$

so that $a_0=2/\sqrt{a_\mathrm{B}}$ and, since $Y_0^0=1/\sqrt{4\pi}$ from the general theory, we find the complete expression for the ground state:

$$\boxed{\displaystyle\psi_{100}(r,\theta,\phi)={1\over\sqrt{\pi a_\mathrm{B}^3}}\exp(-r/a_\mathrm{B})\,.}$$

• The first excited state has $n=2$, so that $j=1-l$ making $l=0$ or $1$.
We start with $l=0$, for which the recursive relation yields:

$$a_1=-a_0\quad\text{and}\quad a_2=0$$

so that $v(\rho)=a_0(1-\rho)$ and thus:

$$R_{20}(r)={a_0\over2a_\mathrm{B}}\left(1-{r\over 2a_\mathrm{B}}\right)e^{-r/2a_\mathrm{B}}\,.$$

Now for $l=1$, the series also self-truncates at the first term, so $v(\rho)$ is constant again and

$$R_{21}(r)={a_0\over4_\mathrm{B}^2}re^{-r/(2a_\mathrm{B})}\,.$$

For each $l$ there is a corresponding number of $m$ solutions. In the general case, from $j=n-l-1$, we find that for each $n$, there are $n$ possible values of $l$, namely:

$$\boxed{l=0, 1, 2, \ldots, n-1}$$

And since for each $l$ we have $(2l+1)$ values of $m$, the degeneracy $d(l)$ (i.e., the number of solutions) of the energy level $E_n$ is

$$d(l)=\sum_{l=0}^{n-1}(2l+1)=n^2\,.$$

As for the series, it turns out that it can be explicited in closed-form in terms of the Laguerre polynomials $L_q$ defined as

$$L_q(x)\equiv e^x\left({d\over dx}\right)^q(e^{-x}x^q)\,.$$

From this, one can define the so-called associated Laguerre polynomials

$$L_{q-p}^p(x)\equiv(-1)^p\left({d\over dx}\right)^pL_q(x)\,,$$

in which case, the function $v(\rho)$ is given by:

$$v(\rho)=L^{2l+1}_{n-l-1}(2\rho)\,.$$

The first few radial wavefunctions can be worked out to be:

\begin{align} R_{10}(r)&=2a_\mathrm{B}^{-3/2}\exp(-r/a_\mathrm{B})\\ R_{20}(r)&={1\over\sqrt{2}}a_\mathrm{B}^{-3/2}\left(1-{1\over2}{r\over a_\mathrm{B}}\right)\exp(-r/(2a_\mathrm{B}))\\ R_{21}(r)&={1\over\sqrt{24}}a_\mathrm{B}^{-3/2}{r\over a_\mathrm{B}}\exp(-r/(2a_\mathrm{B})) \end{align}

The normalization can be found (we are not deriving it but give the result directly):

$$\boxed{\psi_{nlm}(r,\theta,\phi)=\sqrt{\left({2\over na}\right)^3{(n-l-1)!\over2n(n+l)!^3}}\exp\left(-{r\over na_\mathrm{B}}\right)\left({2r\over na_\mathrm{B}}\right)^lL_{n-l-1}^{2l+1}\left({2r\over na_\mathrm{B}}\right)Y_l^m(\theta,\phi)\,.}$$

These functions thus defined form an orthonormal set.

The visualisation of these wavefunctions is a famous problem (see for instance this Wikipedia page). A great piece of code is provided by Jacopo Bertolotti:

\[Alpha]0 = 1;
\[Psi][n_, l_, m_, r_, \[Theta]_, \[Phi]_] := Sqrt[(2/(n \[Alpha]0))^3 (n - l - 1)!/(2 n ((n + l)!))] E^(-r/(n \[Alpha]0)) ((2 r)/(n \[Alpha]0))^l LaguerreL[n - l - 1, 2 l + 1, (2 r)/(n \[Alpha]0)] SphericalHarmonicY[l, m, \[Theta], \[Phi]];
p1 = Flatten@Table[
f = TransformedField["Spherical" -> "Cartesian", \[Psi][n, l, m, r, \[Theta], \[Phi]], {r, \[Theta], \[Phi]} -> {x, y, z}];
DensityPlot3D[Abs[f]^2
, {x, -30, 30}, {y, -30, 30}, {z, -30, 30}, ColorFunction -> Hue, ColorFunctionScaling -> True, Boxed -> False, Axes -> False, PlotLabel -> Style[StringForm["Hydrogen atom orbitals\n |\[Psi]\!$$\*SuperscriptBox[\(|$$, \$$2$$]\) : n= l= m=", n, l, m], Medium, FontFamily -> "DejaVu Serif"], LabelStyle -> {Black, Bold}, RegionFunction -> Function[{x, y, z}, x < 0 || y > 0], PlotLegends -> Automatic], {n, 1, 4}, {l, 0, n - 1}, {m, -l, l}]