Pick'n mix

The Schrödinger equation for realistic systems quickly becomes unwieldy, and analytical solutions are only available for very simple systems - the ones we have described as fundamental systems in this module. Numerical approaches can cope with more complex problems, but are still (and will remain for a good while) limited by the available computer power. Approximations are necessary to cope with real systems.

Within limits, we can use a pick and mix approach, i.e. use linear combinations of solutions of the fundamental systems to build up something akin to the real system. Strictly speaking, linear combinations need to be solutions of the same differential equation to provide an accurate result. This means e.g. that an LCAO approach mixing states from two different types of atom may be a rather crude approximation. There are two mathematical techniques, perturbation and variation theory, which can provide a good approximation along with an estimate of its accuracy. These are described below.

Approximative approaches we have seen so far include:

Perturbation theory vs. variation principle

The Schrödinger equation, \hat{H}\psi=E\psi, gives us two handles to refine a problem to make it more realistic: the Hamiltonian and the wave function. The energy eigenvalues are just scalar values that respond to changes we make to the other terms. The two approaches are compared below.

In both cases (and more generally, too), the energy eigenvalues are found using E=\int_{-\infty}^{\infty}\psi^*\hat{H}\psi{\rm d}x. This prescribes a method of calculation which involves three steps:

  1. Let the Hamiltonian operator act on the wave function: \hat{H}\psi.
  2. Multiply the result with the complex conjugate of the wave function: \psi^*\hat{H}\psi.
    Note that the operator does not act on ψ*. It is important to write it to the left of the operator; otherwise any operations such as differentiation contained in the operator would modify ψ* as well as ψ, which it shouldn't.
  3. Integrate over all space: E=\int_{-\infty}^{\infty}\psi^*\hat{H}\psi{\rm d}x.
    Replace dx with whatever coordinates you are using, e.g. r,θ,φ in the case of spherical coordinates.

The recipe must be followed in this particular order as operators and their operands in general do not commute, i.e. the result is different depending on the order the terms are applied.

Perturbation theory

\hat{H}=\hat{H}_0+\lambda\hat{H}_1

The unperturbed Hamiltonian of a known system is modified by adding a perturbation with a variable control parameter λ, which governs the extent to which the system is perturbed.

The perturbation can affect the potential, the kinetic energy part of the Hamiltonian, or both. As an example, consider a double well potential created by superimposing a periodic potential on a parabolic one. This might apply e.g. to a defect in a crystalline lattice.

Fig.: Double well potential.

The Schrödinger equation for the perturbed system is
\hat{H}\psi_n=E_n\psi_n,
that for the unperturbed (known) system is
\hat{H}_0\psi_n^0=E_n^0\psi_n^0.
The index n just serves to identify a particular wave function (e.g. one which uses the same quantum number both for the perturbed and the unperturbed variant).

Variation principle

\psi=c_1\psi_1+c_2\psi_2

Two (or more) wave functions are mixed by linear combination. The coefficients c1, c2 determine the weight each of them is given. The optimum coefficients are found by searching for minima in the potential landscape spanned by c1 and c2.

Fig.: Potential landscape.

The energy minima are found by finding the differentials \frac{\partial E}{\partial c_1}=\frac{\partial E}{\partial c_2}=0 and setting them to zero.

Given that
E=\int\psi^*\hat{H}\psi{\rm d}x=(c_1\psi_1^*+c_2\psi_2^*)\hat{H}(c_1\psi_1+c_2\psi_2),
we see that the energy eigenvalue has separate contributions coming from ψ1 or ψ2 only:
=c_1^2\int\psi_1^*\hat{H}\psi_1{\rm d}x+c_2^2\int\psi_2^*\hat{H}\psi_2{\rm d}x\,+\cdots
and a cross term known as the overlap integral:
\cdots+\,c_1c_2(\int\psi_1^*\hat{H}\psi_2{\rm d}x+\int\psi_2^*\hat{H}\psi_1{\rm d}x).

The overlap integral causes the difference in energy between bonding and anti-bonding molecular states.

