9.1: Time-Independent Perturbation Theory

If an atom (not necessarily in its ground state) is placed in an external electric field, the energy levels shift, and the wavefunctions are distorted. This is called the Stark effect. The new energy levels and wavefunctions could in principle be found by writing down a complete Hamiltonian, including the external field, and finding the eigenkets. This actually can be done in one case: the hydrogen atom, but even there, if the external field is small compared with the electric field inside the atom (which is billions of volts per meter) it is easier to compute the changes in the energy levels and wavefunctions with a scheme of successive corrections to the zero-field values. This method, termed perturbation theory, is the single most important method of solving problems in quantum mechanics, and is widely used in atomic physics, condensed matter and particle physics.

It should be noted that there are problems which cannot be solved using perturbation theory, even when the perturbation is very weak, although such problems are the exception rather than the rule. One such case is the one-dimensional problem of free particles perturbed by a localized potential of strength $\lambda$. As we found earlier in the course, switching on an arbitrarily weak attractive potential causes the $k=0$ free particle wavefunction to drop below the continuum of plane wave energies and become a localized bound state with binding energy of order $\lambda^2$. However, changing the sign of $\lambda$ to give a repulsive potential there is no bound state, the lowest energy plane wave state stays at energy zero. Therefore the energy shift on switching on the perturbation cannot be represented as a power series in $\lambda$, the strength of the perturbation. This particular difficulty does not in general occur in three dimensions, where arbitrarily weak potentials do not give bound states—except for certain many-body problems (like the Cooper pair problem) where the exclusion principle reduces the effective dimensionality of the available states.

The Perturbation Series

We begin with a Hamiltonian $H^0$ having known eigenkets and eigenenergies: $$H^0|n^0\rangle=E^0_n|n^0\rangle. \label{9.1.1}$$

The task is to find how these eigenkets and eigenenergies change if a small term $H^1$ (an external field, for example) is added to the Hamiltonian, so: $$(H^0+H^1)|n\rangle=E_n|n\rangle. \label{9.1.2}$$

That is to say, on switching on $H^1$, $$|n^0\rangle\to |n\rangle,\;  E^0_n\to E_n. \label{9.1.3}$$

The basic assumption in perturbation theory is that $H^1$ is sufficiently small that the leading corrections are the same order of magnitude as $H^1$ itself, and the true energies can be better and better approximated by a successive series of corrections, each of order $H^1/H^0$ compared with the previous one.

The strategy, then, is to expand the true wavefunction and corresponding eigenenergy as series in $H^1/H^0$. These series are then fed into $(H^0+H^1)|n\rangle=E_n|n\rangle$, and terms of the same order of magnitude in $H^1/H^0$ on the two sides are set equal. The equations thus generated are solved one by one to give progressively more accurate results.

To make it easier to identify terms of the same order in $H^1/H^0$ on the two sides of the equation, it is convenient to introduce a dimensionless parameter $\lambda$ which always goes with $H^1$, and then expand $|n\rangle,\; E_n$ as power series in $\lambda$, $|n\rangle=|n^0\rangle+\lambda|n^1\rangle+\lambda^2|n^2\rangle+\dots$, etc. The ket $|n^m\rangle$ multiplied by $\lambda^m$ is therefore of order $(H^1/H^0)^m$.

This $\lambda$ is purely a bookkeeping device: we will set it equal to 1 when we are through! It’s just there to keep track of the orders of magnitudes of the various terms.

Putting the series expansions for $|n\rangle,\; E_n$ in $$(H^0+\lambda H^1)|n\rangle=E_n|n\rangle \label{9.1.4}$$

we have

$$\begin{matrix} (H^0+\lambda H^1)(|n^0\rangle+\lambda|n^1\rangle+\lambda^2|n^2\rangle+\dots) \\ =(E^0_n+\lambda E^1_n+\lambda^2E^2_n+\dots)(|n^0\rangle+\lambda|n^1\rangle+\lambda^2|n^2\rangle+\dots). \end{matrix} \label{9.1.5}$$

We’re now ready to match the two sides term by term in powers of $\lambda$.

Zeroth Order of Energy

The zeroth-order term ($\lambda =0$), of course, just gives back

$$H^0|n^0\rangle=E^0_n|n^0\rangle$$

Zeroth Order of Eigenstates

$$|n^0\rangle=|n^0\rangle$$

First-Order of Energy

