# The Hydrogen Atom

*Michael Fowler, UVa*

### Factoring Out the Center of Mass Motion

The hydrogen atom consists of two particles, the proton and the electron, interacting via the Coulomb potential $V\left({\overrightarrow{r}}_{1}-{\overrightarrow{r}}_{2}\right)={e}^{2}/r$, where as usual $r=\left|{\overrightarrow{r}}_{1}-{\overrightarrow{r}}_{2}\right|$. Writing the masses of the two particles as ${m}_{1},\text{\hspace{0.17em}}{m}_{2}$ Schrödinger’s equation for the atom is:

$$\left(-\frac{{\hslash}^{2}}{2{m}_{1}}{\overrightarrow{\nabla}}_{1}^{2}-\frac{{\hslash}^{2}}{2{m}_{2}}{\overrightarrow{\nabla}}_{2}^{2}-\frac{{e}^{2}}{r}\right)\psi \left({\overrightarrow{r}}_{1},{\overrightarrow{r}}_{2}\right)=E\psi \left({\overrightarrow{r}}_{1},{\overrightarrow{r}}_{2}\right).$$

But $\text{}{\overrightarrow{r}}_{1},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\overrightarrow{r}}_{2}$ are *not *the
most natural position variables for describing this system: since the potential
depends only on the relative position, a better choice is $\overrightarrow{r},\text{\hspace{0.17em}}\overrightarrow{R}$ defined by:

$$\overrightarrow{r}={\overrightarrow{r}}_{1}-{\overrightarrow{r}}_{2},\text{\hspace{1em}}\overrightarrow{R}=\frac{{m}_{1}{\overrightarrow{r}}_{1}+{m}_{2}{\overrightarrow{r}}_{2}}{{m}_{1}+{m}_{2}}$$

so $\overrightarrow{R}$ is the center of mass of the system. It is convenient at the same time to denote the total mass by $M={m}_{1}+{m}_{2},$ and the reduced mass by $m=\frac{{m}_{1}{m}_{2}}{{m}_{1}+{m}_{2}}.$

Transforming in straightforward fashion to the variables $\overrightarrow{r},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\overrightarrow{R}$ Schrödinger’s equation becomes

$$\left(-\frac{{\hslash}^{2}}{2M}{\overrightarrow{\nabla}}_{R}^{2}-\frac{{\hslash}^{2}}{2m}{\overrightarrow{\nabla}}_{r}^{2}-\frac{{e}^{2}}{r}\right)\psi \left(\overrightarrow{R},\text{\hspace{0.17em}}\overrightarrow{r}\right)=E\psi \left(\overrightarrow{R},\text{\hspace{0.17em}}\overrightarrow{r}\right).$$

Writing the wave function

$\psi \left(\overrightarrow{R},\text{\hspace{0.17em}}\overrightarrow{r}\right)=\Psi \left(\overrightarrow{R}\right)\psi \left(\overrightarrow{r}\right)$

we can split the equation into two:

$\begin{array}{l}\left(-\frac{{\hslash}^{2}}{2M}{\overrightarrow{\nabla}}_{R}^{2}\right)\Psi \left(\overrightarrow{R}\right)={E}_{R}\Psi \left(\overrightarrow{R}\right)\\ \\ \left(-\frac{{\hslash}^{2}}{2m}{\overrightarrow{\nabla}}_{r}^{2}+V\left(\overrightarrow{r}\right)\right)\psi \left(\overrightarrow{r}\right)={E}_{r}\psi \left(\overrightarrow{r}\right)\end{array}$

and the total system energy is $E={E}_{R}+{E}_{r}.$ Note that the motion of the center of mass is
(of course) just that of a free particle, having a trivial plane wave
solution. From now on, we shall only be
concerned with the *relative* motion of
the particles. Since the proton is far
heavier than the electron, we will almost always ignore the difference between
the electron mass and the reduced mass, but it should be noted that the
difference is easily detectable spectroscopically: for example, the lines shift
if the proton is replaced by a deuteron (heavy hydrogen).

We’re ready to write Schrödinger’s equation for the hydrogen atom, dropping the $r$ suffixes in the second equation above, and writing out ${\overrightarrow{\nabla}}^{2}$ explicitly in spherical coordinates:

$$\begin{array}{l}-\frac{{\hslash}^{2}}{2m}\left(\frac{1}{{r}^{2}}\frac{\partial}{\partial r}\left({r}^{2}\frac{\partial \psi}{\partial r}\right)+\frac{1}{{r}^{2}\mathrm{sin}\theta}\frac{\partial}{\partial \theta}\left(\mathrm{sin}\theta \frac{\partial \psi}{\partial \theta}\right)+\frac{1}{{r}^{2}{\mathrm{sin}}^{2}\theta}\frac{{\partial}^{2}\psi}{\partial {\phi}^{2}}\right)-\frac{{e}^{2}}{r}\psi \\ =E\psi .\end{array}$$

### Factoring Out the Angular Dependence: the Radial Equation

Since the potential is spherically symmetric, the Hamiltonian $H$ commutes with the angular momentum operators ${L}^{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{L}_{z}$ so we can construct a common set of eigenkets of the three operators $H,\text{\hspace{0.17em}}\text{\hspace{0.17em}}{L}^{2},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{L}_{z}.$ The angular dependence of these eigenkets is therefore that of the ${Y}_{l}^{m}$ ’s, so the solutions must be of the form

${\psi}_{Elm}\left(r,\theta ,\varphi \right)={R}_{Elm}\left(r\right){Y}_{l}^{m}\left(\theta ,\varphi \right).$

Now, notice that in the Schrödinger equation above, the angular
part of ${\overrightarrow{\nabla}}^{2}$ is *exactly*
the differential operator ${L}^{2}/2m{r}^{2}$,
so operating on ${\psi}_{Elm}\left(r,\theta ,\varphi \right)={R}_{Elm}\left(r\right){Y}_{l}^{m}\left(\theta ,\varphi \right)$ it will give ${\hslash}^{2}l\left(l+1\right)/2m{r}^{2}$.

The spherical harmonic ${Y}_{l}^{m}$ can then be cancelled from the two sides of the equation leaving:

$$\begin{array}{l}-\frac{{\hslash}^{2}}{2m}\left(\frac{1}{{r}^{2}}\frac{d}{dr}\left({r}^{2}\frac{d}{dr}\right)-\frac{l\left(l+1\right)}{{r}^{2}}\right){R}_{El}\left(r\right)-\frac{{e}^{2}}{r}{R}_{El}\left(r\right)\\ =E{R}_{El}\left(r\right)\end{array}$$

it now being apparent that $R\left(r\right)$ cannot depend on $m.$

The radial derivatives simplify if one factors out $1/r$ from the function $R,$ writing

$${R}_{El}\left(r\right)=\frac{u\left(r\right)}{r}$$

and temporarily suppressing the $E$ and $l$ to reduce clutter.

The equation becomes:

$$-\frac{{\hslash}^{2}}{2m}\left(\frac{{d}^{2}}{d{r}^{2}}-\frac{l\left(l+1\right)}{{r}^{2}}\right)u\left(r\right)-\frac{{e}^{2}}{r}u\left(r\right)=Eu\left(r\right).$$

Rearranging,

$$-\frac{{\hslash}^{2}}{2m}\frac{{d}^{2}u\left(r\right)}{d{r}^{2}}+\left(\frac{{\hslash}^{2}}{2m}\frac{l\left(l+1\right)}{{r}^{2}}-\frac{{e}^{2}}{r}\right)u\left(r\right)=Eu\left(r\right).$$

Note that this is the same as the Schrödinger equation for a particle in one dimension, restricted to $r0$, in a potential (for $l\ne 0$ ) going to positive infinity at the origin, then negative and going to zero at large distances, so it always has a minimum for some positive $r.$

We are interested in *bound
states* of the proton-electron system, so $E$ will be a negative quantity. At large separations, the wave equation
simplifies to

$$-\frac{{\hslash}^{2}}{2m}\frac{{d}^{2}u\left(r\right)}{d{r}^{2}}\cong E\text{\hspace{0.05em}}u\left(r\right)\text{(forlarge}r\text{)}$$

having approximate solutions ${e}^{\kappa r},\text{\hspace{1em}}{e}^{-\kappa r}\text{where}\kappa =\sqrt{-2mE/{\hslash}^{2}}.$

The
bound states we are looking for, of course, have exponentially *decreasing* wave functions at large
distances.

### Going to a Dimensionless Variable

To further simplify the equation, we introduce the dimensionless variable

$\rho =\kappa r,\text{\hspace{1em}}\kappa =\sqrt{-2mE/{\hslash}^{2}}$

giving

$$\frac{{d}^{2}u\left(\rho \right)}{d{\rho}^{2}}=\left(1-\frac{2\nu}{\rho}+\frac{l\left(l+1\right)}{{\rho}^{2}}\right)u\left(\rho \right)$$

where (for reasons which will become apparent shortly) we have introduced $\nu $ defined by

$2\nu ={e}^{2}\kappa /E.$

Notice that in transforming from $r$ to the dimensionless variable $\rho $ *the
scaling factor $\kappa $ depends on energy*, so will be different
for different energy bound states!

Consider now the behavior of the wave function near the origin. The dominant term for sufficiently small $\rho $ is the centrifugal one, so

$$\frac{{d}^{2}u\left(\rho \right)}{d{\rho}^{2}}\cong \frac{l\left(l+1\right)}{{\rho}^{2}}u\left(\rho \right)$$

for which the solutions are $u\left(\rho \right)\sim {\rho}^{-l},\text{\hspace{1em}}u\left(\rho \right)\sim {\rho}^{l+1}.$ Since the wave function cannot be singular, we choose the second.

We have established that the wave function decays as ${e}^{-\kappa r}={e}^{-\rho}$ at large distances, and goes as ${\rho}^{l+1}$ close to the origin. Factoring out these two asymptotic behaviors, define $w\left(\rho \right)$ by

$u\left(\rho \right)={e}^{-\rho}{\rho}^{l+1}w\left(\rho \right).$

It is straightforward (if tedious!) to establish that $w\left(\rho \right)$ satisfies the differential equation:

$$\rho \frac{{d}^{2}w\left(\rho \right)}{d{\rho}^{2}}+2\left(l+1-\rho \right)\frac{dw\left(\rho \right)}{d\rho}+2\left(\nu -\left(l+1\right)\right)w\left(\rho \right)=0.$$

Putting in a trial series solution $w\left(\rho \right)={\displaystyle \sum _{k=0}^{\infty}{w}_{k}{\rho}^{k}}$ gives a recurrence relation between successive coefficients:

$$\frac{{w}_{k+1}}{{w}_{k}}=\frac{2\left(k+l+1-\nu \right)}{\left(k+1\right)\left(k+2\left(l+1\right)\right)}.$$

For large values of $k,$ ${w}_{k+1}/{w}_{k}\to 2/k$, so ${w}_{k}\sim {2}^{k}/k!$ and therefore $w\left(\rho \right)\sim {e}^{2\rho}$. This means we have found the diverging radial wavefunction $u\left(\rho \right)\sim {e}^{\rho}$, which is in fact the correct behavior for general values of the energy.

To find the bound states, we must choose energies such that
the series is *not* an infinite
one. As long as the series stops
somewhere, the exponential decrease will eventually take over, and yield a
finite (bound state) wave function. Just
as for the simple harmonic oscillator, this can only happen if for some $k,\text{\hspace{0.17em}}{w}_{k+1}=0.$ Inspecting the ratio ${w}_{k+1}/{w}_{k}$,
evidently the *condition for a bound state
is* that

$\nu =n,\text{aninteger}$

in which case the
series for $w\left(\rho \right)$ terminates at $k=n-l-1.$ *From now on, *since we know that for the
functions we’re interested in $\nu $ is an integer, *we replace $\nu $ by $n.$ *

To find the *energies*
of these bound states, recall $2n=2\nu ={e}^{2}\kappa /E$ and $\kappa =\sqrt{-2mE/{\hslash}^{2}}$, so

$$4{n}^{2}=\frac{{e}^{4}{\kappa}_{n}{}^{2}}{{E}_{n}{}^{2}}=-\frac{{e}^{4}}{{E}_{n}{}^{2}}\frac{2m{E}_{n}}{{\hslash}^{2}},\text{}$$

so

$$\text{}{E}_{n}=-\frac{m{e}^{4}}{2{\hslash}^{2}}\frac{1}{{n}^{2}}=-\frac{13.6}{{n}^{2}}\text{ev=}-\frac{1}{{n}^{2}}\text{Rydberg}.$$

(This defines the *Rydberg*,
a popular unit of energy in atomic physics.)

Remarkably, this is the very same series of bound state
energies found by Bohr from his model!
Of course, this had better be the case, since the series of energies
Bohr found correctly accounted for the spectral lines emitted by hot hydrogen
atoms. Notice, though, that there are
some important differences with the Bohr model: the energy here is determined
entirely by $n,$ called the *principal
quantum number*, but, in contrast to Bohr’s model, $n$ is *not*
the angular momentum. The *true* ground state of the hydrogen atom, $n=1,$ has *zero*
angular momentum: since $n=k+l+1$, $n=1$ means both $l=0$ and $k=0.$ The ground state wave function is therefore
spherically symmetric, and the function $w\left(\rho \right)={w}_{0}$ is just a constant. Hence $u\left(\rho \right)=\rho {e}^{-\rho}{w}_{0}$ and the actual radial wave function is this
divided by $r,$ and of course suitably normalized.

To write the wave function in terms of $r,$ we need to find $\kappa $ (the energy-dependent scaling factor we used in going to a dimensionless variable). Putting together $\rho ={\kappa}_{n}r$, ${\kappa}_{n}=\sqrt{-2m{E}_{n}/{\hslash}^{2}}$ and ${E}_{n}=-\frac{m{e}^{4}}{2{\hslash}^{2}}\frac{1}{{n}^{2}}$ ,

${\kappa}_{n}=\sqrt{2m\frac{m{e}^{4}}{2{\hslash}^{2}}\frac{1}{{n}^{2}}}/\hslash =\frac{m{e}^{2}}{{\hslash}^{2}n}=\frac{1}{{a}_{0}n},$

where the length

${a}_{0}=\frac{{\hslash}^{2}}{m{e}^{2}}=0.529\times {10}^{-10}m.$

is called the *Bohr
radius*: it is in fact the radius of
the lowest orbit in Bohr’s model.

*Exercise*: check this last statement.

It is worth noting at this point that the energy levels can
be written in terms of the Bohr radius *a*_{0}:

$${E}_{n}=-\frac{{e}^{2}}{2{a}_{0}}\frac{1}{{n}^{2}}.$$

(This is actually obvious: remember that the energies ${E}_{n}$ are identical to those in the Bohr model, in which the radius of the ${n}^{\text{th}}$ orbit is ${n}^{2}{a}_{0},$ so the electrostatic potential energy is $-{e}^{2}/n{a}_{0},$ etc.)

Moving on to the excited states: for $n=2,$ we have a choice: *either* the radial function $w\left(\rho \right)$ can have one term, as before, but now the
angular momentum $l=1$ (since $n=k+l+1$ ); *or* $w\left(\rho \right)$ can have *two*
terms (so $k=1$ ), and $l=0.$ Both options give the same energy, -0.25 Ry, since $n$ is the same, and the energy only depends on $n.$

In
fact, there are *four* states at this
energy, since $l=1$ has states with $m=1,\text{\hspace{0.17em}}m=0$ and $m=-1,$ and $l=0$ has the one state $m=0.$ (For the moment, we are not counting the extra
factor of 2 from the two possible *spin*
orientations of the electron.)

For $n=3,$ there are 9 states altogether: $$ gives one, $l=1$ gives 3 and $l=2$ gives 5 different $m$ values. In fact, for principal quantum number $n$ there are ${n}^{2}$ degenerate states. ( ${n}^{2}$ being the sum of the first $n$ odd integers.)

The states can be mapped out, energy vertically, angular momentum horizontally:

The energy $E=-1/{n}^{2},$ the levels are labeled $nl,\text{\hspace{0.17em}}n$ being the principal quantum number and the traditional notation for angular momentum $l$ is given at the bottom of the diagram. The two red vertical arrows are the first two transitions in the spectroscopic Balmer series, four lines of which gave Bohr the clue that led to his model. The corresponding series of transitions to the $1s$ ground state are in the ultraviolet, they are called the Lyman series.

### Wave Functions for some Low-*n*
States

From now on, we label the wave functions with the quantum numbers, ${\psi}_{nlm}\left(r,\theta ,\varphi \right)$, so the ground state is the spherically symmetric ${\psi}_{100}\left(r\right)$.

For this state $R\left(r\right)=u\left(r\right)/r$, where $u\left(\rho \right)={e}^{-\rho}{\rho}^{l+1}w\left(\rho \right)={e}^{-\rho}\rho {w}_{0}$, with ${w}_{0}$ a constant, and $\rho ={\kappa}_{1}r=r/{a}_{0}$.

So, as a function of $r,$ ${\psi}_{100}\left(r\right)=N{e}^{-r/{a}_{0}}$ with $N$ an easily evaluated normalization constant:

${\psi}_{100}\left(r\right)={\left(\frac{1}{\pi {a}_{0}^{3}}\right)}^{1/2}{e}^{-r/{a}_{0}}.$

For $n=2,\text{\hspace{0.17em}}l=1$ the function $w\left(\rho \right)$ is still a single term, a constant, but now $u\left(\rho \right)={e}^{-\rho}{\rho}^{l+1}w\left(\rho \right)={e}^{-\rho}{\rho}^{2}{w}_{0}$, and, for $n=2,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\rho =\kappa r=r/2{a}_{0}$, remembering the energy-dependence of $\kappa .$

Therefore ${\psi}_{210}\left(r,\theta ,\varphi \right)=N\left(\frac{r}{{a}_{0}}\right){e}^{-r/2{a}_{0}}\mathrm{cos}\theta $. Again, evaluating the normalization constant is routine, yielding

${\psi}_{210}\left(r,\theta ,\varphi \right)={\left(\frac{1}{32\pi {a}_{0}^{3}}\right)}^{1/2}\left(\frac{r}{{a}_{0}}\right){e}^{-r/2{a}_{0}}\mathrm{cos}\theta $.

The wave functions for the other $m$ -values, ${\psi}_{21\pm 1}\left(r,\theta ,\varphi \right),$ have the $\mathrm{cos}\theta $ in ${\psi}_{210}$ replaced by $\mp \left(1/\sqrt{2}\right)\mathrm{sin}\theta {e}^{\pm i\varphi}$ respectively (from the earlier discussion of the ${Y}_{l}^{m}$ ’s).

The other $n=2$ state has *$l=0,$* so from $n=k+l+1,$ we have $k=1$ and the series for $w$ has two terms, $k=0$ and $k=1,$ the ratio being

$$\frac{{w}_{k+1}}{{w}_{k}}=\frac{2\left(k+l+1-n\right)}{\left(k+1\right)\left(k+2\left(l+1\right)\right)}=-1$$

for the relevant values: $k=0,\text{\hspace{0.17em}}l=0,\text{\hspace{0.17em}}n=2.$ So ${w}_{1}=-{w}_{0},\text{\hspace{1em}}w\left(\rho \right)={w}_{0}\left(1-\rho \right)$. For $n=2,$ $\rho =r/2{a}_{0}$, the normalized wave function is

$${\psi}_{200}\left(r\right)={\left(\frac{1}{32\pi {a}_{0}^{3}}\right)}^{1/2}\left(2-\frac{r}{{a}_{0}}\right){e}^{-r/2{a}_{0}}.$$

Note that the zero angular momentum wave functions are nonzero and have nonzero slope at the origin. This means that the full three dimensional wave functions have a slope discontinuity there! But this is fine$\u2014$the potential is infinite at the origin. (Actually, the proton is not a point charge, so really the kink will be smoothed out over a volume of the size of the proton$\u2014$a very tiny effect.)

### General Solution of the Radial Equation

In practice, the first few radial functions $w\left(\rho \right)$ can be constructed fairly easily using the method presented above, but it should be noted that the differential equation for $w\left(\rho \right)$

$\rho \frac{{d}^{2}w\left(\rho \right)}{d{\rho}^{2}}+2\left(l+1-\rho \right)\frac{dw\left(\rho \right)}{d\rho}+2\left(n-\left(l+1\right)\right)w\left(\rho \right)=0$

is in fact Laplace’s equation, usually written

$\left(z\frac{{d}^{2}}{d{z}^{2}}+\left(k-1-z\right)\frac{d}{dz}+p\right){L}_{p}^{k}\left(z\right)=0$

where $k,p$ are integers, and ${L}_{p}^{k}\left(z\right)$ is a *Laguerre
polynomial * (Messiah, page 482).

The two equations are the same if *z* = 2*ρ*,
and the solution to the radial equation is therefore

${w}_{nl}\left(\rho \right)={L}_{n-l-1}^{2l+1}\left(2\rho \right).$

Quoting Messiah, the Laguerre polynomials ${L}_{p}^{0}\left(z\right)$,
and the *associated Laguerre polynomials*
${L}_{p}^{k}\left(z\right)$ are given by:

$\begin{array}{l}{L}_{p}^{0}\left(z\right)={e}^{z}\frac{{d}^{p}}{d{z}^{p}}{e}^{-z}{z}^{p}\\ {L}_{p}^{k}\left(z\right)={\left(-1\right)}^{k}\frac{{d}^{k}}{d{z}^{k}}{L}_{p+k}^{0}\left(z\right).\end{array}$

(These representations can be found neatly by solving
Laplace’s equation using $\u2013$ surprise $\u2013$ a

$\underset{0}{\overset{\infty}{\int}}{e}^{-z}{z}^{k}{L}_{p}^{k}}{L}_{q}^{k}dz=\frac{{\left[\left(p+k\right)!\right]}^{3}}{p!}{\delta}_{pq}.$

But what do they look like? The function ${e}^{-z}{z}^{p}$ is zero at the origin (apart from the trivial case $p=0$ ) and zero at infinity, always positive and having nonzero slope except at its maximum value, $z=p$. The $p$ derivatives bring in $p$ separated zeroes, easily checked by sketching the curves generated by successive differentiation. Therefore, ${L}_{p}^{0}\left(z\right)$, a polynomial of degree $p,$ has $p$ real positive zeroes, and value at the origin ${L}_{p}^{0}\left(0\right)=p!$, since the only nonzero term at $z=0$ is that generated by all $p$ differential operators acting on ${z}^{p}$.

The associated Laguerre polynomial ${L}_{p}^{k}\left(z\right)$ is generated by differentiating ${L}_{p+k}^{0}\left(z\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}k$ times. Now
${L}_{p+k}^{0}\left(z\right)$ has $p+k$ real positive zeroes, differentiating it gives
a polynomial one degree lower, with zeroes which must be one in each interval
between the zeroes of ${L}_{p+k}^{0}\left(z\right)$. This argument remains valid for successive
differentiations, so ${L}_{p}^{k}\left(z\right)$ *must
have p real separate zeroes*.

Putting all this together, and translating back from $\rho $ to $r,$ the radial solutions are:

${R}_{nl}\left(r\right)=N{e}^{-r/n{a}_{0}}{\left(\frac{r}{n{a}_{0}}\right)}^{l}{L}_{n-l-1}^{2l+1}\left(\frac{2r}{n{a}_{0}}\right)$

with $N$ the normalization constant.

Here are the three $n=3$ radial wave functions:

The number of nodes, the radial quantum number, is $3-l-1.$ (*Note*:
The relative normalizations are correct here, but not the overall
normalization.)

For higher $n$ values, the wave functions become reminiscent of classical mechanics. for example, for $n=10,$ the highest angular momentum state probability distribution peaks at $r=100{a}_{0}$, the Bohr orbit radius:

whereas for $n=10,\text{\hspace{0.17em}}l=0,$ we find:

Notice this peaks just below twice the Bohr radius. This can be understood from classical mechanics: for an inverse square force law, elliptical orbits with the same semimajor axis have the same energy. The $l=n-1$ orbit is a circle, the $l=0$ orbit is a long thin ellipse (one end close to the proton), so it extends almost twice as far from the origin as the circle. Furthermore, the orbiting electron will spend longer at the far distance, since it will be moving very slowly.

(*Note*: the
normalizations in the above graphs are only approximate.)