. 7
( 14)


of N electrons is given by a Slater determinant of the N lowest-lying single-
particle states. Since we can much more readily solve for these, the Kohn“
Sham scheme provides us with a route to practical calculations. But there is
no free lunch “ the price we pay in the Kohn“Sham scheme is that the
equations have to be solved self-consistently. The e¬ective potential Vs will
turn out to depend on the electron density. In practical calculations, one then
typically starts by assuming an initial density. This gives an input potential
Vs , which can then be used to solve for the single-particle states. From these a
new density is obtained, which gives a new Vs . The equations are then solved
again, and this process is repeated until self-consistency is obtained, i.e., until
the input and output density in one iteration are su¬ciently close to one
another. Much e¬ort has been spent over the years to come up with e¬cient
schemes for such self-consistent calculations.
Let us then start by considering a noninteracting N-electron system in an
external potential Vs . The Hamiltonian Hs of this system is given by

Hs ¼ T þ V s :

We then apply the Hohenberg“Kohn theorem to this system. Accordingly,
there exists a unique energy functional
E s ½nŠ ¼ Ts ½nŠ þ Vs °rÞn°rÞ dr: °5:2:1Þ

We note here that Ts ½nŠ is the kinetic energy functional of a system of N
noninteracting electrons, and is consequently a di¬erent functional from the
T½nŠ that forms part of FHK ½nŠ in Eq. (5.1.6).
The ground-state density of this system is easily obtained. It is simply

ns °rÞ ¼ ji °rÞj2 ; °5:2:2Þ

where we have occupied the N single-particle states, or orbitals, that satisfy
the Schrodinger-like equation
À02 2
r þ Vs °rÞ i °rÞ ¼ E i i °rÞ; E1 E2 ÁÁÁ; °5:2:3Þ
5.2 The Kohn“Sham formulation
and have the N lowest eigenvalues E i . But we are really interested in a system
of N interacting electrons in an external potential Vext , so the question we
would like to answer is the following: can we determine the form that Vs (the
external potential of the noninteracting system) must take in order for the
noninteracting system to have the same ground-state density as the interact-
ing system in the external potential Vext ? The strategy we use is to solve for
the density using the auxiliary noninteracting system, and then insert this
density (which by construction is the same as that for the interacting system)
into an approximate expression for the total energy of the interacting system.
The ¬rst step in this process is to rewrite the energy functional E½nŠ of the
interacting system, which was given in Eq. (5.1.5), as

& '
n°rÞn°r 0 Þ
dr dr 0
E½nŠ ¼ Ts ½nŠ þ T½nŠ À Ts ½nŠ þ V½nŠ À
jr À r 0j
°° °
n°rÞn°r 0 Þ
e2 0
þ dr dr þ n°rÞVext °rÞ dr
jr À r 0 j
°° °
n°rÞn°r 0 Þ
e2 0
 Ts þ dr dr þ n°rÞVext °rÞ dr þ E xc ½nŠ: °5:2:4Þ
jr À r 0 j

Here we have added and subtracted both the kinetic energy functional Ts ½nŠ
of a noninteracting system and the direct, or Hartree, term in the electrostatic
energy. We have then de¬ned the sum of the terms in braces to be the
exchange-correlation energy functional E xc ½nŠ. From Eq. (5.1.6), this func-
tional is

n°rÞn°r 0 Þ
dr dr 0 À Ts ½nŠ:
E xc ½nŠ  FHK ½nŠ À °5:2:5Þ
jr À r 0j

We have thus swept all our ignorance about electron interactions beyond the
Hartree term under the rug that we call E xc ½nŠ. What we gain in writing E½nŠ in
this way is that we can eventually focus on developing reasonable approx-
imations for E xc ½nŠ.
According to the Hohenberg“Kohn theorem, the density n that minimizes
the functional E½nŠ is the ground-state density. Thus by taking the variation
of Eq. (5.2.4) with respect to the particle density we obtain
n°r 0 Þ
E½nŠ Ts ½nŠ
dr 0 þ Vext °rÞ þ vxc ½n°rފ ¼ 0;
¼ þ e2 °5:2:6Þ
n°rÞ n°rÞ jr À r 0j
190 Density functional theory
where we have formally de¬ned the exchange-correlation potential as

E xc ½nŠ
vxc ½n°rފ  :

We now use the auxiliary noninteracting system and its Schrodinger equa-
tion, from which we can similarly show that

Ts ½nŠ
þ Vs °rÞ ¼ 0:

By comparing this result with Eq. (5.2.6) we see that this e¬ective potential
Vs °rÞ must satisfy

n°r 0 Þ
dr 0 þ vxc °rÞ:
Vs °rÞ ¼ Vext °rÞ þ e 2
jr À r 0 j

We are now in a position to implement the self-consistent Kohn“Sham
scheme. We ¬rst choose an initial trial form of the function n°rÞ and sub-
stitute into Eq. (5.2.7) to ¬nd a trial form of Vs . We then solve Eq. (5.2.3) for
the single-particle wavefunctions i °rÞ, and use Eq. (5.2.2) to ¬nd the next
iteration for n°rÞ. When this procedure has been repeated a su¬cient number
of times that no further changes occur, then a solution for n°rÞ has been
found that not only satis¬es the Schrodinger equation for the reference non-
interacting electrons, but also is the correct density for the interacting system.
We close this section by highlighting a few points about the Kohn“Sham
formalism. First of all, it is formally exact, supposing that we can ¬nd the
exact exchange-correlation potential vxc °rÞ. Second, we have cast the solution
of the interacting N-electron problem in terms of noninteracting electrons in
an external potential Vs °rÞ. This is of great practical importance. The ground-
state wavefunction of the noninteracting system is just a Slater determinant
of the N orbitals, the so-called Kohn“Sham orbitals, with the lowest eigen-
values E . It is relatively easy to solve for these single-particle orbitals even
for as many as a few hundred electrons. The Kohn“Sham equations formally
look very much like self-consistent Hartree equations, the only di¬erence
being the presence of the exchange-correlation potential. This makes them
much simpler to solve than the Hartree“Fock equations, in which the poten-
tial is orbital-dependent. In the Kohn“Sham and Hartree equations, the
e¬ective potential is the same for every orbital.
5.3 The local density approximation
5.3 The local density approximation
Before we can actually implement the Kohn“Sham formalism, we have to
devise some workable approximation for the exchange-correlation potential
vxc °rÞ. The ¬rst such approximation to be suggested was the Local Density
Approximation, or LDA. The idea behind the LDA is very simple; it just
ignores the nonlocal aspects of the functional dependence of vxc . The true
form of vxc °rÞ will depend not only on the local density n°rÞ but also on n at
all other points r 0 , and this functional dependence is in general not known.
This di¬culty is avoided with the assumption that vxc depends only on the
local density n°rÞ, and that E xc ½nŠ can thus be written as
E xc ½nŠ ¼ xc ½nŠ n°rÞ dr;

where xc ½nŠ is the exchange-correlation energy per particle of a homogeneous
system of density n. In the LDA, we assume that the density of our (in-
homogeneous) system varies very slowly, so that the exchange-correlation
energy is locally that of a homogeneous system at the local density.
For practical calculations, we must then determine what xc ½nŠ is. Although
no general form is known exactly, the low-density and high-density limits can
be calculated analytically. Usually, the density is expressed in terms of the
dimensionless parameter rs , which is the radius of the sphere that can be
assigned to each electron (so that the volumes of all these spheres add up
to the total volume of the system), measured in units of the Bohr radius a0 .
That is,

rs ¼ °4na3 =3ÞÀ1=3 :