Matching the terms linear in $\lambda$ on both sides: $$H^0|n^1\rangle+H^1|n^0\rangle=E^0_n|n^1\rangle+E^1_n|n^0\rangle. \label{9.1.6}$$

This equation is the key to finding the first-order change in energy $E^1_n$. Taking the inner product of both sides with $\langle n^0|$ : $$\langle n^0|H^0|n^1\rangle+\langle n^0|H^1|n^0\rangle=\langle n^0|E^0_n|n^1\rangle+\langle n^0|E^1_n|n^0\rangle, \label{9.1.7}$$

then using $\langle n^0|H^0=\langle n^0|E^0_n$, and $\langle n^0|n^0\rangle=1$, we find $$E^1_n=\langle n^0|H^1|n^0\rangle. \label{9.1.8}$$

Taking now $\lambda=1$, we have established that the first-order change in the energy of a state resulting from adding a perturbing term $H^1$ to the Hamiltonian is just the expectation value of $H^1$ in that state.

For example, we can estimate the ground state energy of the helium atom by treating the electrostatic repulsion between the electrons as a perturbation. The zeroth-order ground state has the two (opposite spin) electrons in the ground state hydrogen-atom wavefunction (scaled for the doubling of nuclear charge). The first-order energy correction $E^1_0$ is then given by computing the expectation value $\langle e^2/r_{12}\rangle$ for this ground state wavefunction.

First-Order of Eigenstates

The general expression for the first-order change in the wavefunction is found by taking the inner product of the first-order equation with the bra $\langle m^0|,\;  m\neq n$, $$\langle m^0|H^0|n^1\rangle+\langle m^0|H^1|n^0\rangle=\langle m^0|E^0_n|n^1\rangle+\langle m^0|E^1_n|n^0\rangle. \label{9.1.9}$$

The last term is zero, since $\langle m^0|n^0\rangle=0$, and in the first term $\langle m^0|H^0=\langle m^0|E^0_m$, so $$\langle m^0|n^1\rangle=\frac{\langle m^0|H^1|n^0\rangle}{E^0_n-E^0_m} \label{9.1.10}$$

and therefore the wavefunction correct to first order is:

$$\begin{align} |n\rangle &= |n^0\rangle+|n^1\rangle \\[5pt] &=|n^0\rangle+\sum_{m\neq n} \frac{|m^0\rangle\langle m^0|H^1|n^0\rangle}{E^0_n-E^0_m}. \end{align} \label{9.1.11}$$

Second-Order of Energy

To find the second-order correction to the energy, it is necessary to match the second-order terms in

$$\begin{matrix} (H^0+\lambda H^1)(|n^0\rangle+\lambda|n^1\rangle+\lambda^2|n^2\rangle+\dots) \\ =(E^0_n+\lambda E^1_n+\lambda^2E^2_n+\dots)(|n^0\rangle+\lambda|n^1\rangle+\lambda^2|n^2\rangle+\dots) \end{matrix} \label{9.1.12}$$

giving:

$$H^0|n^2\rangle+H^1|n^1\rangle=E^0_n|n^2\rangle+E^1_n|n^1\rangle+E^2_n|n^0\rangle. \label{9.1.13A}$$

Taking the inner product with $\langle n^0|$ yields:

$$\langle n^0|H^0|n^2\rangle+\langle n^0|H^1|n^1\rangle=E^0_n\langle n^0|n^2\rangle+E^1_n\langle n^0|n^1\rangle+E^2_n\langle n^0|n^0\rangle. \label{9.1.13B}$$

The leading terms on the two sides cancel as before. What about the term $E^1_n\langle n^0|n^1\rangle$? Since $|n\rangle=|n^0\rangle+|n^1\rangle$, and both $|n\rangle$,  $|n^0\rangle$ are normalized, $\langle n^0|n^1\rangle+\langle n^1|n^0\rangle=0$ in leading order—that is to say, $\langle n^0|n^1\rangle$ is pure imaginary. That just means that if to this order $|n\rangle$ has a component parallel to $|n^0\rangle$, that component has a small pure imaginary amplitude, and $|n\rangle$ can be written (to this order) as $|n\rangle=e^{i\alpha} |n^0\rangle+ kets \perp |n^0\rangle$, with $\alpha$ small. But the phase factor can be eliminated by redefining the phase of $|n\rangle$, and with that redefinition $|n^1\rangle$ has no component in the $|n^0\rangle$ direction, we can therefore drop the term $E^1_n\langle n^0|n^1\rangle$.

So the second-order correction to the energy is:

$$\begin{align} E^2_n &=\langle n^0|H^1|n^1\rangle \\[5pt] &= \langle n^0|H^1\sum_{m\neq n} \dfrac{|m^0\rangle\langle m^0|H^1|n^0\rangle}{E^0_n-E^0_m} \\[5pt] &=\sum_{m\neq n} \dfrac{|\langle m^0|H^1|n^0\rangle|^2}{E^0_n-E^0_m}. \end{align} \label{9.1.14}$$

Selection Rules

Perturbation theory involves evaluating matrix elements of operators. Very often, many of the matrix elements in a sum are zero—obvious tests are parity and the Wigner-Eckart theorem. These are examples of selection rules: tests to find if a matrix element may be nonzero.

The Quadratic Stark Effect

When a hydrogen atom in its ground state is placed in an electric field, the electron cloud and the proton are pulled different ways, an electric dipole forms, and the overall energy is lowered. The perturbing Hamiltonian from the electric field is

$$H^1=e​ \mathscr{E}​z=e\mathscr{E}​r\cos\theta$$

where $\mathscr{E}$ is the electric field strength, the field is in the z- direction, the electron charge $e$ is negative.

We shall denote the unperturbed eigenenergies of the hydrogen atom by $E_n=E_{nlm}=-1/n^2$, so in particular we denote the ground state energy by $E_1$.

The first-order correction to the ground state energy $E^1_1=\langle 100|e\mathscr{E}z|100\rangle$, where

$$|100\rangle≡\psi_{100}(r)=\left( \frac{1}{\pi a^3_0}\right)^{1/2}e^{-r/a_0}. \label{9.1.15}$$

This first-order term is zero since there are equal contributions from positive and negative z.

The second-order term is $$E^2_1 =\sum_{n\neq 1;l,m} \frac{|\langle nlm|e\mathscr{E}z|100\rangle|^2}{E_1-E_n} \label{9.1.16}$$

where we are now using $|nlm\rangle$ to denote the unperturbed hydrogen atom wavefunctions, and here the $E_n=-1/n^2$ (in Rydbergs) are the unperturbed energies.

Most of the terms in this infinite series are zero—the selection rules help get rid of them as follows: since $e\mathscr{E}z$ is the $m=0$ component of a spherical vector and $|100\rangle$ is a zero angular momentum state, it follows from the Wigner-Eckart theorem that $\langle nlm|$ can only be $\langle n10|$. This reduces the second-order sum over states to: $$E^2_1 =\sum_{n\neq 1} \frac{|\langle n^10|e\mathscr{E}z|100\rangle|^2}{E_1-E_n}. \label{9.1.17}$$

This is still not easy to evaluate, but an upper bound can be found by observing that $|E_1-E_n|\ge |E_1-E_2|$, so $$\begin{matrix} |E^2_0|<\frac{1}{E_2-E_1}\sum_{n\neq 1} |\langle n10|e\mathscr{E}z|100\rangle|^2 \\ =\frac{1}{E_2-E_1}\sum_{n\neq 1;l,m} \langle 100|e\mathscr{E}z|nlm\rangle\langle nlm|e\mathscr{E}z|100\rangle \end{matrix} \label{9.1.18}$$

where we have temporarily restored the full sum over $n,l,m$ that is, we’ve put back all the zero terms. The reason for this seeming backward step is that, having taken the energy-difference denominator outside the sum, we can even include $|100\rangle$ in the $|nlm\rangle$ sum (it’s another zero term) and in fact we can even include the plane-wave (ionized) states as well as the bound states, since the plane waves all have energy greater than zero. At this point, the $\sum_{n,l,m}$ sum becomes a sum over all states, and therefore just becomes the unit operator, $$\sum_{n,l,m} |nlm\rangle\langle nlm|=I, \label{9.1.19}$$

so $$|E^2_1 |<\frac{1}{E_2-E_1}\langle 100|(e\mathscr{E}z)^2|100\rangle. \label{9.1.20}$$

For the ground state hydrogen wavefunction, $\langle 100|z^2|100\rangle=a^2_0,\; E_1=-e^2/2a_0,\;   E_2=E_1/4$, so $$|E^2_1 |<\frac{1}{(\frac{3}{4}e^2/2a_0)}(e\mathscr{E})^2a^2_0=\frac{8}{3}\mathscr{E}^2a^3_0. \label{9.1.21}$$

Furthermore, since all the terms in the series for $E^2_1$ are negative, the first term sets a lower bound on $|E^2_1 |$:

$$|E^2_1 |> \frac{|\langle 210|e\mathscr{E}z|100\rangle|^2}{E_1-E_2}. \label{9.1.22}$$

This can be evaluated in straightforward fashion to find $|E^2_1 |>0.55\times \frac{8}{3}\mathscr{E}^2a^3_0$.

So, even though we have not actually evaluated the second-order correction to the energy explicitly, we have it bracketed between two values, the lower one being more than half the upper one. Other ingenious methods have been developed (see Shankar or Sakurai) to find that the true answer is $|E^2_1 |=\frac{9}{4}\mathscr{E}^2a^3_0$, but in fact the whole problem can be solved exactly using parabolic coordinates.

Degenerate Perturbation Theory: Distorted 2-D Harmonic Oscillator

The above analysis works fine as long as the successive terms in the perturbation theory form a convergent series. A necessary condition is that the matrix elements of the perturbing Hamiltonian must be smaller than the corresponding energy level differences of the original Hamiltonian. If $H^0$ has different states with the same energy, in other words degenerate energy levels, and the perturbation has nonzero matrix elements between these degenerate levels, then obviously the theory breaks down. To see just how it breaks down, and how to fix it, we consider the two-dimensional simple harmonic oscillator:

$$H^0=\frac{p^2_x+p^2_y}{2m}+\frac{1}{2}m\omega^2(x^2+y^2). \label{9.1.23}$$

Recall that for the one-dimensional simple harmonic oscillator the ground state wavefunction is

$$|0\rangle=\left(\frac{m\omega}{\pi \hbar}\right)^{1/4}e^{-m\omega x^2/2\hbar}=\left(\frac{m\omega}{\pi \hbar}\right)^{1/4}e^{-\xi^2/2} \;\; with\;\; \xi =\sqrt{\frac{m\omega}{\hbar}}x, \;\; and\;\; |1\rangle=\left(\frac{m\omega}{\pi \hbar}\right)^{1/4}\sqrt{2}\xi e^{-\xi^2/2}. \label{9.1.24}$$

The two-dimensional oscillator is simply a product of two one-dimensional oscillators, so, writing $\eta =\sqrt{\frac{m\omega}{\hbar}}y$, the ground state is $|0\rangle=\left(\frac{m\omega}{\pi \hbar}\right)^{1/2}e^{-(\xi^2+\eta^2)/2}$, and the two (degenerate) next states, energy $\hbar\omega$ above the ground state, are

$$|1,0\rangle=\left(\frac{m\omega}{\pi \hbar}\right)^{1/2}\sqrt{2}\xi e^{-(\xi^2+\eta^2)/2}, \;\; |0,1\rangle=\left(\frac{m\omega}{\pi \hbar}\right)^{1/2}\sqrt{2}\eta e^{-(\xi^2+\eta^2)/2}. \label{9.1.25}$$

Suppose now we add a small perturbation

$$H^1=\alpha m\omega^2xy, \label{9.1.26}$$

with $\alpha$ a small parameter.

Notice that

$$\langle 0|H^1|0\rangle=\langle 1,0|H^1|1,0\rangle=\langle 0,1|H^1|0,1\rangle=0$$

so according to naïve perturbation theory, there is no first-order correction to the energies of these states.

However, on going to second-order in the energy correction, the theory breaks down. The matrix element $\langle 1,0|H^1|0,1\rangle$ is nonzero, but the two states $|0,1\rangle,  |1,0\rangle$ have the same energy! This gives an infinite term in the series for $E^2_n$.

Yet we know that a small term of this type will not wreck a two-dimensional simple harmonic oscillator, so what is wrong with our approach? It is helpful to plot the original harmonic oscillator potential $\frac{1}{2}m\omega^2(x^2+y^2)$ together with the perturbing potential $\alpha m\omega^2xy$. The first of course has circular symmetry, the second has axes in the directions $x=\pm y$, climbing most steeply from the origin along $x=y$, falling most rapidly in the directions $x=-y$. If we combine the two potentials into a single quadratic form, the original circles of constant potential become ellipses, with their axes aligned along $x=\pm y$.

The problem arises even in the classical two-dimensional oscillator: picture a ball rolling backwards and forwards in a smooth saucer, a circular bowl. Now imagine the saucer is made slightly elliptical. The ball will still roll backwards and forwards through the center if it is released along one of the axes of the ellipse, although with different periods, as the axes differ in steepness. However, if it is released at a point off the axes, it will describe a complex path resolvable into components in the two axis directions having different periods.

