Hydrogen atom

The hydrogen atom is the atom of the element hydrogen, a quantum physical system made up of a central positive charge carrier (a proton) and a negative charge carrier (an electron) surrounding it. It is the simplest case of a broader class of objects called hydrogenic atoms (or hydrogen-like atoms) which may possess a stronger central charge, that is, a full atomic nucleus.

The Hamiltonian for the Schrödinger equation of the hydrogen atom is

where $M$ is the mass of the nucleus, $m$ is the mass of the electron, $\hbar$ is the reduced Planck constant, $Z$ is the atomic number (the number of protons), $\varepsilon_{0}$ is the vacuum permittivity and $r_{e}$ is the distance of the electron from the origin. The last term is the potential energy due to electromagnetic attraction.

Derivation#

First things first, the end goal is finding the wavefunction of the system. We'll solve it in position representation, $\psi(r)=\braket{ r | \psi }$. This wavefunction is of course dependent on the spatial coordinates, but since the electron carries Spin ($s=1/2$ and $m_{s}=\pm1/2$ specifically), this will also be parameter. If we call $q$ the generalized coordinates of the wavefunction, we can write it as $\psi\equiv \psi(q)=\psi(r_{e},R,s,m_{s})$.

However, the Hamiltonian has no spin term. $s$ and $m_{s}$ never appear. This suggests that the wavefunction that we get by solving $\hat{H}\psi=E\psi$ is actually just the spatial part, with the spin part being a separate piece. Using separation of variables, we can split the two:

We can now solve these separately.

The spin part#

Since the Hamiltonian does not affect or depend on spin in any way, the spin part is just the eigenstates of the generic spin operator $\hat{S}$:

Thankfully, we know the spin: it's $s=1/2$, and this is the special case that is already solved in Spin > Spin 1/2. We know the eigenvectors of both $\hat{S}^{2}$ and $\hat{S}_{z}$ in spin 1/2: they're somewhat up to preference, so we'll choose

And that's the spin part done already. The real challenge is solving the spatial part.

The spatial part#

Since we have a central potential, the system lends itself well to be rewritten in spherical coordinates, as the Schrödinger equation has a general solution in this form. To do so, we make a first Coordinate transformation to the frame of reference of the center of mass. We define the reduced mass as usual:

Then, using $\mu$ and expressing the Laplacian in the relative coordinate $\mathbf{r}=\mathbf{r}_{e}-\mathbf{R}$, the equation becomes

