2.7. Hydrogen atom

Here, we solve the time-independent Schrödinger equation for a hydrogenic atom (i.e. an atom with a nucleus of charge \(Ze\) and a single electron).

(2.114)\[\hat{H}\psi(r,\theta,\phi) = E \psi(r,\theta,\phi)\]

If we assume the nucleus fixed (since it is much heavier than the electron), the Hamiltonian is

(2.115)\[\hat{H} = -\frac{\hbar^2}{2m_\mathrm{e}} \nabla^2 + V(r)\]

with the Coulomb potential

(2.116)\[V(r) = -\frac{1}{4\pi\epsilon_0} \frac{Ze^2}{r}\]

Since this is a spherically symmetric potential, the \(\nabla^2\) operator is best expressed in spherical coordinates:

(2.117)\[\nabla^2 = \frac{1}{r} \frac{\partial^2}{\partial r^2} r + \frac{1}{r^2} \left[ \frac{1}{\sin\theta} \frac{\partial}{\partial\theta} \sin\theta \frac{\partial}{\partial\theta} + \frac{1}{\sin^2\theta} \frac{\partial^2}{\partial\phi^2} \right]\]

We abbreviate the angular term in square brackets as \(\hat{\varLambda}^2\) and write

(2.118)\[\nabla^2 = \frac{1}{r} \frac{\partial^2}{\partial r^2} r + \frac{1}{r^2} \hat{\varLambda}^2\]

2.7.1. Separation of variables

The Hamiltonian contains second derivatives with respect to \(r\), \(\theta\), and \(\phi\), but none of them are mixed. Therefore, we can attempt a separation-of-variables approach and use the product ansatz

(2.119)\[\psi(r,\theta,\phi) = R(r)\cdot\varTheta(\theta)\cdot\varPhi(\phi) = R(r)\cdot Y(\theta,\phi)\]

Inserting this into the Schrödinger equation, and pulling \(R\) and \(Y\) in front of the derivatives wherever possible gives

(2.120)\[-\frac{\hbar^2}{2m_\mathrm{e}} \left[ Y \frac{1}{r} \frac{\partial^2}{\partial r^2} \left( rR \right) + \frac{1}{r^2} R \hat{\varLambda}^2 Y \right] + V(r) R Y = E R Y\]

Dividing by \(Y\) gives

(2.121)\[-\frac{\hbar^2}{2m_\mathrm{e}} \left[ \frac{1}{r} \frac{\partial^2}{\partial r^2} \left( rR \right) + \frac{1}{r^2} R \left( \frac{1}{Y} \hat{\varLambda}^2 Y \right) \right] + V(r) = E R\]

The only term that depends on \(\theta\) and \(\phi\) is the one in parentheses. Changing \(\theta\) or \(\phi\) can potentially change the value of this term. However, if this value changed, it would break the equality. Therefore, this term must be constant:

(2.122)\[\frac{1}{Y} \hat{\varLambda}^2 Y = \mathrm{const} \quad \rightarrow \quad \hat{\varLambda}^2 Y = \mathrm{const} \cdot Y\]

We have seen this before - it is the Schrödinger equation of the rigid rotor!

(2.123)\[-\frac{\hbar^2}{2m_\mathrm{e}r^2} \hat{\varLambda}^2 Y = E\cdot Y\]

without the prefactor \(\hbar^2/2m_\mathrm{e}r^2\). Its physical solutions require \(E = (\hbar^2/2m_\mathrm{e}r^2)l(l+1)\) with an integer \(l = 0,1,\dots\). Dividing both sides by \(-\hbar^2/2m_\mathrm{e}r^2\) then gives

(2.124)\[\hat{\varLambda}^2 Y = -l(l+1)Y \quad \rightarrow \quad \frac{1}{Y}\hat{\varLambda}^2 Y = -l(l+1)\]

We next insert this result back into the overall equation. We get a differential equation for the radial function \(R(r)\):