Derivation of the energy correction in a perturbed system

We'll leave the fine detail of the variation technique to the fourth-year module, but will derive here a recipe (for the impatient: it's highlighted at the bottom!) by which we can determine the energy correction due to a perturbation acting on a known system (i.e. one whose Hamiltonian, wave functions and eigenvalues we know already).

The Schrödinger equation of the perturbed system contains the perturbed Hamiltonian (known) and the perturbed wave functions and eigenvalues (as yet unknown):
\hat{H}\psi_n=(\hat{H}_0+\lambda\hat{H}_1)\psi_n=E_n\psi_n.

Since both wave functions and eigenvalues are unknown, they are expanded into Taylor series as an approximation (leaving out the index n for clarity here):
\psi=\psi|_{\lambda=0}+\lambda\left.\frac{\partial\psi}{\partial\lambda}\right|_{\lambda=0}+\frac{\lambda^2}{2!}\left.\frac{\partial^2\psi}{\partial\lambda^2}\right|_{\lambda=0}+\frac{\lambda^3}{3!}\left.\frac{\partial^3\psi}{\partial\lambda^3}\right|_{\lambda=0}+\cdots, and
E=E|_{\lambda=0}+\lambda\left.\frac{\partial E}{\partial\lambda}\right|_{\lambda=0}+\frac{\lambda^2}{2!}\left.\frac{\partial^2E}{\partial\lambda^2}\right|_{\lambda=0}+\frac{\lambda^3}{3!}\left.\frac{\partial^3E}{\partial\lambda^3}\right|_{\lambda=0}+\cdots.
Here, each term is a progressively smaller correction (i.e. the series converges to the true value of ψ or E, respectively). To achieve this, they are weighted with prefactors in progressive powers of λ and progressive inverse factorials -- the prefactors are diminishing very rapidly given that the control parameter λ ranges from zero to one. The differentials take into account how ψ and E respond to changes in λ, i.e. the severity of the perturbation. The notation |λ=0 indicates that all differentials are evaluated in the limit of very small λ.

We can write the two series as sums:
\psi=\sum_{i=0}^{\infty}\frac{\lambda^i}{i!}\left.\frac{\partial^i\psi}{\partial\lambda^i}\right|_{\lambda=0} and E=\sum_{i=0}^{\infty}\frac{\lambda^i}{i!}\left.\frac{\partial^iE}{\partial\lambda^i}\right|_{\lambda=0}.
For ease of notation, we can bundle up the differentials and factorials into abbreviated symbols:
\psi^{(i)}:=\frac{1}{i!}\left.\frac{\partial^i\psi}{\partial\lambda^i}\right|_{\lambda=0} and E^{(i)}:=\frac{1}{i!}\left.\frac{\partial^iE}{\partial\lambda^i}\right|_{\lambda=0},

Then ψ(0) and E(0) are the unperturbed wave functions and eigenvalues, while ψ(i) and E(i) are the changes to the wave functions and energy eigenvalues due to the perturbation, evaluated to the i-th order of perturbation theory. Using this, the two sums can be written as
\psi=\sum_{i=0}^{\infty}\lambda^i\psi^{(i)} and E=\sum_{i=0}^{\infty}\lambda^iE^{(i)}.

With the series expansions, the Schrödinger equation \hat{H}\psi=E\psi becomes (\hat{H}_0+\lambda\hat{H}_1)(\psi^{(0)}+\lambda\psi^{(1)}+\lambda^2\psi^{(2)}+\cdots)=(E^{(0)}+\lambda E^{(1)}+\lambda^2E^{(2)}+\cdots)(\psi^{(0)}+\lambda\psi^{(1)}+\lambda^2\psi^{(2)}+\cdots)

