Hydrogen
Hydrogen is the most common element in the Universe. Even ourselves, being made of about 70% of water, itself consisting of 66% hydrogen (the terminology hydrogen comes from making water), we are basically hydrogen. It is therefore an important object to understand in some details.
The Bohr model
Spectroscopy in the late 19th century discovered that the radiation from hydrogen consisted of spectral lines, i.e., the emission occurs at precise frequencies.
Empirically, it was found (by Rydberg) that the lines occurred at wavelengths given by the difference of inverse square of integers:
$${\displaystyle {\frac{1}{\lambda}}=R_{\text{H}}\left({\frac{1}{n_{1}^{2}}}-{\frac {1}{n_{2}^{2}}}\right),}$$
with $R_\text{H}$ the Rydberg constant. Quantum mechanics first manifested itself through the efforts of Niels Bohr to find a theoretical basis for this result. The main hurdle to those who attempted such a derivation before him was rooted in the fact that an accelerating charge radiates, a fact well known from classical electrodynamics, which was the most successful theory of the time. Now, a circular orbit is an acceleration, not of speed (modulus) but of velocity (vector). Calculations show that the lifetime of an electron in this context is extremely short and this goes against the observed stability of hydrogen. Bohr was therefore looking for a physical picture that would constrain the motion on orbits of fixed radii, with no possibility to spiral inward the nucleus as a result of dissipating energy. The existence of particular orbits would then explain the discrete set of spectral lines. Bohr found a supporting physical picture in de Broglie's theory of matter that states that the momentum of a particle $p$ can be related to a wave with the corresponding wavelength
$$\lambda={h\over p}$$
so that, writing a condition for the orbits of the electron to corresponding to standing waves
$$2\pi r=n\lambda$$
or, equivalently, $r/\lambda=n/(2\pi)$, leads to postulating the quantization of the angular momentum $L=mvr$ (with $v$ the speed of the electron orbiting at a distance $r$ from the nucleus) into integer multiples of Planck's constant $\hbar$ (we remind that $\hbar\equiv h/(2\pi)$ and $p=mv$):
$${mvr}=n{h\over2\pi}\tag{1}$$
where we used de Broglie's $p={h/\lambda}$. In this equation, $n=1, 2, \ldots$ is an integer called the principal quantum number. Bohr thus settled for the following postulates, having already solved the problem of the stability of hydrogen:
- The electron moves in circular orbits around the atom.
- The orbit is such that the angular momentum of the electron is quantized (in units of Planck's constant $\hbar$).
- The electron can change orbits, by jumping from one to the other, releasing or acquiring the energy difference in the process.
The proof for the validity of these postulates is that they allow to derive the main results of hydrogen's spectroscopy, and, compellingly, to link its empirically derived constant to many of the most important fundamental constants of nature. This is achieved as follows. For an electron of mass $m$ and charge $-e$ (opposite as the proton so their product is $-e^2$), and moving in a circular orbit of radius $r$ in Coulomb's potential $V=-e^2/(4\pi\epsilon_0 r)$, Newton's second law $F=ma$ gives:
$${e^2\over 4\pi\epsilon_0 r^2}={mv^2\over r}\tag{2}$$
since $a=v^2/r$ for a central force, with $v$ the speed of the electron, and $F=-|\nabla V|$.
Equations (1) and (2) make two equations for two unknown $r$ and $v$, which can be easily solved:
$$\begin{align} r&=n^2a_\mathrm{B}\,,\\ v&=\alpha c/n\,, \end{align}$$
where we have defined the Bohr radius (with units of distance):
$$\boxed{\displaystyle a_\mathrm{B}\equiv{4\pi\epsilon_0\hbar^2\over me^2}}$$
and the fine-structure constant (which has no dimension and includes the speed of light $c$ so that $v$ can be related to a fundamental constant of speed):
$$\boxed{\displaystyle\alpha\equiv{1\over4\pi\epsilon_0}{e^2\over \hbar c}\approx{1\over 137}}$$
The measured value of the fine-structure constant is $\alpha^{−1}=137.035999173(35)$ (uncertainty in parentheses) with a great agreement to the theoretical value. The small value of this parameter is the basis for the standard model of elementary particles as it allows a perturbative treatment in terms of a series of processes whose $n$th order is weighted by $\alpha^n$. It also shows, since the speed of the electron is $v=\alpha c$, that relativistic effects can be neglected in a first approximation. If we compute the energy as kinetic $K$ (${}=\frac12mv^2$) plus potential $V$ (Coulomb), we find:
$$E={e^2\over 8\pi\epsilon_0r}-{e^2\over4\pi\epsilon_0 r}=-{e^2\over8\pi\epsilon_0r}$$
(where we used (2) again for the kinetic energy). Substituting the radius $r=n^2a_\mathrm{B}$ in this result yields
$$\boxed{\displaystyle E_n={E_1\over n^2}}$$
where
$$E_1=-{m\over2\hbar^2}\left({e^2\over4\pi\epsilon_0}\right)^2=-{mc^2\alpha^2\over2}\approx-13.6\text{[eV]}$$
This is an important constant that is related to Rydberg's constant and was known already at the time to be the ionization energy of Hydrogen. The 3rd Bohr postulate thus explains the discrete set of spectral lines from a transition between an initial $i$ and final state $f$ that yields an energy for the emitted photon $E_\gamma=E_i-E_f$:
$$E_\gamma=E_1\left({1\over n_f^2}-{1\over n_i^2}\right)$$
and from Planck's formula $E_\gamma=h\nu$ relating energy and frequency, one can finally derive Rydberg's formula for the wavelengths $\lambda=c/\nu$ of the hydrogen spectrum. The underlying mechanism is thus captured in this picture:
So with these simple arguments and manipulations, Bohr was clearly onto something.
Schrödinger's equation
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}$$
so that the radial equation reads:
$${d^2u\over d\rho^2}=\left[1-{\rho_0\over\rho}+{l(l+1)\over\rho^2}\right]u\,.\tag{3}$$
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 (3), 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{4}$$
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. (4), 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{5}$$
Let us consider again the asymptotic behaviour. For $j$ large enough, Eq. (5) 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{6}$$
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. (6) can be excluded. This therefore brings a new constrain:
$$2(j+l+1)-\rho_0=0\,.$$
This introduces another integer, the principal quantum number $n$ (already introduced by Bohr):
$$n\equiv j+l+1$$
for $j$ also an integer, that varies, defining the principal quantum number and the quantity
$$\rho_0=2n\,.$$
In this way, we have linked $n$ to the energy, since $\rho_0$ was itself energy dependent. This yields:
$$E_n=-{\hbar^2\kappa^2\over2m}=-{me^4\over 8\pi^2\epsilon_0^2\hbar^2\rho_0^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 only need to explicit the radial part of the solution, which can be worked out recursively for the various quantum numbers $n=1, 2, \ldots$ with $j=n-l-1$ and
$$R_{nl}(r)={1\over r}\rho^{l+1}e^{-\rho}v(\rho)\,.$$
- 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}]