For the quantum oscillator as for the classical one, as soon as the perturbation is introduced, the eigenkets are in the direction of the new elliptic axes. This is a large change from the original x and y axes, and definitely not proportional to the small parameter $\alpha$. But the original unperturbed problem had circular symmetry, and there was no particular reason to choose the x and y axes as we did. If we had instead chosen as our original axes the lines $x=\pm y$, the kets would not have undergone large changes on switching on the perturbation.

The resolution of the problem is now clear: before switching on the perturbation, choose a set of basis kets in a degenerate subspace such that the perturbation is diagonal in that subspace.

In fact, for the simple harmonic oscillator example above, the problem can be solved exactly: $$\begin{matrix} \frac{1}{2}m\omega^2(x^2+y^2)+\alpha m\omega^2xy \\ =\frac{1}{2}m\omega^2 \left[ (1+\alpha)(\frac{x+y}{\sqrt{2}})^2+(1-\alpha)(\frac{x-y}{\sqrt{2}})^2\right] \end{matrix} \label{9.1.27}$$

and it is clear that, despite the results of naïve first-order theory, there is indeed a first order shift in the energy levels, $$\hbar\omega \to \hbar\omega \sqrt{1\pm \alpha}\approx \hbar\omega (1\pm \alpha /2). \label{9.1.28}$$

The Linear Stark Effect

The hydrogen atom, like the two-dimensional harmonic oscillator discussed above, has a nondegenerate ground state but degeneracy in its lowest excited states. Specifically, there are four $n=2$ states, all having energy -1/4 Ryd : $$\begin{matrix} \psi_{200}(r)=\left(\frac{1}{32\pi a^3_0}\right)^{1/2}\left( 2-\frac{r}{a_0}\right)e^{-r/2a_0}, \\ \psi_{210}(r,\theta ,\phi)=\left(\frac{1}{32\pi a^3_0}\right)^{1/2}\left(\frac{r}{a_0}\right)e^{-r/2a_0}cos\theta , \\ \psi_{21\pm 1}(r,\theta ,\phi)=\left(\frac{1}{32\pi a^3_0}\right)^{1/2}\left(\frac{r}{a_0}\right)e^{-r/2a_0}sin\theta e^{\pm i\phi}. \end{matrix} \label{9.1.29}$$

Perturbing this system with an electric field in the z- direction, $H^1=e\mathscr{E}z=e\mathscr{E}​rcos\theta$, note first that naïve perturbation theory predicts no first-order shift in any of these energy levels. However, to second order, there is a nonzero matrix element between two degenerate levels $\langle 200|H^1|210\rangle$. All the other matrix elements between these basis kets in the four-dimensional degenerate subspace are zero, so the only diagonalization necessary is within thetwo-dimensional degenerate subspace spanned by $|200\rangle$,  $|210\rangle$, where $$H^1=\begin{pmatrix} 0&\Delta \\ \Delta&0 \end{pmatrix} \label{9.1.30}$$

with

$$\begin{align} \Delta &=\langle 200|H^1|210\rangle \\[4pt] =e\mathscr{E}\left(\frac{1}{32\pi a^3_0}\right) \int_{0}^{\infty}\left(2-\frac{r}{a_0}\right) \left(\frac{rcos\theta}{a_0}\right)^2 e^{-r/a_0}r^2dr \sin\theta d\theta d\phi \\[4pt] &=-3e\mathscr{E}a_0.\end{align} \label{9.1.31}$$

Diagonalizing $H^1$ within this subspace, then, the new basis states are $(|200\rangle\pm |210\rangle)/\sqrt{2}$ with energy shifts $\pm \Delta$, linear in the perturbing electric field.

The states $|21\pm 1\rangle$ are not changed by the presence of the field to this approximation, so the complete energy map of the $n=2$ states in the electric field has two states at the original energy of -1/4Ryd, one state moved up from that energy by $\Delta$, and one down by $\Delta$.

Notice that the new eigenstates $(|200\rangle\pm |210\rangle)/\sqrt{2}$ are not eigenstates of the parity operator—a sketch of their wavefunctions reveals that in fact they have nonvanishing electric dipole moment $\vec{\mu}$, indeed this is the reason for the energy shift, $\pm \Delta =\mp 3e\mathscr{E}a_0=\mp \vec{\mu}\cdot \vec{\mathscr{E}}$.