By factoring out, we can split this into terms of different order in λ:
\hat{H}_0\psi^{(0)}+\lambda(\hat{H}_1\psi^{(0)}+\hat{H}_0\psi^{(1)})+\lambda^2(\hat{H}_1\psi^{(1)}+\hat{H}_0\psi^{(2)})+\cdots=E^{(0)}\psi^{(0)}+\lambda(E^{(0)}\psi^{(1)}+E^{(1)}\psi^{(0)})+\lambda^2(E^{(0)}\psi^{(2)}+E^{(1)}\psi^{(1)}+E^{(2)}\psi^{(0)})+\cdots
The zeroth-order term corresponds to the unperturbed system, and we can use the first-order term to derive the energy corrections, E(0). The equation can be split into separate equations for each order, which can be solved independently.

Concentrating on the 1st-order equation, \hat{H}_1\psi^{(0)}+\hat{H}_0\psi^{(1)}=E^{(0)}\psi^{(1)}+E^{(1)}\psi^{(0)}, we see that the LHS is the sum of the perturbation applied to the original wave function and the original Hamiltonian applied to the (unknown) 1st-order correction to the wave function.

Generally, as explained at the top of this page, we can find energy eigenvalues by sandwiching the Hamiltonian between the wave function and its complex conjugate and integrating over all space: E=\int_{-\infty}^{\infty}\psi^*\hat{H}\psi{\rm d}x. We can bring the LHS terms in that shape by multiplying from the left with ψm(0)* and integrating:
\int\psi_m^{(0)*}\hat{H}_1\psi_n^{(0)}{\rm d}x+\int\psi_m^{(0)*}\hat{H}_0\psi_n^{(1)}{\rm d}x=\int\psi_m^{(0)*}E_n^{(0)}\psi_n^{(1)}{\rm d}x+\int\psi_m^{(0)*}E_n^{(1)}\psi_n^{(0)}{\rm d}x.
Here the index n reappears, because the wave function we multiply the equation with needn't be the same as the one that's already there -- they could be ones with different values of a quantum number. To distinguish them, we use m and n as indices.

On the RHS, the energies are just scalars and can be taken outside the integrals:
\int\psi_m^{(0)*}\hat{H}_1\psi_n^{(0)}{\rm d}x+\int\psi_m^{(0)*}\hat{H}_0\psi_n^{(1)}{\rm d}x=E_n^{(0)}\int\psi_m^{(0)*}\psi_n^{(1)}{\rm d}x+E_n^{(1)}\int\psi_m^{(0)*}\psi_n^{(0)}{\rm d}x.

Hermitian operators

Since the energy eigenvalue must be a real number rather than a complex one, the result of E=\int_{-\infty}^{\infty}\psi^*\hat{H}\psi{\rm d}x must be real. This is the case if the imaginary parts of ψ and just cancel out, i.e. \int\psi^*\hat{H}\psi{\rm d}x\stackrel{!}{=}\int\psi(\hat{H}\psi)^*{\rm d}x. Any operator that meets this criterion is described as an Hermitian operator, after mathematician Charles Hermite. All Hamiltonians in quantum mechanics are Hermitian, but the mathematical concept is not limited to quantum mechanics.


It is also useful to consider that xy*=(x*y)*, as we can see by factorising the two complex numbers x=a+ib and y=c+id:
xy*=(a+ib)(c-id)=ac-iad+ibc+bd=(ac+bd)+i(bc-ad)
and
(x*y)*={(a-ib)(c+id)}*=(ac+bd)+i(bc-ad)
Thus xy*=(x*y)*.

The second term on the left needs some further attention. Because the Hamilton operator is Hermitian (see box), we can swap the two wave functions:
\int\psi_m^{(0)*}\hat{H}_0\psi_n^{(1)}{\rm d}x=\int\psi_n^{(1)}(\hat{H}_0\psi_m^{(0)})^*{\rm d}x.
Using xy*=(x*y)* (see box), we have =\int(\psi_n^{(1)*}\hat{H}_0\psi_m^{(0)})^*{\rm d}x,
and the integral doesn't interfere with the complex conjugate:
=[\int\psi_n^{(1)*}\hat{H}_0\psi_m^{(0)}{\rm d}x]^*.
Now, H0ψm(0) is just the LHS of the unperturbed Schrödinger equation for state m, and we can replace it with the corresponding RHS:
=[\int\psi_n^{(1)*}E_m^{(0)}\psi_m^{(0)}{\rm d}x]^*,
and since the eigenvalue is constant and real, it can be taken out of the integral and lose its complex-conjugate asterisk:
=E_m^{(0)}[\int\psi_n^{(1)*}\psi_m^{(0)}{\rm d}x]^*.
Finally, we just need to undo the double complex conjugate:
=E_m^{(0)}\int\psi_m^{(0)*}\psi_n^{(1)}{\rm d}x.