In the low-density limit (rs ) 1), the electrostatic potential energy dominates,
and the electrons condense into what is known as a Wigner crystal, the
energy of which can be calculated. While the density of the Wigner crystal
is not strictly uniform, we can still use the energy per electron of this system
to develop an estimate for xc ½nŠ at low densities. In the high-density limit
(rs ( 1), the kinetic energy dominates, and the random-phase approximation
becomes exact. Unfortunately, real metals have rs of the order of unity.
Usually, one then uses one of several interpolation schemes that join the
low- and high-density limits of xc . The most commonly used approximations
combine numerical calculations with Pade approximants, which are ratios of
polynomials, for the interpolations. In local spin density functional calcula-
tions, the exchange-correlation energy is frequently separated into exchange
192 Density functional theory
and correlation parts. The exchange part is in this context just the exchange
energy that we discussed in Chapter 2, and the correlation energy is every-
thing else. These two quantities are then calculated at zero and unit polariza-
tions ( ¼ 0 and  ¼ 1). The exchange energy is straightforward to calculate,
and the result is

3 9 1=3 e2
x °n;  ¼ 0Þ ¼ À ;
4 4 r s a0


x °n;  ¼ 1Þ ¼ 21=3 x °n;  ¼ 0Þ:

A popular approximation for the correlation energy c is based on Monte
Carlo calculations by Ceperley and Alder for certain values of rs for  ¼ 0
and  ¼ 1. These are then parametrized. In one example, Perdew and Zunger
use a Pade approximant in r1=2 for rs ! 1:
´ s

c °n;  ¼ 0; 1Þ ¼ ;
1 þ 1 °Þ rs þ 2 °Þrs


°0Þ ¼ À0:1423e2 =a0 ; 1 °0Þ ¼ 1:0529; 2 °0Þ ¼ 0:3334

°1Þ ¼ À0:0843e2 =a0 ; 1 °1Þ ¼ 1:3981; 2 °1Þ ¼ 0:2611:

This form is then joined smoothly to the high-density form of c (for rs 1),
which is

c °n;  ¼ 0; 1Þ ¼ A°Þ ln rs þ B°Þ þ C°Þrs ln rs þ D°Þrs : °5:3:2Þ

For  ¼ 0, the result is matched to the classical random-phase approximation
result by Gell-Mann and Brueckner, and for  ¼ 1 there exists a scaling
relation (also obtained from the random-phase approximation), which states

c °rs ;  ¼ 1Þ ¼ 1 c °rs =24=3 ;  ¼ 0Þ: °5:3:3Þ

Finally, by requiring that the correlation energy and the resulting correlation
potential, vc °n;  ¼ 0; 1Þ ¼ °d=dnÞ°nc °n;  ¼ 0; 1ÞÞ be continuous at rs ¼ 1,
5.3 The local density approximation
the parameters A°Þ, B°Þ, C°Þ, and D°Þ can be ¬xed. The result (with all
parameters given in units of e2 =a0 ) is

A°0Þ ¼ 0:0311; B°0Þ ¼ À0:048; C°0Þ ¼ 0:0020; D°0Þ ¼ À0:0166
A°1Þ ¼ 0:01555; B°1Þ ¼ À0:0269; C°1Þ ¼ 0:0007; D°1Þ ¼ À0:0048:

It then remains to interpolate exchange and correlation energies to arbi-
trary polarizations. For the exchange energy alone, one can show that there is
an exact expression for the exchange energy of a homogeneous system at
arbitrary polarizations. This expression is usually written

x °n; Þ ¼ x °n;  ¼ 0Þ þ ½x °n;  ¼ 1Þ À x °n;  ¼ 0ފ f °Þ; °5:3:5Þ

where the dimensionless function f °Þ is

°1 þ Þ4=3 þ °1 À Þ4=3 À 2
f °Þ ¼ : °5:3:6Þ
2°21=3 À 1Þ

This function has the value zero at  ¼ 0 and unity at  ¼ 1. In the LDA, one
then just replaces n and  in Eq. (5.3.5) by their local values n°rÞ and °rÞ. The
correlation energy is again a little trickier. There exists no exact interpolation
to arbitrary polarization, even for the homogeneous system, as there does for
the exchange energy. In lieu of better approximations, the ¬rst thing that
comes to mind is to use the same scaling relation for the correlation energy as
for the exchange energy, and write

c °n°rÞ; °rÞÞ ¼ c °n°rÞ; °rÞ ¼ 0Þ þ ½c °n°rÞ; °rÞ ¼ 1Þ À c °n°rÞ; °rÞ ¼ 0ފ f °°rÞÞ;

with the same scaling function f °Þ as for the exchange energy, Eq. (5.3.6).
This turns out to be surprisingly accurate: the error given by this for the
Perdew“Zunger parametrization of the exchange and correlation energies is
at most 3.5% when compared with values calculated by more laborious
numerical schemes.
The ¬rst real calculation using DFT was by Lang and Kohn, who found
the work functions of simple metals in the LDA. Before their calculation, this
had been a di¬cult problem. The Hartree approximation typically gave at
least the right sign, but wrong values, while ˜˜improvements™™ using Hartree“
Fock theory often would yield the wrong sign, meaning that according
to Hartree“Fock theory, metals were not stable! The calculations by Lang
and Kohn not only showed metals to be stable, but gave very good
194 Density functional theory
quantitative agreement with experiments. These calculations thus showed the
importance of the correlation energy. The exchange energy is typically attrac-
tive, but too large to give good results for metals without further corrections.
The inclusion of correlation energy compensates for the large exchange
It may seem surprising that the LDA gives such good results for work
functions. After all, the calculation of work functions involves metal surfaces.
While the LDA is based on the assumption of a slowly varying electron
density, the density near metal surfaces varies very rapidly, and so metal
surfaces should perhaps not be within the region of applicability of the
LDA. Other cases in point are also given by DFT“LDA calculations of
atomic systems, in which the density again varies rapidly on atomic length
scales. Yet, the LDA yields very good results (of the order of one to ten
percent error depending on the quantity). Why is that?
First of all, the LDA is by construction the exchange-correlation energy of
a physical system (an in¬nite homogeneous system). As such, it then satis¬es
many of the relations, the so-called sum rules, that are required of the
exchange-correlation energy and potential. There is also another reason
that has a deep and physical origin, and for which we need to know a little
about something called the exchange hole. Consider an electron system inter-
acting through Hartree and exchange terms alone. While the Hartree term is
blind to the spin of electrons, the exchange term is not. This term, which is
really a manifestation of the Pauli Exclusion Principle, acts to keep electrons
of like spins apart. As a consequence, if we put an up-spin electron at the
origin, there will be a de¬cit of other up-spin electrons in a neighborhood
around the origin. One can actually write down a precise expression for this
de¬cit, which, if integrated out, adds up to precisely one electron. We can
therefore think of each up-spin (or down-spin) electron as it is moving
through the system as being surrounded by a little bubble of de¬cit of up-
spin (or down-spin) electrons. This ˜˜bubble,™™ which moves with the electron,
is the exchange hole. Now, if we go beyond Hartree and Fock terms, and add
electron correlations, there is still a de¬cit of up-spin (or down-spin) charge
around the electron. The correlation e¬ects counteract the exchange term
locally to some extent, but the net de¬cit is still precisely one electron. This
net de¬cit is, quite naturally, called the exchange-correlation hole. If we
consider an electron at a point r, we can write the exchange-correlation
hole density at r 0 as xc °r; r 0 Þ, and one can rigorously show that
n°rÞ V°r À r 0 Þxc °r; r 0 Þ dr dr 0 :
E xc ¼
5.4 Electronic structure calculations
For translationally invariant interactions, xc °r; r 0 Þ ¼ xc °jr À r 0 jÞ. The
exchange-correlation energy is consequently determined by the form of this
spherically symmetric exchange-correlation hole. The LDA amounts to using
a particular approximation for the exchange-correlation hole, and this
approximation is successful in representing the form of the true exchange-
correlation hole. Part of this is due to the fortunate fact that errors from
exchange and correlation terms in the LDA cancel each other to some degree.

5.4 Electronic structure calculations
Density functional theory is the approach that is now almost universally used
in performing electronic structure calculations in condensed matter physics,
and the same is fast becoming the case in quantum chemistry. While many of
the earlier calculations used the Kohn“Sham scheme, there is now also a
large body of work, in particular for large systems, in which the approach
is to focus on minimizing the DFT“LDA expression for the ground-state
energy directly. There are in common use two di¬erent kinds of DFT
ground-state calculation. One is the so-called all-electron calculation, in
which the Coulomb potentials of the fully ionized atoms, i.e., of the bare
nuclei, are used for the external potential. All the electrons are then added to
the system and are considered in the calculations. This approach, however, is
frequently not practical for large systems consisting of hundreds of atoms, as
there are then many thousands of electrons.
This di¬culty was resolved in the previous chapter for the case of non-
interacting electrons when we introduced the concept of the pseudopotential
in Section 4.4. There we argued that the bare Coulomb potential is very
strong, and thus tightly binds the core electrons, which consequently have
little e¬ect on the properties of the system. We developed an analysis in which
the comparatively weak pseudopotential Uk given in Eq. (4.4.4) acted on a
comparatively smooth pseudo-wavefunction k . We can proceed to apply a
similar technique within DFT, but at a cost of further loss of rigor, since
there is in this case no Hohenberg“Kohn theorem that formally allows us to
cast the problem in terms of only the ground-state density. This di¬culty
stems from the fact that pseudopotentials are operators, and hence are non-
local. That is, a pseudopotential from an ion at position r acting on a valence
electron at r 0 depends on r and r 0 separately, and not just on the di¬erence
r À r 0 . A DFT can be constructed for pseudopotentials, but the basic variable
here is the so-called single-particle density“density matrix
°r; r 0 Þ, which
introduces another layer of complexity. For reasons of computational con-
venience, many users of pseudopotential methods within DFT typically
196 Density functional theory
ignore the inconsistency posed by the nonlocality of Uk . A powerful justi¬ca-
tion is that pseudopotential calculations frequently give results that are as
accurate as any all-electron calculations. This may well be due to the fact
that it is in the end necessary to resort to approximations for the exchange-
correlation energy and potential, and these approximations are probably
responsible for most of the errors in the calculations.
Density functional theory, as we have described it here, is a theory for the
ground state of an interacting electron system, although there do exist exten-
sions to thermodynamic equilibria and to excited states. While it is formally
exact and can in principle be used to determine the expectation value of any
observable of the system, practical calculations focus on the exact formula-
tion for the ground-state energy and electron density, and are constructed to
give good approximations for these quantities. It is therefore not surprising
that quantities extracted directly from the ground-state energy and electron
density tend to be more accurate than quantities found by methods for which
there exists no formal justi¬cation. But even if we restrict ourselves to only
the former class, we can, with a little ingenuity, calculate a host of quantities.
First of all, we have the total energy. By calculating this for a variety of
di¬erent possible unit cells for a solid, one can predict the crystalline struc-
ture at various macroscopic densities. Calculated values for the total energy
at di¬erent lattice parameters for each lattice structure for a given material
can then be ¬tted to an approximate equation of state, from which the bulk
modulus can be calculated. By analyzing these equations of state, predictions
can be made about phase transitions as a function of pressure. In fact, a high-
pressure phase (a so-called -tin phase) of silicon was ¬rst predicted from
total-energy DFT calculations, and later found experimentally. Most com-
mon semiconductors, such as Si, Ge, GaAs, GaP, InP, InAs, and InSb, have
been studied extensively in total-energy calculations of this sort. Unit cell
volumes and bond distances can be determined directly from the structure
having the minimum energy, with an accuracy of a few percent for the unit
cell volumes, and up to about one third of one percent (less than 0.01 A) for
bond distances. Total energies (for nonmagnetic crystals) are also well pre-
dicted, with energies accurate to within ten percent of experimental values.
The bulk modulus is found as a second derivative of the energy with respect
to volume, and is correspondingly less accurately predicted, with errors typi-
cally being some tens of percent.
Surfaces have also been the subject of many calculations. For example, the
Lang“Kohn calculations mentioned earlier were the ¬rst DFT calculations
of simple surface properties. Compared with exact calculations of surface
exchange-correlation energies for jellium systems, the LDA performs
5.4 Electronic structure calculations
acceptably well, with errors of the order of two percent. Surface reconstruc-
tions, in particular of semiconductors, have been particularly extensively
studied. In these calculations, a perfect bulk crystal is terminated at a parti-
cular plane surface, and then the atoms are allowed to move to establish new
positions of lower energy. Since the atoms at or near the surface have fewer
near neighbors than those in the bulk, they experience a net force. As these
atoms move to new equilibrium positions, the symmetry of their arrangement
at the surface is typically altered from the bulk symmetry. Perhaps the most
celebrated example is the reconstruction of Si atoms on a (111) surface, the
so-called 7 ‚ 7 reconstruction.
Phonon energies can be determined from ˜˜frozen phonon™™ calculations, in
which the lattice is given a static deformation corresponding to a phonon of a
particular branch of a chosen wavevector q and polarization s, and the total
energy is calculated and compared with the ground-state energy. These ˜˜fro-
zen phonon™™ calculations use the Born“Oppenheimer approximation, which
is based on the fact that electrons are very much lighter than the ions, and so
move much faster. The electrons can thus adjust very quickly to a change in
the ionic positions and at every instant form the ground state of the system
de¬ned by whatever positions the ions are occupying at that time. These
calculations are easiest for special phonon wavevectors q at high symmetry
points in the Brillouin zone, and give results that are typically accurate to
within a few percent.
For atoms and simple molecules, the LDA (or alternatively the Local Spin
Density Approximation, LSDA) gives good results for geometrical quanti-
ties, such as bond lengths, and for electron densities, vibrational frequencies,
and energy di¬erences such as ionization potentials. These results are often
an improvement over results obtained using the Hartree“Fock approxima-
tion. We remember also that the Hartree“Fock approximation is di¬cult and
time-consuming to use in practical calculations, since the e¬ective potential is
orbital-dependent. For open-shell atoms the LSDA tends to overestimate the
ground-state energy in comparison with the best experimental values, and
with the Hartree“Fock approximation. Another important quantity for
molecules is the bond dissociation energy, or the atomization energy. This
is the energy required to break bonds and dissociate a simple molecule into its
atomic constituents. The LDA typically does rather poorly for this quantity,
overestimating it with errors around 20%, and sometimes over 100%.
The band structures of solids are usually calculated by interpreting the
Kohn“Sham eigenvalues E k;n for Bloch states of wavevector k in band n as
being the band energies. Although there is no formal justi¬cation for this
interpretation, it usually works remarkably accurately. However, if one tries
198 Density functional theory
to calculate the band gap of insulators and semiconductors by taking the
di¬erence between the highest occupied LDA Kohn“Sham eigenvalue and
the lowest unoccupied one, the result typically underestimates the real band
gap by as much as 50%. The error has two distinct sources. One is the LDA
approximation, which introduces errors due to the inexact nature of the LDA
exchange-correlation potential. The other is of a more subtle kind. One can
rigorously show that the exchange-correlation potential must have disconti-
nuities in its functional derivative at integer particle numbers. In order to
evaluate the band gap of an N-electron system, one must include this dis-
continuity at N. The necessity of doing this adds such complexity to the
calculations that this correction is often neglected, sometimes with unfortu-
nate consequences. It is known that for some systems, such as Si, omission of
the discontinuity in the derivative of the exchange-correlation potential is
responsible for over 80% of the error in the calculated band gap.

5.5 The Generalized Gradient Approximation
As we have mentioned, the local density approximation (or its extension to
spin-polarized systems) has been the most commonly used approximation for
the exchange-correlation energy. It is simple to implement in calculations, it
gives very reasonable results, and it has the appeal of actually being the
exchange-correlation energy of something real (an in¬nite homogeneous
electron gas) and thus satis¬es sum rules and other constraints. It is,
however, tempting to regard the LDA as only the lowest-order term in an
expansion of the exchange-correlation energy in powers of the ¬rst- and
higher-order gradients of the density.
For a homogeneous electron gas, there is really only one intrinsic length
scale, and that is kÀ1 , which is proportional to nÀ1=3 . For the expansion to be
justi¬ed the length scale of variation of the density must thus be large com-
pared with kÀ1 . We can then formally expand the exchange-correlation
energy in density variations about the constant density of the homogeneous
gas by writing
E xc ½nŠ ¼ ½g00 °nÞ þ g22 °nÞ°rnÞ2 þ g42 °nÞ°r2 nÞ2

þ g43 °nÞ°r2 nÞ°rnÞ2 þ g44 °nÞ°rnÞ4 þ Á Á ÁŠ dr: °5:5:1Þ

This form satis¬es some constraints on rotational invariance, and the ¬rst
term in this expression, g00 °n°rÞ), constitutes the LDA. As terms of higher
5.5 The Generalized Gradient Approximation
order bring in more powers of rn, the coe¬cients that accompany them must
bring in an equal number of powers of °kF nÞÀ1 if the dimensions of the
expression are to be preserved. Because kF / n1=3 we conclude that
Eq. (5.5.1) is an expansion in powers of jrnj=n4=3 .
Early attempts to incorporate the corrections in Eq. (5.5.1) were focused
almost exclusively on the term g22 °n°rÞÞ, mainly because it has the attractive
property of being quadratic in rn. As a result, any variation in the energy
produces a linear response that can be calculated more or less straightfor-
wardly. However, the results obtained were discouraging in that they did not
improve on the LDA. In fact, the converse was often true, and inclusion of
gradient corrections gave worse results than use of just the LDA. Blindly
expanding the exchange-correlation energy in gradients of the density leads
to functionals that violate some sum rules that the LDA has to satisfy. For
example, such expansions often lead to an incorrect behavior for the long-
wavelength contributions to the exchange-correlation energy, and diverge or
increase sharply as k ! 0, whereas the long-wavelength contributions to the
true exchange-correlation energy vanish as k ! 0. The incorrect behavior of
this expansion is particularly easy to see for ¬nite systems, such as atoms and
molecules. For such systems, the density vanishes exponentially at distances
far from the nuclei, and this makes jrnj=n4=3 diverge.
Further analysis of the gradient expansion shows that the resulting
exchange-correlation hole is not negative de¬nite. The question then arises
of whether one can construct a gradient expansion that avoids these short-
comings, so that the resulting exchange-correlation hole satis¬es the most
important sum rules. Such analyses and ˜˜¬xes™™ have been performed, most
notably by Perdew and co-workers. The resulting gradient corrections with
exchange-correlation hole constraints restored are referred to as Generalized
Gradient Approximations (GGAs). Carefully constructed GGAs also satisfy
other physical limits, such as giving the correct exchange-correlation energy
in the homogeneous electron gas limit (rn°rÞ ! 0), coordinate scaling rela-
tions, global bounds on the exchange-correlation energy, and correctly giving
the lowest-order term g22 °n°rÞÞ in the expansion Eq. (5.5.1). The exchange-
correlation energy in the GGA is thus written as
E GGA ¼ f °n" ; n# ; rn" ; rn# Þ dr; °5:5:2Þ