(2.125)\[-\frac{\hbar^2}{2m_\mathrm{e}} \frac{1}{r} \frac{\mathrm{d}}{\mathrm{d}r^2} (rR) + \left[ \frac{\hbar^2}{2m_\mathrm{e}} \frac{l(l+1)}{r^2} - \frac{Ze^2}{4\pi\epsilon_0} \frac{1}{r} \right] R = E\cdot R\]

The first term represent the radial kinetic energy. The term in \([\dots]\) represent the effective potential energy: the first term represents the angular kinetic energy (centrifugal force) and dominates at short distances \(r\) and large \(l\), and the second term is the Coulomb potential energy, which dominates at long \(r\).

Next, we solve this equation for bound states (i.e. \(E <0\)).

2.7.2. The radial equation

First, we simplify the radial equation by rearranging and introducing abbreviations. Multiplication by \(-2m_\mathrm{e}/\hbar^2\) and \(r\) gives

(2.126)\[\frac{\mathrm{d}^2}{\mathrm{d}r^2} (rR) - \left[ \frac{l(l+1)}{r^2} - \frac{2m_\mathrm{e}}{\hbar^2} \frac{Ze^2}{4\pi\epsilon_0} \frac{1}{r} \right] rR = -\frac{2m_\mathrm{e}E}{\hbar^2} r R\]

With the abbreviations \(u = rR\), \(b = (2m_\mathrm{e}/\hbar^2)(Ze^2/4\pi\epsilon_0)\), and \(k=\sqrt{-2m_\mathrm{e}E}/\hbar\) (giving positive \(k\), since \(E\) is always negative), and moving the right-hand-side to the left

(2.127)\[\frac{\mathrm{d}^2 u}{\mathrm{d}r^2} - \left[ k^2+\frac{l(l+1)}{r^2}-\frac{b}{r}\right]u =0\]

Divide by \(4k^2\) and introducing the unitless scaled radius \(\rho = 2kr\) gives

(2.128)\[\frac{\mathrm{d}^2 u}{\mathrm{d}\rho^2} - \left[ \frac{1}{4}+\frac{l(l+1)}{\rho^2}-\frac{b/2k}{\rho}\right]u =0\]

Now we have the equation in a compact form and can start looking for solutions. First, let’s examine the asymptotic behavior. As \(\rho\) goes towards \(\infty\), the \(1/\rho\) and \(1/\rho^2\) terms can be neglected compared to \(1/4\), and we get the simple equation

(2.129)\[\frac{\mathrm{d}^2 u}{\mathrm{d}\rho^2} -\frac{1}{4}u = 0\]

The general solution to this equation is \(u(\rho) = A\mathrm{e}^{-\rho/2} + B\mathrm{e}^{+\rho/2}\). The second term grows exponentially with growing \(\rho\) and is unphysical because it renders the function not square integrable, so we discard it.

As \(\rho\) goes towards zero, the \(1/\rho^2\) term dominates and we can neglect the \(1/\rho\) and \(1/4\) terms, yielding

(2.130)\[\frac{\mathrm{d}^2 u}{\mathrm{d}\rho^2} - \frac{l(l+1)}{\rho^2} u =0\]

The general solution to this equation is \(u(\rho) = C\rho^{l+1}+D\rho^{-l}\), as you can check easily by insertion. The second term of the solution diverges as \(\rho\) approaches zero and renders the function not square integrable. Therefore, it is nonphysical, and we discard it.

With these two asymptotic properties, we can now make the general ansatz

(2.131)\[u(\rho) = L(\rho) \rho^{l+1} \mathrm{e}^{-\rho/2}\]

with a function \(L\) which we still have to determine. Expanding it as a general polynomial

(2.132)\[L(\rho) = \sum_{q=0}^\infty c_q\rho^q = c_0 + c_1\rho + c_2\rho^2 + \cdots\]

and inserting this into the equation gives after some lengthy and tedious algebra the following recursion equation for the coefficients \(c_q\):

(2.133)\[c_{q+1} = c_q \frac{(q+l+1)-b/2k}{(q+1)(q+2l+2)}\]