With the second term of the perturbed Schrödinger equation now simplified, we have:
\int\psi_m^{(0)*}\hat{H}_1\psi_n^{(0)}{\rm d}x+E_m^{(0)}\int\psi_m^{(0)*}\psi_n^{(1)}{\rm d}x=E_n^{(0)}\int\psi_m^{(0)*}\psi_n^{(1)}{\rm d}x+E_n^{(1)}\int\psi_m^{(0)*}\psi_n^{(0)}{\rm d}x.
Note that the second and third integral are the same, so we can combine the two terms on the LHS and put the other two on the right:
(E_m^{(0)}-E_n^{(0)})\int\psi_m^{(0)*}\psi_n^{(1)}{\rm d}x=E_n^{(1)}\int\psi_m^{(0)*}\psi_n^{(0)}{\rm d}x-\int\psi_m^{(0)*}\hat{H}_1\psi_n^{(0)}{\rm d}x.

It is useful to consider what the knowns and unknowns are in this equation. All symbols that have an index (0) are known, because they relate to the original, unperturbed system. H1 is also known, as we've started by defining it as a perturbation on the original Hamiltonian. The only unknowns are ψn(1) and En(1), the corrections to the wave function and the energy eigenvalue, respectively. We can work them out, separately, by considering the two cases where m=n (eigenvalue) and mn (wave function).

For m=n, the LHS is zero because the two energies are the same. On the RHS, the integral \int\psi_n^{(0)*}\psi_n^{(0)}{\rm d}x=1 because wave functions are normalised, so integrating one over all space always gives 1. That leaves:
0=E_n^{(1)}-\int\psi_n^{(0)*}\hat{H}_1\psi_n^{(0)}{\rm d}x. Solve for the unknown:
E_n^{(1)}=\int\psi_n^{(0)*}\hat{H}_1\psi_n^{(0)}{\rm d}x,
in other words, to find the energy correction E(1) in a perturbed system, apply the perturbation H1 to the unperturbed wave function ψ(0) in the same way as you would normally determine the energy eigenvalue.

Note that we do not need to know or work out the perturbed wave function in order to calculate the energy correction! Also, the control parameter λ was necessary to separate the terms of different order, but it has dropped out of the equation a long way up - it doesn't matter how strong the perturbation is.

The mn case can be used to work out the correction to the wave function if required.

An example: the anharmonic oscillator

The harmonic oscillator is based on Hooke's law and the idea of perpetual conversion of kinetic into potential energy and vice versa. In practice, this description is only an adequate approximation of e.g. a chemical bond in the limit of very small oscillations. The anharmonic oscillator model addresses that by adding higher-order terms of the elongation x to account for deviations farther away from x=0. These extra terms are a perturbation of the Hamiltonian for the harmonic oscillator:
\hat{H}=\hat{H}_0+\hat{H}_1=-\frac{\hbar^2}{2m}\frac{\rm d^2}{{\rm d}x^2}+\frac{1}{2}kx^2+cx^3+dx^4,
where the red parts represent the perturbation.

Knowing the wave function of the harmonic oscillator and its eigenvalues, we can calculate the anharmonic energy corrections by applying the perturbation Hamiltonian, H1 to the unpeerturbed wave function of the harmonic oscillator. This is one of the examples included in worksheet 1.

Being able to combine atomic states to build molecular ones, and to add complications by means of perturbation theory, we can move on to make more complex aggregates of atoms: solids.