where the function f is a universal function of the spin-up and spin-down
densities and their gradients. There exist several versions of the GGA, all of
which are parametrized somewhat di¬erently.
200 Density functional theory
The great strength of the GGA lies in the dramatic improvement it gives
over the LDA in calculating such properties as bond dissociation energies,
which the LDA may overestimate by as much as 100%, while the GGA gives
errors typically of the order of ten percent or less. The GGA also gives a great
improvement over the LDA for bulk moduli of solids, with an error of
around ten percent, compared with around 20% for the LDA (or LSDA).
With modern GGA exchange-correlation potentials, atomic and molecular
quantities can now be calculated with chemical accuracy. Before the GGA,
this was only possible using so-called con¬guration-interaction schemes, in
which the true many-body wavefunction is expanded in some small set of
Slater determinants composed of suitable atomic basis functions. This is an
extremely arduous and computationally intensive approach. Density func-
tional theory with the GGA, on the other hand, is essentially no more com-
plicated than Hartree or LDA calculations. The exchange-correlation
potential is only slightly more complicated than in the LDA, and (more
importantly) it is still a multiplicative potential, leading to simple, e¬ective
single-particle equations. The combined advantages of simple equations and
high accuracy have resulted in a revolution in quantum chemistry.
We close this section with a cautionary note. While mean errors, averaged
over many compounds, atoms, or molecules, tend to be improved in the
GGA compared with the LDA, this is no guarantee that the GGA will be
better for some speci¬c calculation. For example, while the LDA on average
gives poorer results for the bulk modulus, it has been found better for GaAs
than some GGA calculations. New forms for the functions used in Eq. (5.5.2)
are continually being suggested, each of which has its own advantages. The
so-called PBE (Perdew“Burke“Ernzerhof) GGA is one of those constructed
to satisfy the constraints and limits mentioned earlier. There are other GGAs
not constructed this way, but obtained by brute-force ¬tting to a large data
set. While such GGAs may give excellent results for quantities included in the
¬tting procedure, they can yield large errors for other quantities. Because
these kinds of ¬tting procedure tend not to satisfy physical constraints and
limits and often violate, for example, the uniform electron gas limit, their use
is limited to specialized applications.

