20 Time-independent perturbation theory
There are a very limited number of quantum-mechanical systems which we can solve analytically from start to finish; to deal with realistic systems, we usually have to resort to some sort of approximation. We’ve seen one such method already, the WKB approximation. However useful WKB is, it can only be applied under fairly special conditions. The general idea of perturbation theory, in which we systematically expand our solutions in some available small parameter, is much more general.
The basic idea of perturbation theory should be familiar to you already - at its core, it just boils down to series expansion in a small parameter. But there are some important subtleties specific to quantum mechanics, especially once time dependence enters in. We’ll start with “time-independent” perturbation theory.
20.1 Rayleigh-Schrödinger perturbation theory
The setup here is very simple: we assume that we are interested in a system described by a Hamiltonian \hat{H} = \hat{H}_0 + \lambda \hat{W}, where \lambda is some small, continuous parameter. We assume that in the \lambda \rightarrow 0 limit, we have a full analytic solution to the system described by \hat{H}_0, i.e. we have solved the eigenvalue problem \hat{H}_0 \ket{n^{(0)}} = E_n^{(0)} \ket{n^{(0)}}. Our new task is to find the eigenvalues and eigenstates of the full Hamiltonian, (\hat{H}_0 + \lambda \hat{W}) \ket{n} = E_n \ket{n}.
We haven’t said anything about what \lambda is, exactly; the nature of \lambda depends strongly on what problem we’re trying to solve. Most commonly, \lambda is a natural expansion parameter, some small dimensionless number which appears in the full Hamiltonian; for example, the fine-structure constant \alpha = e^2 / (\hbar c) \approx 1/137 is an excellent expansion parameter if we’re studying corrections due to electromagnetic effects.
Sometimes, \lambda isn’t even anything physical at all: it can be used as an auxiliary variable, where we take an existing \hat{W}, multiply by arbitrary \lambda, do the expansion and set \lambda = 1 at the end. On the other extreme, \lambda can even be a real tunable parameter, like the strength of an electric field in your experimental apparatus.
To proceed, we search for a power-series solution in \lambda for both the eigenenergies and eigenstates of the full Hamiltonian: \ket{n} = \ket{n^{(0)}} + \lambda \ket{n^{(1)}} + \lambda^2 \ket{n^{(2)}} + ... \\ E_n = E_n^{(0)} + \lambda E_n^{(1)} + \lambda^2 E_n^{(2)} + ... We can take these expansions and simply plug back into the (time-independent) Schrödinger equation: (\hat{H}_0 + \lambda \hat{W}) \left( \ket{n^{(0)}} + \lambda \ket{n^{(1)}} + \lambda^2 \ket{n^{(2)}} + ...\right) \\ = \left( E_n^{(0)} + \lambda E_n^{(1)} + \lambda^2 E_n^{(2)} + ... \right) \left( \ket{n^{(0)}} + \lambda \ket{n^{(1)}} + \lambda^2 \ket{n^{(2)}} + ... \right) Now we simply proceed order by order. At order \lambda^0, we recover the \hat{H}_0 eigenvalue equation in full. First order is more interesting: we find (\hat{H}_0 - E_n^{(0)}) \ket{n^{(1)}} = (E_n^{(1)} - \hat{W}) \ket{n^{(0)}} and continuing to second order, (\hat{H}_0 - E_n^{(0)}) \ket{n^{(2)}} = (E_n^{(1)} - \hat{W}) \ket{n^{(1)}} + E_n^{(2)} \ket{n^{(0)}}. Clearly we can keep going as far as we want, although things start to get messy beyond second order; in cases where you need to go to high order, Brillouin-Wigner perturbation theory (to be described below) is a better choice to set things up.
20.1.1 First-order energy eigenvalues and eigenstates
For now, we assume that all of the eigenvalues E_n^{(0)} of \hat{H}_0 are non-degenerate; we’ll come back to deal with that complication later. We also assume that we’re dealing with a discrete system (although if n becomes a continuous label instead, all the derivations are very similar, replacing sums with integrals as appropriate.)
Now, notice that if we multiply by the zeroth-order eigenstate on the left side, we have \bra{n^{(0)}} (\hat{H}_0 - E_n^{(0)}) \ket{n^{(1)}} = 0, letting \hat{H}_0 act to the left. Therefore, we see that \bra{n^{(0)}} (E_n^{(1)} - \hat{W}) \ket{n^{(0)}} = 0, i.e. the zeroth-order energy has been shifted by the amount \Delta_n^{(1)} = \lambda E_n^{(1)} = \bra{n^{(0)}} \lambda \hat{W} \ket{n^{(0)}}. This proves the result we borrowed before, namely that the first-order energy shift in perturbation theory is simply the expectation value of the perturbation with respect to the old eigenstates. Intuitively, it makes sense that we don’t need the first-order states \ket{n^{(1)}} for this: the effect of \hat{W} on the eigenstates themselves begins at order \lambda, and so corrections in expectation values due to the states shifting will only appear starting at second order.
Now let’s find the eigenstates. Since we presumably know what all of the zeroth-order eigenstates are, a reasonable approach would be to try and expand \ket{n^{(1)}} in terms of the states \ket{m^{(0)}}: \ket{n^{(1)}} = \sum_{m \neq n} c_m \ket{m^{(0)}}. Note that we explicitly exclude the state \ket{n^{(0)}} from the sum, i.e. we take \left\langle n^{(1)} | n^{(0)} \right\rangle = 0. Fixing the normalization of our states \left\langle n | n \right\rangle = 1 leads to this condition as a requirement; we can also notice that from the first-order equation above, we can add an arbitrary amount of \ket{n^{(0)}} to \ket{n^{(1)}} without affecting the equation.
Show that the condition of normalizing the states \left\langle n | n \right\rangle = 1 requires us to have \left\langle n^{(1)} | n^{(0)} \right\rangle = 0.
Answer:
Just taking the expansion of \ket{n} defined above and plugging it in to \left\langle n | n \right\rangle, we start organizing terms order by order in \lambda: \left\langle n | n \right\rangle = \left\langle n^{(0)} | n^{(0)} \right\rangle + \lambda \left[\left\langle n^{(1)} | n^{(0)} \right\rangle + \left\langle n^{(0)} | n^{(1)} \right\rangle \right] \\ + \lambda^2 \left[ \left\langle n^{(0)} | n^{(2)} \right\rangle + \left\langle n^{(1)} | n^{(1)} \right\rangle + \left\langle n^{(2)} | n^{(0)} \right\rangle \right] + ... If we want to have properly normalized states both before and after perturbation, i.e. \left\langle n | n \right\rangle = \left\langle n^{(0)} | n^{(0)} \right\rangle = 1, then all of the correction terms have to vanish identically, i.e. \sum_{k=0}^j \left\langle n^{(j-k)} | n^{(k)} \right\rangle = 0 from order \lambda^j in the expansion above. So things get more complicated at second order (where the condition relates \left\langle n^{(0)} | n^{(2)} \right\rangle to -\left\langle n^{(1)} | n^{(1)} \right\rangle.)
Although we can impose this normalization condition on the states as we go, it’s not vital that we do so; in fact, in a later section below we’re going to choose a completely different condition that makes the state \ket{n} improperly normalized. There is nothing wrong with using unnormalized states: when we come to calculating probabilities or expectation values, we can always just “renormalize” the state by dividing out the norm, \ket{\psi}_N = \frac{1}{\left\langle \psi | \psi \right\rangle} \ket{\psi}.
Going back to our first-order equation and using this expansion gives \sum_{m \neq n} (\hat{H}_0 - E_n^{(0)}) c_m \ket{m^{(0)}} = (E_n^{(1)} - \hat{W}) \ket{n^{(0)}} \\ \sum_{m \neq n} (E_m^{(0)} - E_n^{(0)}) c_m \ket{m^{(0)}} = (E_n^{(1)} - \hat{W}) \ket{n^{(0)}}. If we multiply on the left by another arbitrary eigenstate \ket{k^{(0)}}, we can collapse the sum: \sum_{m \neq n} \bra{k^{(0)}} c_m (E_k^{(0)} - E_n^{(0)}) \ket{m^{(0)}} = E_n^{(1)} \delta_{kn} - \bra{k^{(0)}} \hat{W} \ket{n^{(0)}} Setting k=n make the left side vanish, and just returns our expression for E_n^{(1)}. If instead we require k \neq n, then the E_n^{(1)} term vanishes, the sum on the left collapses, and we end up with the expression c_k = \frac{\bra{k^{(0)}} \hat{W} \ket{n^{(0)}}}{E_n^{(0)} - E_k^{(0)}} and so \ket{n^{(1)}} = \sum_{k \neq n} \ket{k^{(0)}} \frac{W_{kn}}{E_n^{(0)} - E_k^{(0)}}. where I’ve introduced a more compact notation for the matrix elements of \hat{W} with respect to the unperturbed basis, W_{kn} = \bra{k^{(0)}} \hat{W} \ket{n^{(0)}}. Don’t forget that the order matters, in general; W_{nk} = W_{kn}^\star, so flipping the order of the indices comes with a complex conjugation. (Typically \hat{W} will be Hermitian since it’s part of a Hamiltonian, in which case this doesn’t matter.)
It’s easy to remember that the sum above is only over k \neq n; if you try to include that term, it will blow up! Similarly, if any of the energy eigenvalues are degenerate, the sum also blows up, so you can see why we had to make that assumption in deriving this. (If there is such a degeneracy, then the derivation above gives us no information on the corresponding c_k, so we’ll have to find another approach.)
20.1.2 Going to second order in energy
In general, the eigenstates at order n-1 are all we need to compute the energies at order n. We won’t keep going indefinitely, but now that we’ve found the first-order eigenstates, it’s easy enough to continue and get the energy corrections at second order. Returning to our above equation for the second-order perturbation and again dotting in \ket{n^{(0)}} on the left, we see that E_n^{(2)} = \bra{n^{(0)}} (\hat{W} - E_n^{(1)}) \ket{n^{(1)}}. The second term is zero since as we argued above, \ket{n^{(0)}} and \ket{n^{(1)}} are orthogonal. Expanding the first term out, E_n^{(2)} = \bra{n^{(0)}} \hat{W} \left( \sum_{k \neq n} \ket{k^{(0)}} \frac{\bra{k^{(0)}} \hat{W} \ket{n^{(0)}}}{E_n^{(0)} - E_k^{(0)}} \right) \\ = \sum_{k \neq n} \frac{|W_{nk}|^2}{E_n^{(0)} - E_k^{(0)}}.
An easy way to remember the sign of the denominator in the second-order correction is to remember the fact that the correction to the ground-state energy n=0 is always negative. This is not simply a coincidence; if we look at the expectation value of the full Hamiltonian with respect to the zeroth-order eigenstates, we see \bra{n^{(0)}} \hat{H} \ket{n^{(0)}} = \sum_{m} \left\langle n^{(0)} | m \right\rangle \bra{m} \hat{H} \ket{m} \left\langle m | n^{(0)} \right\rangle, where \ket{m} are the eigenstates of the full Hamiltonian. But that means that the matrix elements are just the true energy values of \hat{H}: \bra{n^{(0)}} \hat{H} \ket{n^{(0)}} = \sum_m E_m |\left\langle n^{(0)} | m \right\rangle|^2. The squared inner products under the sum are exactly the probabilities for observing zeroth-order eigenstate \ket{n^{(0)}} in full eigenstate \ket{m}; the sum over all of them will be 1, in particular. Meanwhile, the left-hand side is just the first-order energy, i.e. E_n^{(0)} + \lambda \bra{n^{(0)}} \hat{W} \ket{n^{(0)}} = \sum_m E_m |\left\langle n^{(0)} | m \right\rangle|^2. Since the first-order energy eigenvalues are given by a weighted sum over the full energy eigenvalues, with positive weights, we conclude that the distribution of all the approximate energy eigenvalues must be more closely spaced than that of the true eigenvalues. This explains is why the next correction is negative for the ground-state; it tends to increase the overall spread in all energy eigenvalues as the true solution is approached.
Let’s collect our results together and summarize:
Given a Hamiltonian of the form \hat{H} = \hat{H}_0 + \hat{W}, in terms of the energy eigenvalues \ket{E_n^{(0)}} of \hat{H}_0, the corrected energy to second order in perturbation theory is given by E_n = E_n^{(0)} + \lambda W_{nn} + \lambda^2 \sum_{k \neq n} \frac{|W_{nk}|^2}{E_n^{(0)} - E_k^{(0)}} + ... while the first-order corrected energy eigenstates are \ket{n} = \ket{n^{(0)}} + \sum_{k \neq n} \frac{\lambda W_{kn}}{E_n^{(0)} - E_k^{(0)}} + ...
The sign of the total second-order energy correction is always negative.
20.2 Controlled example: the harmonic oscillator
The best way to test out a new method is of course to apply it to a system where you already know the answer. So let’s take the one-dimensional harmonic oscillator, and add a “perturbation” to it of the same quadratic form: \hat{H} = \frac{\hat{p}{}^2}{2m} + \frac{1}{2} m\omega^2 \hat{x}^2 + \lambda \hat{x}^2. Of course, we don’t need perturbation theory to deal with this: for any value of \lambda, we can just rewrite this as W(\hat{x}) = \frac{1}{2} m (\omega^2 + \frac{2\lambda}{m}) \hat{x}^2 = \frac{1}{2} m \omega'{}^2 \hat{x}^2 where \omega' = \sqrt{\omega^2 + \frac{2\lambda}{m}} \approx \omega + \frac{\lambda}{m \omega} - \frac{\lambda^2}{2m^2 \omega^3} + ... Note that \lambda has the same units as m \omega^2, which is a quick check that our expansion is sensible.
Now, we go back and identify \hat{H}_0 = \frac{\hat{p}{}^2}{2m} + \frac{1}{2} m\omega^2 \hat{x}^2, \\ \lambda \hat{W} = \lambda \hat{x}^2 The unperturbed (zeroth-order) energies are then E_{n}^{(0)} = \hbar \omega \left(n + \frac{1}{2} \right), corresponding to the standard eigenstates \ket{n^{(0)}}. The first-order energy correction is then E_n^{(1)} = \bra{n^{(0)}} \hat{x}{}^2 \ket{n^{(0)}}. We can work this out using the ladder operators: E_n^{(1)} = \frac{\hbar}{2m\omega} \bra{n^{(0)}} (\hat{a} + \hat{a}^{\dagger})^2 \ket{n^{(0)}} \\ = \frac{\hbar}{2m\omega} \bra{n^{(0)}} (\hat{a} \hat{a}^{\dagger} + \hat{a}^\dagger \hat{a}) \ket{n^{(0)}} \\ = \frac{\hbar}{m\omega} (n+\frac{1}{2}). Plugging this back in above, we find that E_n \approx E_n^{(0)} + \lambda E_n^{(1)} = \hbar \left( \omega + \frac{\lambda}{m \omega} \right) \left(n + \frac{1}{2}\right), matching our first-order expansion of \omega' above.
Note that we’re using the ladder operators constructed explicitly for the unperturbed Hamiltonian \hat{H}_0 here. There’s nothing inconsistent about this approach, since they were defined explicitly from a rearrangement of the simple harmonic oscillator Hamiltonian, \hat{H}_0. In general, when we add a perturbation we will simply find that \hat{N} = \hat{a}^\dagger \hat{a} defined in this way will still commute with \hat{H}_0, but not with the full \hat{H}. In our current example, it’s easy to see that \hat{a}^\dagger \hat{a} = \frac{\hat{H}_0}{\hbar \omega} - \frac{1}{2} \\ \Rightarrow \hat{H} = \hbar \omega (\hat{N} + \frac{1}{2}) + \lambda \hat{x}^2 and [\hat{N}, \hat{x}^2] = \frac{\hbar}{2m\omega}[\hat{N}, (\hat{a} + \hat{a}^{\dagger})^2] \\ = \frac{\hbar}{2m\omega} \left( [\hat{N}, \hat{a} + \hat{a}^\dagger] (\hat{a} + \hat{a}^\dagger) + (\hat{a} + \hat{a}^\dagger) [\hat{N}, \hat{a} + \hat{a}^\dagger) \right) \\ = \frac{\hbar}{2m\omega} \left( (-\hat{a} + \hat{a}^\dagger) (\hat{a} + \hat{a}^\dagger) + (\hat{a} + \hat{a}^\dagger) (-\hat{a} + \hat{a}^\dagger) \right) \\ = \frac{\hbar}{m\omega} \left( \hat{a}^\dagger \hat{a}^\dagger - \hat{a} \hat{a} \right) which is clearly not zero. In this special case, we can of course define new ladder operators \hat{a}' which use the modified frequency \omega', and then everything will commute and we can label the states of the full Hamiltonian with the new number operator \hat{N}'. However, this is only because we can map the perturbative problem back onto the full SHO here. Normally, when we deal with perturbations, we want to work in terms of the operators and states of the analytically solvable \hat{H}_0 only.
Let’s keep going, and look at the second-order energy shift. To evaluate it, we will need the matrix elements of \hat{W} between arbitrary zeroth-order energy eigenstates. We already evaluated the diagonal case, now we need the explicitly off-diagonal terms k \neq n: \bra{k^{(0)}} \hat{W} \ket{n^{(0)}} = \frac{\hbar}{2m\omega} \bra{k^{(0)}} \hat{a}{}^2 + (\hat{a}^\dagger)^2 \ket{n^{(0)}} \\ = \frac{\hbar}{2m\omega} \left[ \sqrt{n(n-1)} \delta_{k,n-2} + \sqrt{(n+1)(n+2)} \delta_{k,n+2} \right] So the sum over all states for the second-order energy correction collapses: E_n^{(2)} = \sum_{k \neq n} \frac{|W_{nk}|^2}{E_n^{(0)} - E_k^{(0)}} \\ = \frac{\hbar^2}{4m^2 \omega^2} \left[ \frac{n(n-1)}{\hbar \omega (n - (n-2))} + \frac{(n+1)(n+2)}{\hbar \omega (n-(n+2))}\right] \\ = \frac{\hbar}{8m^2 \omega^3} \left[ n(n-1) - (n+1)(n+2) \right] \\ = \frac{-\hbar}{2m^2 \omega^3} \left(n + \frac{1}{2} \right). Multiplying by \lambda^2, we see that our second-order formula becomes E_n \approx E_n^{(0)} + \lambda E_n^{(1)} + \lambda^2 E_n^{(2)} \\ = \hbar \left( \omega + \frac{\lambda}{m \omega} - \frac{\lambda^2}{2m^2 \omega^3} \right) \left(n + \frac{1}{2} \right), again nicely matching our series expansion of the exact solution.
20.2.1 Another controlled example: the shifted harmonic oscillator
We could use the formula we derived above to solve for the eigenstates at first order as well; but since the position wavefunction of the harmonic oscillator doesn’t have a simple dependence on the parameter \omega, it would be a rather complicated exercise. Instead, let’s turn to another form of perturbed simple harmonic oscillator: V(\hat{x}) = \frac{1}{2} m \omega^2 \hat{x}^2 + \lambda \hat{x}. so \hat{W} = \hat{x}. Notice that once again, we can turn this into an ordinary SHO with a change of variables: completing the square, the potential above is equal to V(\hat{x}) = \frac{1}{2} m \omega^2 \left( \hat{x} + \frac{\lambda}{m \omega^2} \right)^2 - \frac{\lambda^2}{2m \omega^2} This rewritten harmonic oscillator has the same frequency \omega, but its position has been translated by -\lambda / (m \omega^2), and the overall energy has been shifted by -\lambda^2 / (2m\omega^2). The energy levels of the SHO don’t depend on where the origin is, so the energy shift is exactly proportional to \lambda^2 and independent of n: \Delta E_n = -\frac{\lambda^2}{2m \omega^2}. This means that if we attempt to calculate the perturbation at first order, we should find zero, and indeed we see immediately that E_n^{(1)} = \sqrt{\frac{\hbar}{2m\omega}} \bra{n^{(0)}} (\hat{a} + \hat{a}^\dagger) \ket{n^{(0)}} = 0 for all n. Proceeding to second order, the off-diagonal matrix elements are \bra{k^{(0)}} \hat{W} \ket{n^{(0)}} = \sqrt{\frac{\hbar}{2m\omega}} \bra{k^{(0)}} (\hat{a} + \hat{a}^\dagger) \ket{n^{(0)}} \\ = \sqrt{\frac{\hbar}{2m\omega}} (\sqrt{n} \delta_{k,n-1} + \sqrt{n+1} \delta_{k,n+1}).
so \lambda^2 E_n^{(2)} = \sum_{k \neq n} \frac{|W_{nk}|^2}{E_n^{(0)} - E_k^{(0)}} \\ = \frac{\hbar \lambda^2}{2m\omega} \left(\frac{n}{\hbar \omega (n-(n-1))} + \frac{n+1}{\hbar \omega (n-(n+1))} \right) \\ = \frac{\lambda^2}{2m\omega^2} (n-(n+1)) = -\frac{\lambda^2}{2m\omega^2} independent of n as expected.
Now let’s look at the energy eigenstates themselves; since our perturbation can be rewritten as a translation, we expect a relatively simple change to the energy eigenstate wavefunctions. Let’s focus on the ground state \ket{0} explicitly: we have for the unperturbed Hamiltonian \left\langle x | 0^{(0)} \right\rangle = \left(\frac{m \omega}{\pi \hbar}\right)^{1/4} \exp \left[ -\frac{m\omega}{2\hbar} {x}^2 \right]. After shifting, the eigenstate of the full Hamiltonian is \left\langle x | 0 \right\rangle = \left( \frac{m \omega}{\pi \hbar}\right)^{1/4} \exp \left[ -\frac{m \omega}{2\hbar} \left({x} + \frac{\lambda}{m\omega^2}\right)^2 \right] \\ \approx \left( \frac{m \omega}{\pi \hbar}\right)^{1/4} \left(1 - \frac{\lambda {x}}{\hbar \omega} + ... \right) \exp \left[ -\frac{m\omega}{2\hbar} {x}^2 \right] at first order in \lambda. Recalling the perturbative formula for the first-order eigenket, we have \left\langle x | 0^{(1)} \right\rangle = \sum_{k \neq n} \ket{k^{(0)}} \frac{W_{kn}}{-\hbar \omega k} \\ = -\left\langle x | 1^{(0)} \right\rangle \frac{1}{\omega} \sqrt{\frac{1}{2m\hbar\omega}} \\ = -\frac{1}{\omega} \sqrt{\frac{1}{2m\hbar \omega}} \left( \frac{m \omega}{\pi \hbar}\right)^{1/4} \sqrt{\frac{2m\omega}{\hbar}} {x} \exp \left[ -\frac{m \omega}{2\hbar} {x}^2 \right] \\ = - \left(\frac{m \omega}{\pi \hbar}\right)^{1/4} \frac{x}{\hbar \omega} \exp \left[ - \frac{m \omega}{2\hbar} x^2 \right], exactly matching the expansion above for \ket{0} \approx \ket{0^{(0)}} + \lambda \ket{0^{(1)}} + ...
Now that we’ve had some practice, let’s move on to a more formal treatment of perturbation theory.
20.3 Formal solution: Brillouin-Wigner perturbation theory
This approach is mainly a different way to write the expansion we’ve done already, but it’s a way which makes continuing to higher orders in the expansion more systematic. The eigenvalue equation that we’re trying to solve can be written in the form (\hat{H}_0 - E_n) \ket{n} = -\lambda \hat{W} \ket{n}. or rearranging slightly, (E_n^{(0)} - \hat{H}_0) \ket{n} = (\lambda \hat{W} - \Delta_n) \ket{n}, where we’ve defined \Delta_n = E_n - E_n^{(0)}, the total energy shift due to the perturbation. Now, if we can find an inverse to the operator on the left-hand side, then the solution can be written \ket{n} = (E_n^{(0)} - \hat{H}_0)^{-1} (\lambda \hat{W} - \Delta_n) \ket{n}. We could then solve this equation iteratively; starting with \ket{n^{(0)}} on the right will determine the order-\lambda correction \ket{n^{(1)}}, and so forth. However, even assuming that all of the energy states are non-degenerate, we still have a problem; the inverse operator doesn’t exist because its action on the state \ket{n^{(0)}} is undefined.
To escape this dilemma, we define the complementary projection operator to state \ket{n^{(0)}}: \hat{\phi}_n \equiv 1 - \ket{n^{(0)}}\bra{n^{(0)}} = \sum_{k \neq n} \ket{k^{(0)}}\bra{k^{(0)}}. The combination of this operator with the inverse above is always (with no degeneracy in the energies) well-defined: \frac{1}{E_n^{(0)} - \hat{H}_0} \hat{\phi}_n = \sum_{k \neq n} \frac{1}{E_n^{(0)} - E_k^{(0)}} \ket{k^{(0)}} \bra{k^{(0)}}. Applying this projected operator to (\lambda \hat{W} - \Delta_n) will then give us the general solution to the equation above…almost. If we blindly plug back in to the equation above, we will find that \ket{n} is given by a superposition of unperturbed eigenstates \ket{k^{(0)}}, but only for k \neq n. But we know that in the limit \lambda \rightarrow 0, we should recover \ket{n} = \ket{n^{(0)}}. So we know that \ket{n} must have some component in the \ket{n^{(0)}} direction, which means we should re-write the equation for \ket{n} as \ket{n} = c_n(\lambda) \ket{n^{(0)}} + \frac{1}{E_n^{(0)} - \hat{H}_0} \hat{\phi}_n (\lambda \hat{W} - \Delta_n) \ket{n}, where c_n(\lambda) is an unknown function, but we know that it must satisfy \lim_{\lambda \rightarrow 0} c_n(\lambda) = 1. Fortunately, we only missed the coefficient of one of the unperturbed energy eigenstates, which means that the value of c_n(\lambda) is simply determined by the normalization condition for this state.
We could impose the standard normalization \left\langle n | n \right\rangle = 1, as we saw before, but a more convenient choice turns out to be c_n(\lambda) = \left\langle n^{(0)} | n \right\rangle = 1 for all \lambda. We’re free to do this because c_n(\lambda) is determined only by the overall normalization of our state; any other choice of c_n(\lambda) simply amounts to a multiplicative rescaling of all of the unperturbed eigenstates, which doesn’t affect any of our physical results (as long as we enforce normalization when we compute probabilities.)
With this normalization condition (known as intermediate normalization) imposed, our equation for the state \ket{n} becomes \ket{n} = \ket{n^{(0)}} + \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} (\lambda \hat{W} - \Delta_n) \ket{n}.
In general, the notation above where we write a fraction involving two operators is ambiguous, because the order of operators matters. In this case, we are free to adopt the fraction notation since the projection operator can be applied on either side of the inverse of E_n^{(0)} - \hat{H}_0: \frac{1}{E_n^{(0)} - \hat{H}_0} \hat{\phi}_n = \hat{\phi}_n \frac{1}{E_n^{(0)} - \hat{H}_0} = \hat{\phi}_n \frac{1}{E_n^{(0)} - \hat{H}_0} \hat{\phi}_n.
Before we do anything else, notice that this formulas makes a key feature of perturbation theory very obvious: the perturbed state \ket{n} can be decomposed into two vector components, one pointing along the original ket \ket{n^{(0)}}, and one orthogonal component which is small (because it’s multiplied by the perturbation and by \lambda.)
Now, if we just plug in and match terms order by order, we’ll end up back at the same expressions we had for the Rayleigh-Schrödinger expansion above. Instead, we note that because this equation is exactly correct, we can plug it into itself. First for compactness, we define \lambda \hat{R} \equiv \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} (\lambda \hat{W} - \Delta_n). Then we have \ket{n} = \ket{n^{(0)}} + \lambda \hat{R} \ket{n} \\ = \ket{n^{(0)}} + \lambda \hat{R} \left( \ket{n^{(0)}} + \lambda \hat{R} \ket{n} \right) \\ = \ket{n^{(0)}} + \lambda \hat{R} \ket{n^{(0)}} + \lambda^2\hat{R}^2 \ket{n^{(0)}} + ...
from which we simply read off that \ket{n^{(j)}} = \hat{R}^{j} \ket{n^{(0)}}. What about the energies? If we go all the way back to the second equation in this section and multiply on the left by \bra{n^{(0)}}, the left-hand side vanishes, and we find \bra{n^{(0)}} (\lambda \hat{W} - \Delta_n) \ket{n} = 0. Now you can see why the intermediate normalization comes in handy; with \left\langle n^{(0)} | n \right\rangle = 1, the energy shift \Delta_n is simply given by \Delta_n = \bra{n^{(0)}} \lambda \hat{W} \ket{n}. All of this nicely sets up the possibility of an iterative solution. We start with \ket{n} = \ket{n^{(0)}} as our initial estimate, plug in to the equations above for the state and the energy and get the first-order solution, then plug that in and keep iterating until we get to the desired order.
This framework also gives us a systematic way to recover the Rayleigh-Schrödinger series in a much more systematic way, simply by plugging in for both \ket{n} and E_n as power series in \lambda and then matching order by order.
Notice that the series expansion for \Delta_n starts at order \lambda:
\Delta_n = \lambda E_n^{(1)} + \lambda^2 E_n^{(2)} + ...
If we set \ket{n} =\ket{n^{(0)}} and \Delta_n = 0 on the right-hand side of the equation for \ket{n} above, we find
\ket{n} \approx \ket{n^{(0)}} + \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} (\lambda \hat{W}) \ket{n^{(0)}}
and so
\ket{n^{(1)}} = \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} \hat{W} \ket{n^{(0)}} \\
= \sum_{k \neq n} \frac{1}{E_n^{(0)} - E_k^{(0)}} \ket{k^{(0)}} \bra{k^{(0)}} \hat{W} \ket{n^{(0)}},
recovering the formula we derived in Rayleigh-Schrödinger above.
At second order, we see that
E_n^{(2)} = \bra{n^{(0)}} \hat{W} \ket{n^{(1)}} \\
= \bra{n^{(0)}} \hat{W} \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} \hat{W} \ket{n^{(0)}}.
Higher order states and energies can be built up in succession just by applying more copies of the projection operator and \hat{W}. In general, we can write the equation for \ket{n} as a power series on both sides:
\ket{n^{(1)}} + \lambda \ket{n^{(2)}} + ... = \\
\frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} (\hat{W} - E_n^{(1)} - \lambda E_n^{(2)} - ...)
\times (\ket{n^{(0)}} + \lambda \ket{n^{(1)}} + ...)
where I’ve cancelled the zeroth-order terms. We can read off the second-order correction to the state:
\ket{n^{(2)}} = \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} \left[ (\hat{W} - E_n^{(1)}) \ket{n^{(1)}} - E_n^{(2)} \ket{n^{(0)}} \right] \\
= \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} \hat{W} \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} \hat{W} \ket{n^{(0)}} \\
- \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} \bra{n^{(0)}} \hat{W} \ket{n^{(0)}} \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} \hat{W} \ket{n^{(0)}}
where the third term which would arise from E_n^{(2)} vanishes, since it contains the projector \hat{\phi}_n acting directly on the state \ket{n^{(0)}}.
How about the third-order state? This would be an enormous pain to work out in the Rayleigh-Schrödinger approach we used last time, but in this operator construction we can write it down immediately: \ket{n^{(3)}} = \frac{\hat{\phi}_n}{E_n^{(0)} - \hat{H}_0} \left[ (\hat{W} - E_n^{(1)}) \ket{n^{(2)}} - E_n^{(2)} \ket{n^{(1)}} \right], and we just plug in all of our lower-order expansions from above to express this in terms of operators acting on \ket{n^{(0)}}.
All of the expansions we’ve now studied above, particularly the one for the states, make it clear that the rigorous criterion for convergence of the perturbative series is that the matrix elements in the numerator should be small compared to the energy differences in the denominator, i.e. we should have in general \frac{\lambda |W_{ij}|}{E_i^{(0)} - E_j^{(0)}} \ll 1 for rapid convergence. (To actually derive a radius of convergence, we would need to know more about the system that we’re studying.)
20.3.1 Wavefunction renormalization
As we have noted, the states \ket{n} derived in this approach are not properly normalized. We can formally “renormalize” the states just by rescaling them: \ket{n}_N = Z_n^{1/2} \ket{n}, where Z_n = 1/\left\langle n | n \right\rangle. Notice that from the equation above, we also have the condition Z_n^{1/2} = \left\langle n^{(0)} | n \right\rangle_N. With the state \ket{n}_N now properly normalized, we can interpret Z_n itself as a probability of finding the perturbed eigenstate \ket{n}_N in the associated unperturbed state \ket{n^{(0)}}. We can clarify this somewhat by expanding in \lambda: Z_n^{-1} = \left\langle n | n \right\rangle = \left\langle n^{(0)} | n^{(0)} \right\rangle + \lambda^2 \left\langle n^{(1)} | n^{(1)} \right\rangle + ... where I’ve once again invoked the fact that \left\langle n^{(0)} | n^{(1)} \right\rangle = 0. Inverting and plugging in the expression for \ket{n^{(1)}} from above, we have Z_n = 1 - \lambda^2 \sum_{k \neq n} \frac{|W_{kn}|^2}{(E_n^{(0)} - E_k^{(0)})^2} + ... The renormalization factor Z_n is less than 1, as it must be since we have a probabilistic interpretation of it. For non-zero \lambda, the probability of observing state \ket{n^{(0)}} is reduced, with the second term representing the sum over probabilities to find any other state as an outcome.
20.4 More examples of time-independent perturbation theory
This has been very formal, so let’s digress and see a couple more examples of how perturbation theory is used in practice. Let’s start with a simple example, but one where we don’t have an exact result to compare to
20.4.1 Example: particle in a warped box
For this first example, we start with the one-dimensional particle in a box, V(x) = \begin{cases} 0, & 0 \leq x \leq a; \\ \infty, & {\rm elsewhere}. \end{cases} Now we add a perturbation of the form \lambda W(x) = \lambda E' \cos \left( \frac{2\pi x}{a} \right) where now \lambda is dimensionless.
The eigenenergies of the unperturbed system are E_n^{(0)} = \frac{\hbar^2 \pi^2 n^2}{2ma^2} = n^2 E_1^{(0)} and the states can be written simply as \left\langle x | n^{(0)} \right\rangle = \sqrt{\frac{2}{a}} \sin \left( \frac{n\pi x}{a} \right). Let’s evaluate the general matrix elements of \hat{W} with respect to the unperturbed energy eigenstates; since we’re working in position space this will require an integral, W_{mn} = \bra{m^{(0)}} \hat{W} \ket{n^{(0)}} \\ = \frac{2E'}{a} \int_0^a dx\ \sin \left( \frac{m\pi x}{a} \right) \cos \left( \frac{2\pi x}{a} \right) \sin \left( \frac{n \pi x}{a} \right). To evaluate, we can make use of the product-sum formula \sin u \cos v = \frac{1}{2} \left[ \sin (u+v) + \sin (u-v) \right], giving W_{mn} = \frac{E'}{a} \int_0^a dx\ \sin \left( \frac{m \pi x}{a} \right) \left[ \sin \left( \frac{(n+2) \pi x}{a} \right) + \sin \left( \frac{(n-2) \pi x}{a} \right) \right] \\ = \frac{E'}{2} (\delta_{m,n+2} + \delta_{m,n-2}) using the orthogonality of the trig functions on this interval. Actually, this formula is only valid for n > 2; we have to treat n=1 and n=2 separately. It turns out that the product-sum formula actually works for negative values too, but we want to express everything in terms of the states with n>1. For n=2, the second term in the formula goes to \sin ((2-2) \pi x/a) = 0, so we simply have W_{m2} = \frac{E'}{2} \delta_{m, 4}. For n=1, the second term gives instead \sin(-\pi x/a) = - \sin(\pi x/a), and the integration gives W_{m1} = \frac{E'}{2} (\delta_{m,3} - \delta_{m,1}). The first-order energy shifts are given by the diagonal components of this matrix; we see that the only state which overlaps with itself is n=1, and thus E_n^{(1)} = \begin{cases} -E'/2, & n=1; \\ 0, & n>1. \end{cases} The corrected kets at first order are easy to work out as well: \ket{n^{(1)}} = \sum_{k \neq n} \ket{k^{(0)}} \frac{W_{kn}}{E_n^{(0)} - E_k^{(0)}} \\ = \frac{E'}{2E_1} \left[ \frac{\ket{(n-2)^{(0)}}}{n^2 - (n-2)^2} + \frac{\ket{(n+2)^{(0)}}}{n^2 - (n+2)^2} \right] or for the special cases n \leq 2, \ket{2^{(1)}} = -\frac{E'}{24E_1} \ket{4^{(0)}}, \\ \ket{1^{(1)}} = -\frac{E'}{16E_1} \ket{3^{(0)}}.
20.4.2 Example: nucleus finite-size correction
Let’s try a more real-world example. For a hydrogenic atom (one electron), we assume that the electron moves subject to the Coulomb potential, V(r) = -\frac{Ze^2}{r}. where Ze is the charge of the nucleus. This potential is valid at large r, but if we look at very short distances, we will realize that the nucleus is not a point charge, but has some finite size associated with it. If we treat the nucleus as a uniform ball of charge with finite radius r_0, then the potential is corrected to V(r) = \begin{cases} -\frac{Ze^2}{r}, & r \geq r_0; \\ \frac{Ze^2}{2r_0} \left[ \left( \frac{r}{r_0}\right)^2 - 3 \right], & r \leq r_0. \end{cases} We can thus rewrite our Hamiltonian as \hat{H} = \hat{H}_0 + \lambda \hat{W}, where \hat{W} = \begin{cases} 0, & r \geq r_0; \\ \frac{Ze^2}{2r_0} \left[ \left( \frac{r}{r_0}\right)^2 - 3 + \frac{2r_0}{r} \right], & r \leq r_0. \end{cases} where this time I am using \lambda purely as an auxiliary constant, which we’ll set to 1 at the end. Why is it valid to treat this correction perturbatively? We know that the nucleus is very small, and so we expect that the overlap of the electron wavefunction with the region we’re changing the potential in is quite small. We will find that the matrix elements of \hat{W} depend on the ratio r_0 / a_0, where a_0 is the (rescaled) Bohr radius a_0 = \frac{\hbar^2}{Z m_e e^2}. The small size of this ratio validates our perturbative expansion; for hydrogen, a_0 \sim 10^{-10} m, while the charge radius of the proton is closer to 10^{-15} m.
To compute the energy shift due to this finite-size effect, we have to integrate using hydrogenic wavefunctions. Let’s just look at the ground state \ket{n,l,m} = \ket{1,0,0}. Here the angular dependence vanishes, and we simply have \psi_{100}(r, \theta, \phi) = \frac{2}{\sqrt{4\pi a_0^3}} e^{-r/a_0}. The first-order energy correction is then E_{100}^{(1)} = \bra{100} \hat{W} \ket{100} \\ = \int_0^{r_0} dr\ r^2 \frac{2Z e^2}{r_0 a_0^3} \left[ \left( \frac{r}{r_0} \right)^2 - 3 + \frac{2r_0}{r} \right] e^{-2r/a_0}. (the 4\pi vanished due to the angular integral we skipped over.) This is a straightforward integral, which you can do in Mathematica or with integral tables. The result, as promised, depends on the small ratio r_0 / a_0. Although taken at face value, the energy shift will contain higher powers of r_0 / a_0, we know that we’re only working at first order here we should really only keep the first term in a series expansion, which gives E_{100}^{(1)} = \frac{2Ze^2 r_0^2}{5 a_0^3} = \frac{4}{5} E_{100}^{(0)} \left( \frac{r_0}{a_0}\right)^2. The correction is therefore of order (r_0 / a_0)^2, i.e. a 10^{-10} effect for the hydrogen atom; very small!
This is, in fact, the largest correction factor out of all of the states of hydrogen; we know that for higher n or higher l, the wavefunction of the electron will be localized further away from the origin, and they will turn out to be suppressed by further powers of (r_0 / a_0). This is a tiny effect, but can be relevant in very precise atomic physics experiments. In particular, for higher-Z nuclei the nuclear radius increases, while the Bohr radius decreases, and so the size of the perturbation can be substantially larger. Replacing the electron with a muon, the electron’s heavier cousin, can also lead to appreciable effects by greatly reducing the orbital radius.
20.4.3 Example: the quadratic Stark effect
Let’s consider another example of a perturbative effect in a hydrogenic atom, this time in response to a uniform external electric field \vec{E}. This time, the perturbation is \vec{E} itself, which we assume that we’re controlling with a dial in our laboratory so that it remains “small”, in the sense that matrix elements involving \vec{E} will give small energy corrections relative to the existing energy splittings of our atom. This means we don’t have to bother with a dimensionless \lambda at all, although if you want to be more neat about units, you could introduce it and write something like \vec{E} = \lambda \vec{E}_0, where maybe \vec{E}_0 is the full-strength magnetic field in your lab for example. We thus write our perturbation as \hat{W} = -e|\textbf{E}| \hat{z}.
There is a subtlety to this particular perturbation; because it is linear in z, we find \hat{W} \rightarrow -\infty as z \rightarrow -\infty. This means that technically, it’s possible for our formerly bound electron to escape the hydrogen atom completely due to the applied field. We can ignore this complication, partly because in reality the electric field doesn’t actually have infinite extent and in the region where it is present, it won’t be strong enough for appreciable tunnelling to happen if |\vec{E}| is small enough.
We don’t have to ignore this effect either - another application of perturbation theory here that could be interesting is calculating the lifetime of the electron states with non-zero |\mathbf{E}|.
In general, the treatment of an electric-field perturbation on the hydrogen atom requires degenerate-state perturbation theory, since all of the hydrogen energy levels are n^2-fold degenerate. However, we’re going to focus explicitly just on the 1s ground state (and ignore the electron spin), which means that we can treat this as a non-degenerate perturbation theory problem since there is only one state to consider.
The other reason to focus on the ground state is the particular physics of the quadratic Stark effect, which we can understand just in terms of classical electrodynamics. A hydrogen atom is a neutral object, with no net charge for an applied electric field to couple to. The next leading possible effect for an object without charge would be an electric dipole, which would couple to the electric field linearly. But in the 1s state, the electron orbital is perfectly spherical, which means it can’t have a dipole moment (or any higher charge moment, for that matter.)
However, this doesn’t mean that the electric field has no effect at all. The hydrogen atom isn’t just a neutral sphere, it consists of charged components. In the presence of an electric field, the proton will feel a force along the direction of \vec{E}, while the electron will be pulled in the opposite direction. This induces a dipole moment, by deforming the charge distribution along the direction of the electric field, and the resulting dipole moment then couples to the electric field. This is why the effect is quadratic: the interaction in the Hamiltonian has the form \vec{d} \cdot \vec{E}, but \vec{d} is itself proportional to \vec{E}.
We’ll see this happen explicitly just arguing based on parity in a moment, but the underlying physics can be helpful to keep in mind. We would expect that when we consider excited states of hydrogen with angular momentum, there will be dipole moments present as the orbitals become more complicated, and we’ll reveal a linear coupling to \vec{E} (the linear Stark effect, which we will see soon.)
To second order, the correction to the ground-state energy is given by \Delta = -e |\mathbf{E}| z_{11} + e^2 |\mathbf{E}|^2 \sum_{j \neq 1} \frac{|z_{1j}|^2}{E_1^{(0)} - E_j^{(0)}} where the indices j are standing in for the triplet of labels \ket{nlm}, with index 1 corresponding to the ground state \ket{100}.
The first term in this expansion, \left\langle \hat{z} \right\rangle, is easily seen to vanish identically; although the electric field breaks parity, the original Hamiltonian is parity invariant, and so the ground state \ket{1^{(0)}} is a parity eigenstate. Since the operator \hat{z} is itself odd, it’s easy to see that \left\langle \hat{z} \right\rangle = 0. Physically, this corresponds to the fact that the spherically-symmetric hydrogen atom does not have an electric dipole moment, as we argued before.
However, as we mentioned the electric field itself breaks parity, causing polarization of the charge distribution. This effect results in an energy shift which depends quadratically on the electric field, since the induced dipole moment is itself proportional to |\mathbf{E}|: \Delta E = -\frac{1}{2} \alpha |\mathbf{E}|^2. The coefficient \alpha is the electric polarizability, which we’ve encountered before briefly in the ammonia maser example. The resulting energy splitting is known as the quadratic Stark effect, quadratic for obvious reasons and Stark after Johannes Stark, who discovered it back in 1913. We can read the polarizability of the atomic ground state off from the expansion above: \alpha = -2e^2 \sum_{k>1} \frac{ |\bra{k^{(0)}} \hat{z} \ket{1^{(0)}}|^2 }{E_1^{(0)} - E_k^{(0)}} where again, the label 1 is standing in for unperturbed state \ket{100}. (Incidentally, the set of higher states \ket{k^{(0)}} includes all of the unbound states of the electron as well. If we really wanted to treat things rigorously here, we would probably write this as a sum plus an integral.)
Instead of trying to evaluate this precisely, let’s resort to some tricks to get a reasonable estimate. Since we’re considering the ground state, notice that we have an upper bound on the sum, \sum_{k>1} \frac{ |\bra{k^{(0)}} \hat{z} \ket{1^{(0)}}|^2 }{E_1^{(0)} - E_k^{(0)}} < \frac{1}{E_1^{(0)} - E_2^{(0)}} \sum_{k>1} |\bra{k^{(0)}} \hat{z} \ket{1^{(0)}}|^2 since all of the other energy differences are larger than the one we’ve factored out. But notice that this second sum is much easier to evaluate: we can write \sum_{k>1} |\bra{k^{(0)}} \hat{z} \ket{1^{(0)}}|^2 = \sum_{\textrm{all}\ k} |\bra{k^{(0)}} \hat{z} \ket{1^{(0)}}|^2 \\ = \left\langle \hat{z}^2 \right\rangle where we’ve just used completeness in the last step. Now we invoke rotational invariance: for the ground state, we must have \left\langle \hat{z}^2 \right\rangle = \left\langle \hat{x}^2 \right\rangle = \left\langle \hat{y}^2 \right\rangle = \frac{1}{3} \left\langle \hat{r}^2 \right\rangle. This last expectation value we can evaluate directly from the wavefunction: \left\langle \hat{r}^2 \right\rangle = \int_0^{\infty} dr\ r^2 |R_{10}(r)|^2 r^2 \\ = \int_0^{\infty} dr\ \frac{4r^4}{a_0^3} e^{-2r/a_0} \\ = \frac{a_0^2}{8} \int_0^{\infty} d\rho\ \rho^4 e^{-\rho}. We can evaluate this integral by integrating by parts repeatedly, taking derivatives of the polynomial term; the boundary terms at 0 and \infty always vanish, until we end up with \left\langle \hat{r}^2 \right\rangle = \frac{a_0^2}{8} \int_0^{\infty} d\rho\ (24 e^{-\rho}) = 3a_0^2, and so \left\langle \hat{z}^2 \right\rangle = a_0^2. (Note: the corresponding formula in Sakurai has a typo.) Writing E_1^{(0)} - E_2^{(0)} = \frac{3}{4}\ {\rm Ry} = -\frac{3e^2}{8a_0}, we can combine our results to find \alpha < \frac{16a_0}{3} (a_0^2) \approx 5.3 a_0^3. We can also find a lower bound on \alpha just by evaluating the first term in the sum, since all of the terms are positive-definite. However, although these bounds are relatively simple to find, it is in fact possible to carry out the complete sum; this is far too elaborate to go through here, but the exact result turns out to be \alpha = 4.5 a_0^3, which agrees well with experimental determinations.