The benefit is that now we only depend on one degree of freedom ($r$) instead of two ($r_{e},R$). (Effectively, it's a single particle located in $r$ of mass $\mu$ subject to the potential). Then, the spherical coordinate transformation from Cartesian coordinates is

With this we can express the Laplacian in spherical coordinates. Writing the potential as $V(r)$ for brevity, the Hamiltonian operator in spherical coordinates is

We now introduce what is arguably the key component of the system: the quantum angular momentum $\hat{L}=\hat{r}\times \hat{p}=-i\hbar (\hat{r}\times \nabla)$¹. It is so important because if we express its Gradient in spherical coordinates, then take the square, we get the equation

We can see that $\hat{L}^{2}$ comes up exactly in the Hamiltonian, so with a simple substitution

With this, the spherical Schrödinger equation for the hydrogen atom is

where the angular dependence of the system is entirely contained inside the angular momentum. As usual, the square angular momentum, the momentum on one axis ($z$) and the Hamiltonian all commute, $[\hat{L}^{2},\hat{L}_{z}]=[\hat{L}^{2},\hat{H}]=[\hat{L}_{z},\hat{H}]=0$, so they all share simultaneous eigenstates. This sets the foundation of our final solution.

For a general solution of this equation in spherical coordinates, using a generic potential $V$, see Equazione di Schrödinger > Soluzione tridimensionale. In there, we set $\psi(r,\theta,\phi)=R(r)Y(\theta,\phi)$, find that the angular component $Y(\theta,\phi)$ of the solution is universal, completely independent on the potential and described by the Spherical harmonics. The radial component instead depends on $V$ and we must solve it uniquely for each the potential. As such, our next goal is to solve the radial component for our $V$. But first, a refresher on the angular part.

Angular part#

The angular components depend on the eigenstates of $\hat{L}^{2}$ and $\hat{L}_{z}$ and are found using the spherical harmonics $Y_{l,m}(\theta,\phi)$:

Using spherical harmonics on $\hat{L}_{z}$ is a bit overkill, since it's really only dependent on $\phi$. In the general derivation, one finds that the harmonics themselves can go through separation of variables as $Y_{l,m}(\theta,\phi)=\Theta_{l,m}(\theta)\Phi_{m}(\phi)$. The $\Phi_{m}(\phi)$ function is pretty easy to solve and is found to be $\Phi_{m}(\phi)=e^{im\phi}$. The $\Theta_{l,m}(\theta)$ function is more complicated is requires invoking the associated Legendre functions $P_{l,m}(\cos \theta)$ to get $\Theta(\theta)=AP_{l,m}(\cos \theta)$ ($A$ is some constant). As such, we can rewrite the $\hat{L}_{z}$ eigenstates as

In the context of the hydrogen atom, the numbers $l$ and $m$ are given special names: $l$ is the azimuthal quantum number and $m$ is the magnetic azimuthal quantum number.

Still, spherical harmonics are generally complex functions. We know very well that only real numbers can be measured from a system, so we need to bridge the gap. We do the usual: take the square modulus

Luckily, the modulus of $\Phi_{m}(\phi)$ is very easy: $\lvert \Phi_{m}(\phi) \rvert^{2}=\lvert e^{im\phi} \rvert^{2}=1/2\pi$. (Notably, it is completely independent of $m$.) Unluckily, the square modulus of $\Theta(\theta)$ isn't. They can, for instance, be found computationally and plotted in graphs or animations. Analytically, we'll remain with

An important detail: since the spherical harmonics for a given $l$ is spherically symmetrical over all $m$, so is the charge distribution (see Spherical harmonics > Symmetry).

Radial part#

The radial part of the equation is is

where we used $\hat{L}^{2}\psi=\hbar ^{2}l(l+1)\psi$ and rewrote the partial derivatives as regular derivatives since $r$ is now the only variable. We set $u(r)\equiv rR(r)$ to simplify things:

Notice how there is a term dependent on the angular momentum number $l$ even in the radial term. What we do is to combine it with the potential term in a single, effective potential $V_\text{eff}$ defined as²

which reduces the equation even further to

or equivalently

This potential looks something like this:

Note that for $l=0$, the potential is purely attractive: nothing is stopping the electron from falling and crashing into the nucleus. But if we add a rotation ($l>0$), suddenly there is a spike of potential beyond a certain point which explains why electrons don't just collide with the nucleus. The value of $r$ at the bottom of each potential well for $l>0$ should, in principle at least, correspond to the "allowed orbits" that the Bohr model identified. As always, the system is bound if $E<0$ and scattering if $E>0$.

Now, solving the radial part is not exactly trivial. First, to solve a (1D) differential equation like the Schrödinger equation we need a boundary condition, which we have not set up until now. We'll set $u(0)=0$, as that prevents the wavefunction from blowing up to infinity when $r=0$³.

For convenience, we introduce a couple of dimensionless variables

This leads to

and so

If we substitute these in the equation (with $V_\text{eff}$ expanded), we get

Noticing $\lambda$ in the middle leads us to

A much simpler form, but still quite abstract. Instead of diving head-first into solving it, we can start by analyzing asymptotic behavior. When $\rho\to \infty$, the equations greatly simplifies to

This a linear second order ODE of the form $\ddot{u}- \frac{1}{4}u=0$ and the solution is an exponential. Thus, for high $\rho$, we have

Using the fact that $u(\rho)$ must be bounded everywhere, we keep only the negative sign. Since we have an asymptotic solution, we can claim that the general solution behaves in the same manner as the asymptote, but scaled by some other function $f(\rho)$. Thus, in general, in and out of the asymptote:

By definition, $f(\rho)$ must contain all of the information about $E$ and $l$, since those are not present in the asymptotic term: $f(\rho)\equiv f_{E,l}(\rho)$. Thus, since

our new equation to solve is

Our next problem is finding $f(\rho)$.

Finding $f(\rho)$#

To start, we again look at asymptotic behavior:

This can only be true if $f(\rho)\sim \rho^{l+1}$. In fact, using the definition of $f(\rho)\sim \rho^{s+1}$:

But only the dominant term $\rho^{s-1}$ matters in the limit so this is like saying

Since $\rho^{s-1}>0$ for all $\rho$, this is true only if $s=l$, which proves that $f(\rho)\sim \rho^{l+1}$.

Outside of the asymptote, $f(\rho)$ must behave something like $f(\rho)=\rho^{l+1}g(\rho)$ where $g(\rho)$ is some function. As an added note, this means that our current understanding of $u(\rho)$ is $u(\rho)=e^{-\rho/2}\rho^{l+1}g(\rho)$. We now substitute this in $(1)$ to get

This equation can be solved analytically by expanding $g(\rho)$ in its power series

Taking the derivatives of this power series yields

Split the sum like so

and the first sum is further split as

Changing the index from $k$ to $k+1$ in the first yields the equivalent form

Plug this in, combine the sums again and collect $\rho^{k}$ to get

This is true for all $\rho$, which means that each component must be zero termwise:

which yields the recurrence relation

For large $k$, the ratio is approximately

As it happens, this ratio is the same as that of the series for $\rho^{p}e^{\rho}$ with $p \in \mathbb{R}$. This implies that $g(\rho)$ has this very behavior

Combined with our previous result $u(\rho)=e^{-\rho/2}\rho^{l+1}g(\rho)$, this implies that the asymptotic behavior of $u(\rho)$ is something of the form $u(\rho)\sim\rho^{l+1+p}e^{\rho/2}$. But this can't be, for it would break our previous assumption of $u(\rho)\sim e^{-\rho/2}$ (and by extension the normalizability of the wavefunction). Thus the series must terminate early, which is to say that it must be a polynomial in $\rho$. We shall call $n_{r}$ the rank of the polynomial, the radial quantum number, and it is a nonnegative integer (possibly zero). This means that any terms of index $k>n_{r}$ must be zero and that the highest term is $c_{n_{r}}\rho^{n_{r}}$.

Using the recurrence relation we can state

Multiply both sides by the denominator to get

but we know that $c_{n_{r}}\neq 0$, so it must be that

We'll rename this to $n$

and we shall call it the principal quantum number. Since $n_{r}=0,1,2,\ldots$ and $l=0,1,2,\ldots$, we have $n=1,2,3,\ldots$.

Eigenvalues#

Now that we finally know what $\lambda$ is, we can invert its definition to extract $E$, substitute $\lambda$ for $n$ and get the energy eigenvalues:

where we used the Bohr radius $a_{0}$ and modified Bohr radius $a_{\mu}=a_{0} (m/\mu)$. Thus, after all this derivation, we have the expression for the allowed energy levels of the hydrogen atom:

Notably, the energy eigenvalues are not dependent on $l$ (not directly, at least). This means that the energy eigenstates are degenerate in $l$ too, not just $m$. For each $n$ there are $n-1$ values of $l$, and for each $l$ there are $2l+1$ values of $m$, which leads to a total eigenvalue degeneracy of

The degeneracy with respect to $m$ is universal for spherically symmetric systems, but the degeneracy in $l$ is actually specific to the Coulomb potential, or rather potentials that go like $V\sim 1/r$. In fact, were the potential to behave differently, the $l$ degeneracy would be lifted and we'd have $n$ distinct states for each principal level $n$. This happens, for instance, in the outer electrons of an atom with more than one electron shell, which are also subject to the repulsion of inner electrons, changing the effective potential to a different shape than $\sim1/r$. You can see this in the Many-electron atom article. In a sense, the $l$ degeneracy has little to do with the shape of an atom specifically: anything around a point charge would work the same.

The levels themselves are infinite and tend towards zero as you get further from the nucleus: this shouldn't be surprising considering that the electromagnetic potential of a point charge (the thing that is responsible for keeping the electron anchored to the nucleus) falls off with distance. They also get progressively closer to one another for larger $n$. The difference between $n=1$ and $n=2$ is bigger than that between $n=2$ and $n=3$ and so on.

Furthermore, these energy levels are identical to those obtained from the Bohr model of the atom. Despite being a (semi)classical model, it already got the quantum levels correct down to the equation. But what it couldn't do, however, are the eigenstates themselves.

Eigenstates#

To compute the eigenstate, we look back at $(2)$. Earlier, we looked only at the asymptotic behavior, which led us to the definition of the radial and principal quantum numbers and $\lambda=n$. We now want to find the actual solutions themselves. This is a second order Ordinary differential equation in $\rho$ and luckily, a similar equation has already been solved in the past.

The (generalized) Laguerre's differential equation is the following second order ODE:

whose solutions $L_{q}^{p}(\rho)$ are the associated Laguerre polynomials

for some nonnegative integers $q$ and $p$. This should ring a bell: it is basically just $(2)$ with slightly different variables. In fact, if we equate the differences we get

Thus, the physically valid solutions of $(2)$ are

where $A$ is some Normalization constant. Substituting $r$ for $\rho$ and substituting $g(\rho)$ into $R(r)$ (after unwinding all the functions we've added)

where in the last step we have absorbed $na_{\mu}/2Z$ inside $A$ since they are both constants. Calculating the constant (see below) yields the radial wavefunction:

Spatial wavefunction#

This completes our proof: we've assembled all the pieces we need to determine the spatial wavefunction of the hydrogenic atom, which is the product of the the radial part and the angular part:

where

Use the definition of $\rho$ to substitute

$$$$

Now, this integral is very hard⁴. Fortunately, someone else did the math already: in the Mathematical Handbook of Formulas and Tables, a special result for the associated Laguerre polynomials (Eq. 30.50) precisely solves this integral. After due substitutions⁵, the result is

\int_{0}^{\infty}\rho^{2l+2}e^{-\rho}[L_{n-l-1}^{2l+1}(\rho)]^{2}d\rho&=\frac{[2(n-l-1)+2l+2](n-l-1+2l+1)!}{(n-l-1)!} \\ &=\frac{2n(n+l)!}{(n-l-1)!} \end{align}

Results#

Each of these wavefunctions represents a possible state of the electron around the nucleus, bringing along the baggage of all the probabilities of where it might be at any given time while "orbiting" the nucleus. In fact, this classical concept of "orbit" remained so prevalent throughout the years that it gave the special name we use to refer to these wavefunctions: orbitals. Common in chemistry too, we use a specific naming scheme to specify which orbital we are talking about, called spectroscopic notation.

As is, the orbitals have limited intepretability. But as with all quantum mechanical system, the square modulus of $\psi_{nlm}$ gives the Probability of finding the system in a given state. Specifically, $\lvert \psi_{nlm}(r,\theta,\phi) \rvert^{2}d\mathbf{r}$ is the probability of finding the electron in an infinitesimal volume element $d\mathbf{r}$ and is thus effectively the (shape of the) volume charge distribution of the hydrogen atom. Meanwhile, the probability of finding it at a distance $r$ anywhere around the nucleus is found by integrating $\lvert \psi_{nlm}(r,\theta,\phi) \rvert^{2}$ over all angles. This leads to what's known as a density function

since $Y_{lm}(\theta,\phi)$ are normalized. The first few radial wavefunction and their corresponding densities look like this.

Other visualizations of the wavefunctions are the $xy$-plane section plot below, where you can see the black lines representing the nodes (i.e. the zeros) of the wavefunctions.

Also see this website, which provides an interactive simulation of the orbitals: https://www.falstad.com/qmatom/.

Total wavefunction#

Now that we painstakingly solved both the spin and spatial parts, we can put them together to officially have the full wavefunction of the hydrogen atom:

where $m$ was renamed to $m_{l}$ to better distinguish it from the spin magnetic number $m_{s}$. It is dependent on five quantum numbers, of which four variable:

• $n$ is the principal quantum number ($n=1,2,3,\ldots$)
• $l$ is the azimuthal quantum number ($l< n$)
• $m_{l}$ is the magnetic azimuthal quantum number ($-l\leq m_{l}\leq l$)
• $s$ is the spin quantum number ($s=1/2$ for electrons)
• $m_{s}$ is the magnetic spin quantum number ($s=\pm 1$ for electrons)

Notably, the spatial part of the wavefunction is a usual complex valued function. However, the spin part is a vector, with two components. This means that the operation we're doing here is scalar multiplication of the spin vector by the spatial wavefunction: in other words, the total wavefunction is a two-dimensional vector, not a number.

For the curious, this is what the wavefunction looks like with all nested functions written out explicitly

The fine structure#

Although the wavefunction above provides good results, it still uses the Schrödinger equation as its starting point. The issue with this is that the Schrödinger equation does not consider the effects of relativity, and since electrons move with average speeds in the order of 1% of the speed of light, this leads to minor, but still significant discrepancies when comparing theory with experiment.

The formally correct way of handling this would be to start from relativistic quantum mechanics and the Dirac equation. However, that is beyond the scope of this section, and it is also historically misplaced. Before the Dirac equation was even developed, relativistic corrections to the hydrogen atom were already discovered and calculated by using perturbation theory, specifically time-independent perturbation theory. Namely, three corrections were developed: relativistic kinetic energy, spin-orbit coupling and the Darwin term. Each of these made the eigenvalues just a little closer to reality, matching the differences that we see in the laboratory. Together, these corrections are known as the fine structure of the hydrogen atom.

This is a very qualitative description of the fine structure. Please refer to e.g. Appendix 7 of Bransden & Joachaim for a rigorous derivation.

Our perturbed Hamiltonian will be

where

is our usual Hamiltonian (the quantity, not the operator, note the lack of hat on $H$) for a charged particle (our electron) in a central Coulomb potential and

will be the perturbation due to the three aforementioned corrections.

Relativistic kinetic energy#

In relativity, the energy of the body is not just kinetic. We must use the full relativistic energy $E=\sqrt{ p^{2}c^{2}+m^{2}c^{4} }$, from which we get the kinetic energy as

If $p^{2}/mc^{2}\ll 1$, which is our case, we can expand the square root in a Taylor series of $p^{2}/mc^{2}$:

Plugging this into the previous equation gives an approximate kinetic energy of

The first term is very clearly the classical kinetic energy, whereas the second is a term that is unique to relativity. This is the relativistic correction to kinetic energy:

The correction to the eigenvalues is given by the average of the perturbation in the generic orbital (i.e. the eigenstates of $H_{0}$)

This equation can be further developed to get

where $\alpha$ is a constant aptly known as the fine-structure constant.

Spin-orbit correction#

The spin-orbit term takes the effect of the interaction between spin and angular momentum of the electron into account. The easiest way to reach it is by considering the frame of reference in which the electron is not moving and the nucleus is rotating around it. In such a frame, the electric and magnetic fields are

The magnetic field can be developed as

The intrinsic magnetic dipole moment of the electron is the Bohr magneton:

Thus, the energy is given by the usual formula for a magnetic dipole:

Evaluating the derivative gives the spin-orbit correction:

This correction does not commute with neither $L_{z}$ nor $S_{z}$, so to find the eigenvalue correction we must first identify some states that diagonalize $H_{2}'$. To do so, we introduce the total angular momentum

which associated quantum numbers $j$ and $m_{j}$. This allows us to write the Scalar product as

A good Basis is the one given by the states

The scalar product in the sum are a well-known set of numbers called the Clebsch-Gordan coefficients for $s=1/2$. In these states we have the eigenvalue equations

Since $s=1/2$, the values for $j$ are $j=s\pm1/2$, while $m_{j}=-j,-j+1,\ldots,j-1,j$ as with all magnetic spin numbers. The energy correction can now be calculated:

It can be proven that

Using

this leads to three possible values:

Darwin term#

The last correction is due to the Heisenberg inequality. Since the concept of "size" is not well-defined in quantum mechanics, the electron is non-localized and is not very well described by a point charge. The Darwin term consists of treating the electron as a cloud of charge described by a volume charge distribution. To do so, we starting by recalling that the usual electric potential for a point charge $e$ goes like

We now introduce a new vector $\mathbf{u}$ which identifies the position of a point in the cloud, starting from the center of the cloud. If $\mathbf{r}$ is the center of the cloud, then $\mathbf{r}+\mathbf{u}$ is a point $\mathbf{u}$ away from the center. The corrected potential energy $\tilde{V}$ will then be a function of $\mathbf{r}+\mathbf{u}$:

where we are integrating over the entire charge. We introduce the volume charge distribution we mentioned before as

This makes our corrected potential energy into

We can expand $\tilde{V}$ in a Taylor series in $\mathbf{r}+\mathbf{u}$ around $\mathbf{r}$ to get