5.6 More acronyms: TDDFT, CDFT, and EDFT
We ¬nish our discussion of density functional theory with a brief tour of some
later developments. We start with time-dependent density functional theory
(TDDFT), which extends the reach of stationary density functional theory in
a very powerful way. Not only are strongly time-dependent phenomena
5.6 More acronyms: TDDFT, CDFT, and EDFT
accessible to computations, but TDDFT also provides a natural way to
calculate the excitation energies of a system.
We have previously seen how the Hohenberg“Kohn theorem establishes a
one-to-one correspondence between external potentials and electron densities
for time-independent systems. It is natural to ask under what circumstances
something similar holds true for systems in which the external potential
depends on time, so that Vext ¼ Vext °r; tÞ. The answer is given by the
Runge“Gross theorem, which we will here only state and not prove. This
theorem establishes a one-to-one correspondence between density n°r; tÞ and
its governing time-dependent external potential v°r; tÞ, but there is a catch:
the correspondence can be established only for a speci¬ed initial many-body
state É0 , and consequently the functional relationships between density and
potential implicitly depend on É0 . However, this will not cause any problem
if É0 is a nondegenerate ground state. To state the theorem, we begin with a
system initially in an eigenstate É0 of the usual homogeneous N-electron

H0 ¼ T þ Vint :

At some time t0 we turn on a time-dependent, scalar potential

Vext °tÞ ¼ vext °ri ; tÞ:

Here the only requirements we put on vext °r; tÞ are that it be ¬nite and Taylor-
expandable about t0 . We do not require É0 to be an eigenstate of Vext °t0 Þ. The
theorem then states that the densities n°r; tÞ and n 0 °r; tÞ that evolve from
a common initial state É0 under the in¬‚uence of two potentials vext °r; tÞ
and vext °r; tÞ are di¬erent provided that the potentials di¬er by more than
a purely time-dependent but spatially uniform function, so that
vext °r; tÞ 6¼ vext °r; tÞ þ c°tÞ:

The proof of the theorem proceeds by ¬rst demonstrating that the current
densities resulting from the application of vext °r; tÞ and vext °r; tÞ necessarily
must di¬er, given the property of Taylor-expandability of the potentials.
Once this is established, one can then use the continuity equation to show
that the densities n°r; tÞ and n 0 °r; tÞ must also di¬er. Therefore, the density
determines the potential, and so the Hamiltonian, and thus the expectation
value of any operator, is uniquely determined by the time-dependent density.
202 Density functional theory
However, we must bear in mind that all these relationships implicitly depend
on the choice of initial state É0 .
Once the Runge“Gross theorem is established, one can then continue and
derive Kohn“Sham equations, just as for the time-independent case, by con-
sidering a noninteracting system with the same density as the interacting one.
Formally, the equations look just as one would expect for the time-dependent
Kohn“Sham orbitals ™°r; tÞ:

@™j °r; tÞ 02 2
¼À r þ vs ½nŠ°r; tÞ ™j °r; tÞ:
@t 2m

Here, the e¬ective potential vs °r; tÞ is given by

vs °r; tÞ ¼ vxc °r; tÞ þ vext °r; tÞ þ vH °r; tÞ;

with the exchange-correlation and Hartree potentials now being time-
Once the time-dependent formalism is established, one can derive formally
exact expressions for the density“density response functions, which
describe how the density of the system changes in response to an external
potential that couples to the density. As an example we consider the
linear response function obtained by expanding the density response to
¬rst order in the applied potential. We suppose an external potential of the

vext °r; tÞ ¼ v0 °rÞ þ v1 °r; tÞ°t À t0 Þ;

where °xÞ is the Heaviside step function, which vanishes for x < 0 and is
unity for x ! 0. This form of external potential is allowed by the Runge“
Gross theorem. For t < t0 , the external potential is v0 °rÞ, which we take to be
the potential of the ion cores of the system, and the system is in its ground
state corresponding to v0 °rÞ with density n0 °rÞ. At a time t0 an additional,
time-dependent perturbation v1 °r; tÞ is applied. According to the Runge“
Gross theorem, the time-dependent density n°r; tÞ is then a functional of
the external potential (since we had ¬xed the state of the system prior to
turning on the time-dependent perturbation). We can then expand the density
in a Taylor series about n0 °rÞ:

n°r; tÞ ¼ n0 °rÞ þ n1 °r; tÞ þ Á Á Á
5.6 More acronyms: TDDFT, CDFT, and EDFT
The ¬rst-order response of the density is by de¬nition linear in v1 °r; tÞ, and
thus can be written
°r; t; r 0 ; t 0 Þv1 °r 0 ; t 0 Þ dr 0 dt:
n1 °r; tÞ ¼

The response function °r; t; r 0 ; t 0 Þ is the amount by which the density n°r; tÞ
varies when a change in external potential is applied that is localized in both
time and space, and of the form °r À r 0 Þ°t À t 0 Þ. It can be expressed as a
functional derivative, but of a type more complicated than those encountered
in Section 5.2, and written as