In this expression, we have two free choices: \(E\) (which enters through \(k\)) and \(c_0\). With those, we can calculate \(c_1\), and then \(c_2\), and so on. We now have a full solution, but still need to check it for physicality. Let’s examine what happens to \(L(\rho)\) at large \(q\). Then we can neglect everything except \(q\) in each of the expressions in parentheses, and also neglect \(b/2k\). In this limit, we get

(2.134)\[c_{q+1} \approx c_q\frac{q}{(q+1)q} = \frac{c_q}{q+1}\]

This implies \(c_q \approx c_0/q!\) and the polynomial asymptotically represents an exponential

(2.135)\[L(\rho)\approx c_0\sum_{q=0}^\infty\frac{1}{q!}\rho^q= c_0 \mathrm{e}^\rho\]

What does this mean for \(u(\rho)\)?

(2.136)\[u(\rho) = L(\rho) \rho^{l+1}\mathrm{e}^{-\rho/2} \approx c_0 \rho^{l+1} \mathrm{e}^{+\rho/2}\]

This diverges with large \(\rho\) because of the exponential. Therefore, the polynomial \(L\) must terminate. The sum cannot be infinite. If we denote the largest \(q\) with non-zero \(c_q\) as \(q_\mathrm{max}\), then we must have \(c_{q_\mathrm{max}+1}=0\). This implies

(2.137)\[(q_\mathrm{max}+l+1) - b/2k = 0\]

Now let’s abbreviate \(n = q_\mathrm{max}+l+1\). This number must be \(l+1\) or larger, since the smallest possible \(q_\mathrm{max}\) is zero. (Note also that this means that \(b/2k\) is an integer.)

Rearranging the equation gives \(k^2 = b^2/(4n^2)\). Now, reinserting the definitions of \(k\) and \(b\) gives

(2.138)\[-\frac{2m_\mathrm{e}}{\hbar^2} = \frac{1}{4} \left( \frac{2m_\mathrm{e}}{\hbar^2} \frac{Ze^2}{4\pi\epsilon_0} \right)^2\]

which finally gives a condition for \(E\)

(2.139)\[E_n = -\frac{1}{4} \left(\frac{2m_\mathrm{e}}{\hbar}\right) \left(\frac{Ze^2}{4\pi\epsilon_0}\right)^2 \cdot \frac{1}{n^2} = -\frac{e^2}{8\pi\epsilon_0 a_0} \cdot \frac{Z^2}{n^2}\]

Taking into account our initial substition \(u = rR\), the physically valid functions are

(2.140)\[R(r) = \frac{1}{r}L(\rho)\rho^{l+1}\mathrm{e}^{-\rho/2}\]

where \(L\) are now not infinite power series, but polynomials of finite degree. These are called the associated Laguerre polynomials. The following table lists them all for \(n\) up to 4.

Symbol	\(n\)	\(l\)	\(L^{2l+1}_{n-l-1}(\rho)\)
1s	1	0	\(L_0^1=1\)
2s	2	0	\(L_1^1=2(2-\rho)\)
2p	2	1	\(L_0^3=6\)
3s	3	0	\(L_2^1=3(6-6\rho+\rho^2)\)
3p	3	1	\(L_1^3=24(4-\rho)\)
3d	3	2	\(L_0^5=120\)
4s	4	0	\(L_3^1=4(24-36\rho+12\rho^2-\rho^3)\)
4p	4	1	\(L_2^3=60(20-10\rho+\rho^2)\)
4d	4	2	\(L_1^5=720(6-\rho)\)
4f	4	3	\(L_0^7=5040\)

If we add the normalization requirement

(2.141)\[\int_0^\infty R_{n,l}^2(r)r^2\mathrm{d}r = 1\]

then the functions are

(2.142)\[R_{n,l}(r) = N_{n,l}L_{n-l-1}^{2l+1}(\rho)\rho^l\mathrm{e}^{-\rho/2}\]

with the normalization constant given by

(2.143)\[N_{n,l} = \sqrt{ \left( \frac{2Z}{na_0} \right)^3 \frac{(n-l-1)!}{2n[(n+l)!]^3} }\]