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

X

N

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

i¼1

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Þ

2m

189

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 Þ

e2

dr dr 0

E½n ¼ Ts ½n þ T½n À Ts ½n þ V½n À

jr À r 0j

2

°° °

n°rÞn°r 0 Þ

e2 0

þ dr dr þ n°rÞVext °rÞ dr

jr À r 0 j

2

°° °

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

2

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 Þ

e2

dr dr 0 À Ts ½n:

E xc ½n FHK ½n À °5:2:5Þ

jr À r 0j

2

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Þ :

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:

n°rÞ

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

°5:2:7Þ

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.

191

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 :

0

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

and

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Þ ¼ ;

p¬¬¬¬

1 þ 1 °Þ rs þ 2 °Þrs

with

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

°5:3:1Þ

°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

that

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

2

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,

193

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:

°5:3:4Þ

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

energy.

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

°°

1

n°rÞ V°r À r 0 Þxc °r; r 0 Þ dr dr 0 :

E xc ¼

2

195

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

197

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

F

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

F

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

199

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Þ

xc

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

201

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

Hamiltonian,

H0 ¼ T þ Vint :

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

X

N

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

i

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Þ

0

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

0

vext °r; tÞ 6¼ vext °r; tÞ þ c°tÞ:

The proof of the theorem proceeds by ¬rst demonstrating that the current

0

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Þ:

i0

@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-

dependent.

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

form

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Þ þ Á Á Á

203

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Þ

F½°rÞ

lim :

°r 0 Þ

!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

205

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

!2

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

X

p

Ds ¼ dk jÈk ihÈk j;

k¼1

where the weights dk of the Slater determinants Èk are positive numbers that

P

k dk ¼ 1: This leads to a ground-state ensemble

sum to unity, so that

density

X X X

p

ns °rÞ ¼ Tr fDs n°rÞg ¼ dk hÈk jn°rÞjÈk i ¼ j™i °rÞj þ

2

fi j™i °rÞj2 ;

i:i <F i:i ¼F

k¼1

207

Problems

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 :

X

q

fj ¼ dk kj :

k¼1

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.

Problems

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

0s

208 Density functional theory

and 2 constants. Estimate the value of rs at which crystallization might

occur.

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# / " #

209

Problems

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

8

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

a30

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

X X

y

0!qs ay aqs ;

H0 ¼ E k c k ck þ qs

q;s

k

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

210

211

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

0

ei°k ÀkÞ Á yl ™ 1 þ i°k 0 À kÞ Á yl

X

À1=2 0

eiq Á l yq :

¼ 1 þ iN °k À kÞ Á

q

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

X

VÀg cy ck ;

¼N kÀg

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

becomes

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,

X X X

y y

0!q ay aq Mkk 0 °ay þ aq Þck ck 0 ;

H¼ E k c k ck þ þ °6:1:1Þ

Àq

q

0

q

k k;k

where the electron“phonon matrix element is de¬ned by

s¬¬¬¬¬¬¬¬¬¬¬¬¬

N0

jk 0 À kjVkÀk 0 ;

¼i °6:1:2Þ

Mkk 0

2M!q

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-

Àq

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.

213

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Þ

q

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

0

k;k

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

q

by unity, and so ¬nd

@E

0!°pÞ ¼

q

@hnq i

X 2hnk i°E k À E k 0 Þ

¼ 0!q þ jMkk 0 j2 :

°E k À E k 0 Þ2 À °0!q Þ2

k

215

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

X

°pÞ

2 jMkk 0 j2 hnk i°E k 0 À E k ÞÀ1 ;

0!q ¼ 0!q À °6:2:4Þ

k

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

q

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 :

2

@qx

k

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

anomaly.

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

q

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

217

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Þ

q

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

q

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

219

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

X

@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

k0

X 20!q hnk 0 i

1

¼ Ek þ jMkk 0 j2 À : °6:4:1Þ

E k À E k 0 À 0!q °E k À E k 0 Þ2 À °0!q Þ2

k0

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

2d

"

"

0vk ™ 1 À 20!D°ÞjM j d :

@k "

dE k 2 À °0!Þ2

Since

dhn°E k þ Þi

¼ À°E k þ À Þ

dE k

we ¬nd

"