n½vext Š°r; tÞ
°r; t; r 0 ; t 0 ; n0 Þ  ;
vext °r 0 ; t 0 Þ

where we have made the dependence of  on the initial density n0 explicit.
The formal de¬nition of the functional derivative F½°rފ=°r 0 Þ of the scalar
functional F½°rފ with respect to the scalar function °r 0 Þ is a function of r
and r 0 , and is de¬ned as

F½°rÞ þ °r À r 0 ފ À F½°rފ
 lim :
°r 0 Þ 

We note that the dimensions of functional derivatives are not what they at
¬rst appear to be, since  has the dimensions of  multiplied by a volume. The
functional derivative n½vext Š°r; tÞ= vext °r 0 ; t 0 Þ thus has dimensions of (number
density)/(energy ‚ volume ‚ time).
Similarly, we can consider the response of the noninteracting auxiliary
system to the external potential vext °r; tÞ, which yields its response function

n½vs Š°r; tÞ
s °r; t; r 0 ; t 0 Þ  : °5:6:1Þ
vext °r 0 ; t 0 Þ

It is possible to show that when this response function is Fourier transformed
with respect to the time di¬erence °t À t 0 Þ and is given in terms of the static,
unperturbed Kohn“Sham orbitals ™i °rÞ one ¬nds

X ™i °rÞ™*°rÞ™*°r 0 Þ™j °r 0 Þ
0 j i
s °r; t; r ; !Þ ¼ ° fj À fi Þ : °5:6:2Þ
0! À °i À j Þ þ i
i; j
204 Density functional theory
Here the i are the Kohn“Sham eigenvalues,  is a positive in¬nitesimal to ¬x
the poles of s , and fi are the occupation numbers (typically 0 or 1) of the
Kohn“Sham orbitals. The sum in Eq. (5.6.2) vanishes unless one of the
orbitals is occupied and one is empty. This can be understood in terms of
electrons making virtual transitions from occupied to unoccupied states. This
means that any practical calculations of the response function must include
enough of the unoccupied orbitals to capture the principal contributions
to s .
What is now left is to relate the response function  of the interacting
system to that of the noninteracting system. First, we note that, by assump-
tion, the densities of the interacting and the noninteracting systems are iden-
tical. Therefore, we can think of n°r; tÞ in Eq. (5.6.1) as the density of the
noninteracting system. This density is a functional of the e¬ective potential
vs °r; tÞ of the noninteracting system, and vs is in its turn a functional of the
external potential vext °r; tÞ. Putting that all together using the chain rule for
functional di¬erentiation, we obtain
n°r; tÞ vs °r 00 ; t 00 Þ 00 00
0 0
°r; t; r ; t Þ ¼ dr dt : °5:6:3Þ
vs °r 00 ; t 00 Þ vext °r 0 ; t 0 Þ

If we then write the exchange-correlation potential as

vxc ½nŠ°r; tÞ ¼ vs ½nŠ°r; tÞ À vext °r; tÞ À vH ½nŠ°r; tÞ

we can use this relation between vs , vext , vH , and vxc to rewrite the second
functional derivative in Eq. (5.6.3) as
°t À t 00 Þ vxc °r; tÞ n°r 00 t 00 Þ
vs °r; tÞ
¼ °r À r 0 Þ°t À t 0 Þ þ dr 00 dt 00 :
vext °r 0 ; t 0 Þ jr À r 00 j n°r 00 ; t 00 Þ vext °r 0 t 0 Þ

We then insert this into Eq. (5.6.1), and use Eq. (5.6.2) to arrive ¬nally at
°t 00 À t 000 Þ vxc °r 00 ; t 00 Þ
°r; t; r 0 ; t 0 Þ ¼ s °r; t; r 0 ; t 0 Þ þ 00 00
s °r; t; r ; t Þ þ
jr 00 À r 000 j n°r 000 ; t 000 Þ
‚ °r 000 ; t000 ; r 0 ; t 0 Þ dr 00 dt 00 dr 000 dt 000 : °5:6:4Þ

Equation (5.6.4) is a formally exact equation for the density“density response
function of the interacting system. Note that this expression (being a
linear response quantity) depends only on the ground-state properties of
the system, and can be calculated if the ground-state Kohn“Sham orbitals
and their eigenvalues have been calculated. By multiplying Eq. (5.6.4) by the
5.6 More acronyms: TDDFT, CDFT, and EDFT
perturbing potential under consideration, and integrating over r 0 and t 0 , we
can also derive a formally exact expression for the linear density response to
the time-dependent potential of our choosing. Moreover, after a little
thought we realize that the exact density“density response function °r; r 0 ; !Þ
has poles at the exact excitation energies of the interacting system. Thus, the
time-dependent linear response theory gives us a practical way to calculate
the excitation energies of an interacting system. This calculation would be exact
if we were fortunate enough to have an exact expression for the quantity

vxc °r; tÞ
fxc °r; t; r 0 ; t 0 Þ  :
n°r 0 ; t 0 Þ

This quantity is called the exchange-correlation kernel. Perhaps the most
natural, and also the most frequently used, approximation for it is one
made in the spirit of the LDA, and is called the adiabatic local density
approximation (ALDA). Here, we consider a homogeneous electron gas
with a time-dependent perturbation that varies su¬ciently slowly that the
system is ˜˜adiabatic,™™ meaning that the exchange-correlation energy per
particle at r at any instant t is that of a homogeneous electron gas in an
external potential vext °r; t 0 Þ°t 0 ¼tÞ . This implies that there is no frequency
dependence of fxc in the ALDA. Just as for the LDA, it is found that the
ALDA works reasonably well, especially for low-frequency perturbations,
and, just as for the LDA, attempts to improve on it turn out to be di¬cult.
The reason for this is very deep: any simple attempt to include a frequency
dependence in a local approximation will violate something called the
Harmonic Potential Theorem, which is related to Galilean invariance. The
conclusion is, somewhat discouragingly, that there is no local, frequency-
dependent exchange-correlation kernel. Instead, to build a theory that does
not run into trouble with Galilean invariance, one has to consider the current
density as the fundamental variable, not the particle density. Such a theory is
called current density functional theory (CDFT).
Current density functional theory was developed to deal with magnetic
systems for which the external magnetic ¬eld is strong enough that the full
canonical kinetic energy
1 e
T¼ pÀ A
2m c

has to be considered. As a consequence, there will in general be diamagnetic
currents in response to the external magnetic ¬eld. The Hohenberg“Kohn
206 Density functional theory
theorem remains perfectly valid, even in the presence of a magnetic ¬eld, so
we can formally still extract everything we want from ˜˜standard™™ density
functional theory. However, experience has shown that practical calculations
are usually more accurate if the broken symmetry can be explicitly included
in the formalism. This is, for example, the case with spin density functional
theory. While we could in principle use regular density-only functional the-
ory, the approximations we make when we consider the spin densities sepa-
rately lead to more accurate calculations. It was thought, therefore, that
building the current density into the formulation would likewise lead to
more accurate computational approaches. While this essentially turned out
not to be the case, current density functional theory has given us a way to
formulate more re¬ned approximations for the exchange-correlation kernel
than can be done using time-dependent density functional theory, and which
are not encumbered by ¬‚agrant violations of Galilean invariance.
Finally, in this discourse on more acronyms, we brie¬‚y discuss ensemble
density functional theory (EDFT), which helps us deal with degeneracies.
The origin of EDFT is simple enough to relate. Consider a time-independent
system, for which we are solving the Kohn“Sham equations. It does not
matter whether we are using the LDA or GGA, or even some exact
exchange-correlation energy. According to the Kohn“Sham scheme, we ¬ll
all the orbitals having the lowest energy eigenvalues. The occupation
numbers fi for the N lowest-lying orbitals are then unity while fi ¼ 0 for
all the others. Suppose now that, while iterating our equations to reach self-
consistency, we ¬nd that as we ¬ll the orbitals in order of ascending energy, p
eigenvalues at the highest occupied level F are degenerate. Then what are we
supposed to do? The answer is formally given by ensemble density functional
theory. Here, the density is not constructed from a single Slater determinant
of Kohn“Sham orbitals, but by an ensemble of Slater determinants Èk . We
do this by writing the density matrix operator as

Ds ¼ dk jÈk ihÈk j;

where the weights dk of the Slater determinants Èk are positive numbers that
k dk ¼ 1: This leads to a ground-state ensemble
sum to unity, so that

ns °rÞ ¼ Tr fDs n°rÞg ¼ dk hÈk jn°rÞjÈk i ¼ j™i °rÞj þ
fi j™i °rÞj2 ;
i:i <F i:i ¼F
where now the occupation numbers fi of the orbitals on the ˜˜Fermi surface™™
are fractional, 0 fi < 1. Formally, these fractional occupation numbers are
given as follows. Consider a particular orbital ™j . This orbital does or does
not appear in a particular Slater determinant Èk , and we can assign an
occupation number kj ¼ 0; 1 to this orbital in that Slater determinant. The
occupation number fj for this orbital is then given by a weighted average of
the kj , with the weights given by the ensemble weights dk :

fj ¼ dk kj :

One can show that EDFT has very general validity, and in fact formally
avoids some problems that occasionally crop up for ˜˜regular™™ density func-
tional theory. These problems refer to what is known as v-representability,
and are concerned with under what circumstances one can use the auxiliary
noninteracting system to represent the density of the interacting one.
The most profound example of applications of EDFT comes in the frac-
tional quantum Hall e¬ect, about which we shall learn more in Chapter 10.
These are strongly correlated systems (that is, the physics is determined
entirely by the interactions), and there is no systematic perturbation theory
that applies. Furthermore, it is known that the ground-state energy has cusps
at the densities at which the fractional quantum Hall e¬ect occurs. This
implies that the LDA exchange-correlation potential has discontinuities,
and so applications of DFT to fractional quantum Hall systems ¬rst have
to deal with these discontinuities. If we assume that this can be done, and
then na±¨ vely apply the Kohn“Sham scheme we ¬nd that in general all Kohn“
Sham eigenvalues are degenerate. This is apparently a situation where we
have to take EDFT very seriously. A practical calculation scheme has been
devised for these systems, and applications show that EDFT in the LDA
gives very good results, and provides perhaps the only practical scheme with
which to study inhomogeneous fractional quantum Hall systems.

5.1 A Wigner crystal forms when the potential energy gained by separating
the electrons outweighs the cost of the extra kinetic energy needed to
localize them. It is not easy to calculate the change in electrostatic energy
per electron, but on dimensional grounds it must be of the form
À 1 e2 =a0 rs , while the localization energy must be 2 02 =ma2 r2 , with 1
208 Density functional theory
and 2 constants. Estimate the value of rs at which crystallization might

5.2 Explicit calculation of even the kinetic energy functional Ts ½nŠ of a non-
interacting electron gas is very di¬cult. Consider the e¬ect of a weak
external potential V cos gx for which g ) kF on n°rÞ and on the kinetic
energy, and hence evaluate the function h22 in the expansion of Ts ½nŠ in
powers of rn that is analogous to the g22 °nÞ de¬ned in Eq. (5.5.1).

5.3 Now consider the calculation of the kinetic energy functional Ts ½nŠ of a
noninteracting electron gas in the opposite limit from Problem 5.2.
Consider the e¬ect of a weak external potential V cos gx for which g (
kF as causing a periodic variation in the radius of the Fermi surface. The
electron density n°rÞ is again perturbed from its original uniform value
n0 °rÞ, and this time the kinetic energy functional can be written as
Ts ½nŠ ¼ Ts ½n0 Š 1 þ °hn2 i À n2 Þ=n2 :
0 0

Evaluate .

5.4 In a collinear spin-polarized electron gas, the total electron density n can
be separated into components n ¼ n" þ n# , where n" and n# are the
densities of electrons with hsz i ¼ þ 0=2 and hsz i ¼ À 0=2, respectively.
In this case, the Hohenberg“Kohn theorem can be formulated in terms
of n" and n# , or, equivalently, in terms of n and the polarization
 ¼ °n" À n# Þ=n. In particular, the exchange energy can then be con-
sidered a functional E x ½n" ; n# Š of n" and n# . Show that this functional
must satisfy the relation

E x ½n" ; n# Š ¼ 1 E x ½n" ; n" Š þ 1 E x ½n# ; n# Š;
2 2

where E x ½n ; n Š means that the functional is evaluated with both argu-
ments set equal to n . You may use the fact that the exchange interaction
acts only between electrons with the same value of hsz i. Since the
exchange energy in the local density approximation is a simple function
of electron density, we then have in this approximation in d dimensions
d d r n þ n :
E x ½n" ; n# Š / " #
Find the value of by using dimensional analysis. Then, by changing
variables to n and  from n" and n# , ¬nd the form of the interpolating
function f °Þ for a two-dimensional electron gas [cf. Eq. (5.3.6)].

5.5 An old interpolation formula, but one that is simply expressed in closed
form, for the correlation energy of a uniform unpolarized electron gas of
N electrons is due to Nozieres and Pines, who wrote

E c =N ¼ 0:0311 ln rs À 0:115:

The unit of energy here is the hartree, which is equal to about 27.2 eV.
Use this form of the correlation energy to estimate the ground-state
energy of the helium atom. Use a simple singlet Ansatz wavefunction
for the two electrons, the spatial part of which is

É°r1 ; r2 Þ ¼ exp ½À2°r1 þ r2 Þ=a0 Š;

where r1 and r2 are the coordinates of the two electrons and a0 is the
Bohr radius. First, evaluate the expectation value of the kinetic and
potential energies. To evaluate the Hartree and correlation energy, you
need to derive an expression for the electron density using the Ansatz
wavefunction. [You do not need to consider the exchange energy, since
there is only one electron of each spin orientation, and it cannot
exchange with itself!] Compare your result to the experimental value of
À98.975 eV. How large is the correlation correction? Then calculate the
electron interaction energy in ¬rst-order perturbation theory by forming
the quantity hÉjV°r1 À r2 ÞjÉi. Discuss the di¬erence you ¬nd from the
previous result.

5.6 In Problem 2.3 you estimated the density at which the Hartree“Fock
approximation predicts that the uniform electron gas has a paramagnetic
to ferromagnetic transition. This critical density will change when the
correlation energy is included. Estimate this revised critical density by
using the scaling relation (5.3.3) and the Nozieres“Pines expression for
the correlation energy for an unpolarized electron gas given in Problem
5.5. You will need to estimate the ground-state energy (kinetic, exchange,
and correlation) of both the unpolarized and of a fully polarized ( ¼ 1)
electron gas at a ¬xed density, and ¬nd when they become equal.
Chapter 6
Electron“phonon interactions

6.1 The Frohlich Hamiltonian
In this chapter we shall consider some of the consequences of the interaction
of phonons with electrons, and in particular with the electrons in a simple
metal. The subject is a complicated and di¬cult one, in that we need to call
on most of the knowledge that we have of the behavior of the electron gas
and of lattice vibrations. A complete calculation should really start with the
Hamiltonian of a lattice of bare ions whose mutual interaction would include
the long-range Coulomb potential. One would then add the electron gas,
which would shield the potential due to the ions in the manner indicated in
Section 1.7. It is, however, possible to explore many of the consequences of
the electron“phonon interaction by use of a simpler model. In this model we
take for granted the concept of screening, and assume that the ions interact
with each other and with the electrons only through a short-range screened
potential, and we treat the electrons themselves as independent fermions. For
a Bravais lattice our unperturbed Hamiltonian is then simply
0!qs ay aqs ;
H0 ¼ E k c k ck þ qs

the phonon frequencies !qs being proportional to q as q ! 0. To this we add
the interaction, H1 , of the electrons with the screened ions. We assume that at
any point the potential due to a particular ion depends only on the distance
from the center of the ion “ an assumption sometimes known as the rigid-ion
approximation “ so that in second-quantized notation
X y
hkjV°r À l À yl Þjk 0 ic k ck 0
H1 ¼
k;k 0 ;l
X 0 y
ei°k ÀkÞ Á °lþyl Þ VkÀk 0 ck ck 0 :
k;k 0 ;l

6.1 The Fro¨hlich Hamiltonian
Here V°rÞ is the potential due to a single ion at the origin, and VkÀk 0 its
Fourier transform. With the assumption that the displacement yl of the ion
whose equilibrium position is l is su¬ciently small that °k 0 À kÞ Á yl ( 1 we
can write
ei°k ÀkÞ Á yl ™ 1 þ i°k 0 À kÞ Á yl
À1=2 0
eiq Á l yq :
¼ 1 þ iN °k À kÞ Á

Then H1 can be split into two parts,

H1 ¼ HBloch þ HeÀp ;

the ¬rst term, HBloch , being independent of the lattice displacements. We have
X 0 y
ei°k ÀkÞ Á l VkÀk 0 ck ck 0
HBloch ¼
k;k 0 ;l
VÀg cy ck ;
¼N kÀg

where the g are reciprocal lattice vectors, and
X 0 y
HeÀp ¼ iN À1=2 ei°k ÀkþqÞ Á l °k 0 À kÞ Á yq VkÀk 0 ck ck 0
k;k 0 ;l;q
X y
°k 0 À kÞ Á ykÀk 0 VkÀk 0 ck ck 0 :
¼ iN 1=2
k;k 0

In terms of the annihilation and creation operators de¬ned in Eq. (3.6.4) this

X 1=2
N0 y y
°k 0 À kÞ Á sVkÀk 0 °a k 0 Àk;s þ akÀk 0 ;s Þc k ck 0 ;
HeÀp ¼ i
2M!kÀk 0 ;s
k;k 0 ;s

where the summation now also includes the three polarization vectors, s, of
the phonons. (Because s and Às represent the same phonon mode, care is
necessary in using this expression. With a di¬erent convention for the direc-
tion of s, the plus sign in this expression would become a minus sign.) For
simplicity we shall assume the phonon spectrum to be isotropic, so that the
phonons will be either longitudinally or transversely polarized. Only the
longitudinal modes, for which s is parallel to k 0 À k, then enter HeÀp “ a
212 Electron“phonon interactions
fact in accord with the semiclassical viewpoint of Section 1.7. We shall also
neglect the e¬ects of HBloch , the periodic potential of the stationary lattice.
With these simpli¬cations we are left with the Fro¨hlich Hamiltonian,
y y
0!q ay aq Mkk 0 °ay þ aq Þck ck 0 ;
H¼ E k c k ck þ þ °6:1:1Þ
k k;k

where the electron“phonon matrix element is de¬ned by
jk 0 À kjVkÀk 0 ;
¼i °6:1:2Þ
Mkk 0

with the phonon wavenumber q equal to k À k 0 , reduced to the ¬rst Brillouin
zone if necessary.
The interaction HeÀp can be considered as being composed of two parts “
y y
terms involving ay c k ck 0 and terms involving aq ck ck 0 . These may be repre-
sented by the diagrams shown in Figs. 6.1.1(a) and 6.1.1(b), respectively. In
the ¬rst diagram an electron is scattered from k 0 to k with the emission of a
phonon of wavenumber k 0 À k. The total wavenumber is then conserved, as
is always the case in a periodic system, unless the vector k 0 À k lies outside
the ¬rst Brillouin zone, so that q ¼ k 0 À k þ g for some nonzero g. Such
electron“phonon Umklapp processes do not conserve wavenumber, and
are important in contributing to the electrical resistivity of metals.
In the remainder of this chapter we shall study some of the consequences of
the electron“phonon interaction for the equilibrium properties of materials.
We shall see that not only are the electron and phonon excitation spectra
modi¬ed, but that there can also appear an e¬ective attractive interaction
between electrons caused by the electron“phonon interaction.

Figure 6.1.1. The Frohlich Hamiltonian includes an interaction term in which an
electron is scattered from k 0 to k with either emission (a) or absorption (b) of a
phonon. In each case the total wavenumber is conserved.
6.2 Phonon frequencies and the Kohn anomaly

6.2 Phonon frequencies and the Kohn anomaly
The e¬ect of the electron“phonon interaction on the phonon spectrum may
be seen by using perturbation theory to calculate the total energy of the
system described by the Frohlich Hamiltonian (6.1.1) to second order in
HeÀp . We have

E ¼ E 0 þ hÈjHeÀp jÈi þ hÈjHeÀp °E 0 À H0 ÞÀ1 HeÀp jÈi;

with E 0 the unperturbed energy of the state È having nq phonons in the
longitudinally polarized mode q and nk electrons in state k. The ¬rst-order
term vanishes from this expression, since the components of HeÀp act on È
either to destroy or create one phonon, and the resulting wavefunction must
be orthogonal to È. In second order there is a set of nonvanishing terms, as
the phonon destroyed by the ¬rst factor of HeÀp to act on È can be replaced
by the second factor of HeÀp , and vice versa. We then ¬nd the contribution,
E 2 , of the second-order terms to be

X y
Mkk 0 °ay þ aq Þck ck 0 °E 0 À H0 ÞÀ1
E 2 ¼ hÈj Àq
k;k 0
X y
Mk 00 k 000 °ay 0 þ aq 0 Þck 00 ck 000 jÈi
‚ Àq
k 00 ;k 000
X y y
jMkk 0 j2 ½ay ck ck 0 °E 0 À H0 ÞÀ1 aÀq ck 0 ck
¼ hÈj Àq
k;k 0

y y
þ aq ck ck 0 °E 0 À H0 ÞÀ1 ay c k 0 ck ŠjÈi; °6:2:1Þ

all other terms having zero matrix element. The ¬rst term in the brackets in
(6.2.1) can be represented as in Fig. 6.2.1(a). An electron is ¬rst scattered
from k to k 0 with the absorption of a phonon of wavenumber Àq ¼ k 0 À k.
The factor °E 0 À H0 ÞÀ1 then measures the amount of time the electron is
allowed by the Uncertainty Principle to stay in the intermediate state k 0 .
In this case the energy di¬erence between the initial and intermediate
states is E k þ 0!Àq À E k 0 , and so a factor of °E k þ 0!Àq À E k 0 ÞÀ1 is con-
tributed. The electron is then scattered back into its original state with
the re-emission of the phonon. We can represent the second term in
Eq. (6.2.1) by Fig. 6.2.1(b), and there ¬nd an energy denominator of
E k À 0!q À E k 0 . A rearrangement of the a™s and c™s into the form of number
214 Electron“phonon interactions

Figure 6.2.1. These two processes contribute to the energy of the electron“phonon
system in second-order perturbation theory.

operators then gives
X hnÀq i hnq þ 1i
E2 ¼ jMkk 0 j2 hnk °1 À nk 0 Þi þ : °6:2:2Þ
E k À E k 0 þ 0!Àq E k À E k 0 À 0!q
k;k 0

Here hnk i and hnk 0 i are electron occupation numbers while hnq i and hnÀq i
refer to phonon states. It may be assumed that !q ¼ !Àq , and hence that in
equilibrium hnq i ¼ hnÀq i. One may then rearrange Eq. (6.2.2), to ¬nd
" #
X 2°E k À E k 0 Þhnq i 1 À hnk 0i
E ¼ E0 þ jMkk 0 j2 hnk i þ ; °6:2:3Þ
°E k À E k 0 Þ2 À °0!q Þ2 E k À E k 0 À 0!q

the term in hnk nk 0 nq i cancelling by symmetry.
The e¬ect of the electron“phonon interaction on the phonon spectrum is
contained in the term proportional to hnq i in Eq. (6.2.3). We again identify
the perturbed phonon energy, 0!°pÞ , with the energy required to increase hnq i
by unity, and so ¬nd

0!°pÞ ¼
@hnq i
X 2hnk i°E k À E k 0 Þ
¼ 0!q þ jMkk 0 j2 :
°E k À E k 0 Þ2 À °0!q Þ2
6.2 Phonon frequencies and the Kohn anomaly

Figure 6.2.2. This alternative way of considering the process of Fig. 6.2.1(a) suggests
that a phonon spends part of its time as a virtual electron“hole pair.

If we neglect the phonon energy in the denominator in comparison with the
electron energies we have
2 jMkk 0 j2 hnk i°E k 0 À E k ÞÀ1 ;
0!q ¼ 0!q À °6:2:4Þ

where, as before, k 0 ¼ k À q. One may picture the origin of this change in
phonon frequency by redrawing Fig. 6.2.1(a) in the form of Fig. 6.2.2, in
which the ¬rst interaction is represented, not as the scattering of an electron,
but as the creation of an electron“hole pair. One can then say that it is the
fact that the phonon spends part of its time in the form of an electron“hole
pair that modi¬es its energy.
One interesting consequence of Eq. (6.2.4) occurs in metals when q has a
value close to the diameter, 2kF , of the Fermi surface. Let us suppose q to be
in the x-direction and of magnitude 2kF , and evaluate 0@!°pÞ =@qx . If we
neglect the variation of Mkk 0 with q the electron“phonon interaction contri-
butes an amount

X @E kÀq
jMkk 0 j2 hnk i°E kÀq À E k ÞÀ2 :

On substituting for E kÀq one ¬nds the summation to contain the factors
hnk i°kx À kF ÞÀ2 . These cause a logarithmic divergence when the summation
is performed, and thus indicate that the phonon spectrum has the form
indicated qualitatively in Fig. 6.2.3. The kink in the spectrum when q ¼ 2kF
re¬‚ects the in¬nite group velocity of the phonons at that point, and consti-
tutes the Kohn anomaly. Its importance lies in the fact that even for very
complex metals there should always be such an image of the Fermi surface in
the phonon spectrum. One should then in principle be able to gain informa-
tion about the electronic structure of metals by studies of the phonon spec-
trum alone. This is important in that neutron di¬raction experiments can be
performed to determine !q at high temperatures and in impure samples
216 Electron“phonon interactions

Figure 6.2.3. The kink that appears in the phonon dispersion curve when the phonon
wavenumber q is equal to the diameter of the Fermi surface is known as the Kohn

where many other techniques are not useful. Such e¬ects have been seen most
clearly in lead, in which the electron“phonon interaction is very strong.

6.3 The Peierls transition
When we calculated the change in phonon energy caused by the electron“
phonon interaction, the answer we found contained a sum over terms of the
form °E kÀq À E k ÞÀ1 . Even though these could become in¬nite for certain
values of k, the sum in Eq. (6.2.4) fortunately remained ¬nite. It is only
the derivative @!°pÞ =@q that diverges, giving rise to the Kohn anomaly. We
are saved from ¬nding an in¬nite perturbation to the phonon frequency by
the fact that we perform the sum in three-dimensional k-space. The integra-
tion proceeds over both the magnitude of k and its direction, and it is the
directional integration that rescues us from the embarrassment of ¬nding an
in¬nite negative perturbation to the phonon frequency.
However, before congratulating ourselves on this good fortune we might
ponder the fact that there are a number of physical systems that approximate
one-dimensional conductors. These include various organic materials, and in
particular certain families of polymers. A polymer is a long chain molecule of
covalently connected units (monomers), and a polymeric crystal can be
thought of as an assembly of parallel chains. Most everyday polymers, like
6.3 The Peierls transition
the polyethylene shown in Fig. 6.3.1(a), for example, have completely ¬lled
electronic bands, and are thus electrical insulators.
Polyacetylene, however, di¬ers from polyethylene in having only one
hydrogen atom on each carbon of the chain backbone as in Fig. 6.3.1(b).
The double bond between every second pair of carbon atoms is the chemist™s
representation of the extra pair of electrons not used in bonding a hydrogen
atom. There is then e¬ectively only one conduction electron per unit cell if we
think of the polymer as a one-dimensional lattice of spacing a.
The band structure for this Bravais lattice has the qualitative form shown
in Fig. 6.3.2(a). The wavenumber at the Fermi level (in one dimension the
Fermi surface has become a pair of Fermi points!) is Æ=2a, as the band is

Figure 6.3.1. The chemical structure of polyethylene (a) makes it an insulator, but
polyacetylene (b) can be a conductor through the motion of the defects (c) known as
charged solitons.

Figure 6.3.2. The half-¬lled band (a) of polyacetylene would make it a metallic
conductor, but dimerization opens up an energy gap.
218 Electron“phonon interactions
half-¬lled. This causes the phonon wavenumber q that connects the two
Fermi points to be =a, which is just half a reciprocal lattice vector. This
phonon, of wavelength 2a, is thus the mode of highest unperturbed fre-
quency, and vibrates in such a way that adjacent monomers move in anti-
phase with each other. To estimate the perturbation to this frequency, we put
q ¼ °=aÞ°1 þ Þ and perform the sum in the one-dimensional version of
Eq. (6.2.4) to ¬nd the result

0!°pÞ À 0!q $ ln : °6:3:1Þ

As expected, this tends to À1 as  tends to zero.
The resolution of this di¬culty comes when we realize that, on its way to
negative in¬nity, 0!°pÞ must pass through zero, and a zero frequency phonon
is simply a static distortion of the lattice. The lattice has thus dimerized, so
that now the e¬ective unit-cell dimension is 2a. The Fermi ˜˜points™™ and the
¬rst Brillouin zone thus coincide. Bragg scattering of the electron states at the
Fermi points thus causes an electronic band gap to open at k ¼ Æ=2a, as in
Fig. 6.3.2(b). If we now repeat the sum in Eq. (6.2.4), but with the new form
for E k , the in¬nity is removed. This process of dimerization is known as the
Peierls transition.
We could have reached a similar conclusion without the use of many-body
theory by examining the total electronic energy as a function of the dimer-
ization distortion. Displacing every second monomer by a distance u opens
up a band gap 2V that is proportional to u. As only the lower band is
occupied, the electronic energy of a k-state is lowered by an amount equal
to E °newÞ À E °oldÞ . From Eq. (4.3.2) we can evaluate the sum of these energies.
k k
The result is that the total electronic energy is lowered by an amount
$ V 2 ln V 2 , which is always negative for small V, and hence for small dimer-
ization parameter u. Adding an elastic energy of distortion, which will be
proportional to u2 , and hence to V 2 , does not change this fact, and so the
Peierls transition is always favored at low temperatures.
Because the energy of the distorted system is an even function of u, the
ground state is degenerate. This corresponds to the chemical picture in which
the double bond can form either between monomers 2n and 2n þ 1, or
between monomers 2n and 2n À 1. If the symmetry of the unperturbed
state is broken di¬erently in two di¬erent parts of the chain, then a defect
occurs at the boundary of these regions, as shown in Fig. 6.3.1(c). This defect
is mobile, since it just requires the displacement of one monomer to move the
defect a distance 2a along the chain, and may also be charged, as it can
represent a region with added electron concentration. The name soliton is
6.4 Polarons and mass enhancement
used to describe this type of excitation. It bears some qualitative similarity to
the solitary waves that we encountered in Section 1.3 when we looked at wave
propagation in the Toda chain.

6.4 Polarons and mass enhancement
Just as the energy of the phonons in a crystal is altered by interaction with the
electrons, so also does the converse process occur. We examine Eq. (6.2.3) in
the limit of low temperatures, when hnq i vanishes for all q. The perturbed
energy of an electron “ once again the energy needed to ¬ll an initially empty
unperturbed state “ is given by
@E 1 À hnk 0 i hnk 0 i
¼ Ek þ jMkk 0 j2 À
@hnk i E k À E k 0 À 0!q E k 0 À E k À 0!q
X 20!q hnk 0 i
¼ Ek þ jMkk 0 j2 À : °6:4:1Þ
E k À E k 0 À 0!q °E k À E k 0 Þ2 À °0!q Þ2

The ¬rst term in the brackets is independent of nk 0 , and is thus a correction to
the electron energy that would be present for a single electron in an insulating
crystal. Indeed, in an ionic crystal the e¬ect of this term may be so great as to
change markedly the e¬ective mass of an electron at the bottom of the con-
duction band. It then becomes reasonable to use the term polaron to describe
the composite particle shown in Fig. 6.2.1(b) that is the electron with its
attendant cloud of virtual phonons. The name arises because one considers
the positive ions to be attracted towards the electron, and thus to polarize the
lattice. If this polarization is too great then second-order perturbation theory
is inadequate, and di¬erent methods must be used.
The second term in the brackets in Eq. (6.4.1) expresses the dependence of
the electron energy on the occupancy of the other k-states. In a metal it has
the e¬ect of causing a kink in the graph of energy against wavenumber for
the perturbed electron states, as shown in Fig. 6.4.1. This kink occurs at the
Fermi wavenumber kF , and leads to a change in the group velocity vk of the
electron. We ¬nd an expression for 0vk by di¬erentiating Eq. (6.4.1) with
respect to k. Instead of the free-electron expression we ¬nd

dX 20!q hnk 0 i
@E k
0vk ¼ 1À jMkk 0 j2 :
@k °E k À E k 0 Þ2 À °0!q Þ2
dE k k 0
220 Electron“phonon interactions

Figure 6.4.1. The electron“phonon interaction changes the e¬ective electron energy
in such a way that the velocity is lowered in the vicinity of the Fermi surface. This
gives rise to an increase in the observed electronic speci¬c heat.

We then argue that the major contribution to the derivative of the summa-
tion comes from the rapid variation of hnk 0 i at the Fermi surface. We change
the summation over k 0 to an integral over E k 0 and make the approximations
" "
of replacing Mkk 0 and !q by their average values M and !, and D°E k 0 Þ by its
value at the Fermi energy . With E k 0 À E k written as  we then have
@E k hn°E k þ Þi
0vk ™ 1 À 20!D°ÞjM j d :
@k "
dE k 2 À °0!Þ2


dhn°E k þ Þi
¼ À°E k þ  À Þ
dE k

we ¬nd


. 7
( 14)