Fundamentals of Condensed Matter Physics (PDFDrive) (104-160)
Fundamentals of Condensed Matter Physics (PDFDrive) (104-160)
Fundamentals of Condensed Matter Physics (PDFDrive) (104-160)
M2
states D(ω)
Density of
Infinite
slope
ω
ω2
Figure 4.17 Density of states in three dimensions near an M2 van Hove singularity.
M3
states D(ω)
Density of
Infinite
slope
ω
ω3
Figure 4.18 Density of states in three dimensions near an M3 van Hove singularity.
Status D(ω)
Density of
ω
ω0 ω1 ω2 ω3
Figure 4.19 The total density of states obtained by including one of each of the critical points. This can be viewed as the simplest
D(ω) for N = 3 if the M1 and M2 singularities are each taken to be threefold degenerate.
For the Debye spectrum with a minimum in the acoustic branch consistent with
ω(q) = vs q, (4.126)
where the velocity of sound vs can be different in different directions. Hence the rules
above do not apply in this case. Therefore, for the low-frequency acoustic modes, our
earlier result D(ω) ∝ ω2 applies.
89 4.5 Critical points and van Hove singularities
LA
TA2
TA1
q
qBZ
Figure 4.20 The phonon dispersion curve for Al in the (1,1,0) direction (TA1 and TA2 are the transverse acoustic modes, and
LA is the longitudinal acoustic mode).
D(ω)
TA1
ω
TA2
ω
LA
ω
LA
Figure 4.21 Sketch of the density of states for the three different phonon branches of Al and the total density of states.
90 4 Lattice vibrations and phonons
For a real crystal example, we can consider Al, which is an fcc crystal with one atom
per unit cell. Therefore, there are three phonon branches, as shown schematically in
Fig. 4.20. The corresponding D(ω) appears in Fig. 4.21, with van Hove singularities easily
recognizable in the spectrum.
Part I Problems
I.1. Crystal structure of MgB2 . MgB2 is a superconductor with a high transition tem-
perature of 39K and has multiple values for the superconducting energy gap. The
material has a rather simple structure as given by the model below in Fig. I.1 – the
boron atoms (dark balls) form graphitic honeycomb layers (stacked identically on top
of each other) and the magnesium atoms (light balls) are located above the center of
the B hexagons exactly in between the layers.
Let a be the B–B distance in the layer and c be the distance between the boron
layers.
(a) Construct the primitive lattice vectors, the reciprocal lattice, the Wigner–Seitz
cell, and the Brillouin zone for MgB2 .
(b) What is the volume of the primitive unit cell? What is the volume of the Bril-
louin zone? Give the locations of the basis atoms in terms of the primitive lattice
vectors.
(c) In general, crystal structures have many symmetry operations; those operations
can be one of four types: rotation, inversion, reflection, or translation (to find
out more about them, read Tinkham’s Group Theory and Quantum Mechanics).
Find, for MgB2 , all symmetry operations which leave the crystalline Hamiltonian
invariant. (You may neglect spin–orbit interaction and other relativistic effects.)
(d) Draw the irreducible part of the Brillouin zone.
I.2. The GaN Crystal. The 2014 Nobel Prize in Physics was given for work on the GaN
crystal, specifically for using it in blue light-emitting diodes (LEDs). A model of its
crystal structure is given in Fig. I.2, where the Ga atoms are denoted by the light balls
and the nitrogen atoms by the dark balls. It has a hexagonal lattice, with a1 = a2 = a,
a3 = c, and the angle between the vectors a1 and a2 is 120◦ .
91
92 Part I Problems
a3
a2
a1
(a) Write down the explicit expressions of a set of primitive reciprocal lattice vectors
b1 , b2 , and b3 in Cartesian coordinates (in terms of the values a and c). Sketch
the first Brillouin zone. What is the volume of the first Brillouin zone in terms
of a and c?
(b) How many atoms of each kind are there in the basis of this crystal structure?
(c) In addition to translations, there are many symmetry operations under which
the above structure remains invariant. Name any six specific operations among
them.
(d) GaN is used in blue LEDs as it is a wide-bandgap semiconductor. Why is it not
surprising that this material is an insulator? Explain your answer in terms of the
band theory.
I.3. Born–Oppenheimer approximation for a molecule. Consider a diatomic molecule
with linear dimension a and reduced mass M. Show that the electronic, vibrational,
and rotational energy levels can be estimated as terms which are proportional to suc-
cessively higher orders in the small ratio m/M (m is the electron mass). How is this
related to the Born–Oppenheimer approximation discussed in the text?
I.4. Hartree–Fock approximation. The Hartree–Fock approximation is a simple yet
important model for understanding electron–electron interaction in crystals.
(a) Derive the Hartree–Fock self-consistent field equations from the variational prin-
ciple using a Slater determinant of single-particle orbitals as a trial many-electron
wavefunction.
(b) Consider the total electronic energy of a system (e.g. a molecule) in the Hartree–
Fock approximation. Calculate the change in the energy of the system if an
electron is promoted from an occupied ith orbital to an unoccupied jth orbital,
assuming that the orbitals are unchanged after the excitation. How is this excita-
tion energy related to the eigenvalues of the Hartree–Fock equations? This kind
of excitation is called a neutral excitation, which occurs for example in a photo-
excitation process, since no electrons are added or removed from the system. The
result is known as Koopmans’ theorem.
I.5. Born–von Karman boundary condition. The use of the Born–von Karman bound-
ary condition makes possible a simple and effective description of infinitely extended
93 Part I Problems
crystals. With this boundary condition, the number of k-states is finite and equiv-
alent in the successively higher Brillouin zone. This leads to different choices of
representation for the electron band structure.
(a) Describe and justify the three possible representations of the band structure
(that is, the electron energy vs. wavevector) for a crystal: reduced zone scheme,
repeated zone scheme, and extended zone scheme.
(b) Show that there are exactly the same number of allowed wavevectors in the first
Brillouin zone as there are unit cells in the crystal.
I.6. Energy bands of elemental solids. Name three elements in each category which
are:
(a) semiconductors (and give their crystal structure);
(b) semimetals (and give their crystal structure);
(c) metals in the fcc structure;
(d) metals in the bcc structure;
(e) metals in the hcp structure.
Consider only elements lighter than Rn. Distinguish metals which are transition
metals in (c), (d), and (e).
I.7. Point group. Prove the general theorem that only rotations that are multiples of 60◦ ,
90◦ , 120◦ , and 180◦ can be symmetry elements of a crystal. (Hint: consider rotation
by θ . By writing it as a matrix, show that 1 + 2 cos θ = integer.)
I.8. Lattice sums. It is often necessary to do “lattice sums” and sums over “q-space.” For
lattice vectors Rn , find the results for the following sums:
iq·R
(a) e n,
n
iq·R
2
(b)
e n
,
iq·R
n
(c) e n .
q
I.9. Free electron gas. Treat a metal as a free-electron gas and calculate:
(a) the density of states as a function of energy,
(b) the density of states at the Fermi level in terms of n and EF ,
(c) the electronic heat capacity.
I.10. Heat capacity from electrons in a semimetal. Find an expression for the electronic
contribution to the heat capacity of a semimetal at very low temperatures for electrons
of mass me and holes of mass mh . Let Ee be the energy of the bottom of the electron
band and let Eh be the energy of the top of the hole band.
I.11. Another form of the Kronig–Penney model. Consider an electron in one dimension
in the presence of the potential
U(x) = U0 (x − ma)(ma + b − x), (I.1)
m
where (x) is a step function (= 0 if x < 0, = 1 if x > 0), a is the lattice constant, b
is the width of the potential well, and U0 is the potential height/depth.
(a) Focus on a single unit cell and write down the boundary condition on the
Schrödinger equation that leads to Bloch states in this unit cell.
94 Part I Problems
(b) Solve the Schrödinger equation in this cell by taking summation of planewaves
and imposing suitable boundary conditions at the positions 0, b, and a. The result
is a condition on the Bloch index k.
2 a−2
(c) Take the limit b → 0, U0 → ∞, and U0 b → W0 a h̄ m . The condition on the
Bloch index should become
W0
cos ka = sin Ka + cos Ka, (I.2)
Ka
where K = 2mε/h̄2 , with ε the energy of the state.
(d) Produce plots of the two lowest-energy bands following from Eq. (I.2) for a = 1,
m = 1, and W0 = 1. Display the bands in the reduced zone scheme and the
extended zone scheme.
I.12. Structure factor. An energy gap is not necessarily opened up at the Brillouin zone
(BZ) edges, due to the vanishing Fourier components of the crystal potential. One
example is nearly-free electrons in a hexagonal lattice. The empirical pseudopotential
method introduced in the text is closely related to the Fourier components of the
crystal field.
(a) Define the term “structure factor” and evaluate it for an hcp lattice. Show that
its value leads to a vanishing energy gap for states on the hexagonal face of the
Brillouin zone for this structure in the nearly-free electron model (NFEM).
(b) Show explicitly that the first three pseudopotential form factors needed for a
band calculation for Si are V(G2 ) for G2 = 3, 8, 11 (in units of (2π/a)2 ). What
are the values of the G2 ’s for the next five relevant form factors?
I.13. Tight-binding model. Using the tight-binding method, and assuming one spherically
symmetric atomic orbital/atom and only nearest-neighbor interactions,
(a) find E(k) for the simple cubic, bcc, and fcc crystal structures,
(b) find E(k) for a two-dimensional (2D) square lattice,
(c) plot E(k) = constant in the first BZ in part (b).
In (a), (b), and (c), treat the hopping integral as a parameter and the overlap integrals
equal to zero for neighboring orbitals.
I.14. Diatomic chain. Calculate the band structure of a diatomic chain with the
Hamiltonian
†
H= 1 c†2i c2i + 2 c†2i+1 c2i+1 + (−t) c2i c2i+1 + h.c. , (I.3)
i i
where 1 > 2 , and c†i and ci are the creation and destruction operators of elec-
tronic orbitals at the ith site, respectively. Suppose t (1 − 2 ) and expand (k) in
t/(1 − 2 ).
I.15. Tight-binding calculation in 2D. In high Tc cuprate superconductors, there are lay-
ers of Cu–O planes, as shown in Figure I.3. Consider the above 2D square Cu–O
crystal with the distance between copper atoms given by a, and the oxygen atoms
95 Part I Problems
Cu O Cu O Cu
O O O
Cu O Cu O Cu
O O O
Cu O Cu O Cu
located midway between the copper atoms. Using the tight-binding method and con-
sidering only nearest-neighbor interactions (assuming orthogonal orbitals), do the
following calculations.
(a) Find E(k), assuming that there is one d-orbital of (x2 − y2 ) symmetry on each
copper atom and one p-orbital on each oxygen atom that is pointing along the
nearest-neighbor direction. How many electronic bands are there for this model
Hamiltonian?
(b) Sketch the energy bands along the [10] and [11] directions.
I.16. Schrödinger equation of crystal electron. Find the Schrödinger equation for the
periodic part of the Bloch function (i.e. uk (r)).
I.17. Dirac equation and spin–orbit interaction. Derive the spin–orbit contribution to
the one-electron Hamiltonian starting with the Dirac equation. How does this affect
Bloch’s theorem? Discuss translational invariance, rotational invariance, inversion
symmetry, and time-reversal invariance for this Hamiltonian.
I.18. Orthogonality of Bloch functions. Show that two Bloch functions eik·r uk (r) and
eik ·r uk (r) are orthogonal if k − k is not a reciprocal lattice vector.
I.19. Parity of Bloch waves. Consider a crystal potential with inversion symmetry.
(a) Show that for most wavevectors in the first Brillouin zone, the Bloch functions
cannot have a definite parity in spite of the crystal inversion symmetry.
(b) There are specific wavevectors for which the Bloch functions can have definite
parity. Find the wavefunctions and describe their locations in the first Brillouin
zone.
I.20. f -sum rule. Examine the problem of deriving the f-sum rule for the case of degener-
ate bands. Assume E(k0 ) is s-fold degenerate. Show that for a crystal with inversion
symmetry, the first-order contributions of perturbation theory are zero and one must
consider second-order effects (Van Vleck perturbation theory). Discuss the analytic
properties of the bands for the s = 2 case and define an effective mass tensor (if
possible) for this case.
I.21. Delta function potential. Solve the Schrödinger equation for a single one-
dimensional delta function potential. Under what conditions does a bound state
96 Part I Problems
exist? What is the energy of the bound state? Do the same calculation for two one-
dimensional delta function potentials separated by a distance a. Describe (without
repeating the previous Kronig–Penney calculations) what happens as one goes to the
many delta function problems.
I.22. High symmetry points. Look up the names (given in Greek letters) of the high
symmetry points in the first Brillouin zone for the common lattices given below.
Compute and plot E(k) for the free-electron model (FEM) in the following symmetry
directions in the first Brillouin zone:
(a) simple cubic, to R;
(b) bcc, to N;
(c) fcc, to X;
(d) hcp, to A.
I.23. Deep core levels. Estimate the effect of other atoms on the deep core levels of an
atom in a crystal in the following way. Consider a model simple cubic crystal with
a = 5Å and one atom per unit cell of atomic number z = 6. Calculate roughly the
width of the “energy band” formed by the lowest-lying s electrons, and calculate
numerically the ratio between this bandwidth and the energy of the atomic s states.
Make appropriate approximations in evaluating any integrals that may appear in your
answer.
I.24. Buckyballs. Buckyballs are molecules with 60 carbon atoms. They crystallize in an
fcc lattice. (Assume no band overlap.)
(a) How many filled electron bands are there?
(b) How many partially filled electron bands do you expect?
When this system is doped with potassium, one can form K3 C60 in an fcc lattice with
the three potassium atoms inserted in the primitive cell.
(c) How many filled electron bands are there for K3 C60 ?
(d) How many partially filled electron bands do you expect in K3 C60 ?
I.25. Density of states.
(a) Show explicitly that the density of states in three dimensions has the form
ds
D(ω) = . (I.4)
(2π ) 3 |∇ q ω|
(b) For a 2D lattice, derive the expression for the behavior of the density of states
near critical points giving:
(i) a minimum,
(ii) a maximum,
(iii) a saddle point.
√
I.26. Velocity of sound. The classical formula for the velocity of sound is vs = R/ρ,
where R is the bulk modulus and ρ is the atomic mass density.
(a) Evaluate vs for a free electron gas using the electron bulk modulus and density.
(b) Are the sound waves the same as ordinary sound waves?
(c) Evaluate this vs for a gas of He-3 atoms, neutrons, muons, positrons, and protons,
assuming in each case a concentration of n fermions/cm3 .
(d) What is kF for each group of fermions in (c)?
97 Part I Problems
α β α
m M m M
a b
I.27. Vibrational modes in graphene. Consider a graphene crystal and let us ignore the
degree of freedom along the direction that is perpendicular to the carbon atomic
plane, i.e. viewing this system as a strictly 2D system.
(a) How many acoustic phonon branches and optical phonon branches are there for
this 2D crystal? Why?
(b) Using a Debye model for the phonon density of states, obtain the temperature
dependence of its phonon-heat capacity at low temperature.
(c) What is the minimum number of van Hove singularities for each optical phonon
branch of such a 2D crystal? Sketch the functional form of an M0 van Hove
singularity in the density of states for an optical branch of this system.
I.28. Debye model. Consider a three-dimensional crystal with one atom per unit cell in
the Debye approximation.
(a) Find the mean square displacement of an atom at absolute zero temperature.
(b) Derive a formula for the contribution of the zero-point motion to the total energy
of the crystal.
(c) Determine an expression for the phonon entropy.
(d) Evaluate (a), (b), and (c) numerically for the case of solid aluminum.
I.29. Lattice dynamics of a molecular chain. Consider the following model (see Fig. I.4)
for a chain of diatomic molecules, where a = 2b, M and m are masses of the atoms,
and α and β are harmonic force constants. In this problem, consider the longitudinal
phonons only.
(a) Construct the dynamical matrix and calculate the phonon frequencies, ωλ (k).
(b) Sketch the phonon dispersion relation and the phonon density of states.
(c) Apply the procedure discussed in the text on second quantization to the system
and show that with a suitable choice of a and a† the Hamiltonian becomes
1
Ĥ = n̂λk + h̄ωλk , (I.5)
2
λk
I.31. Phonons in 2D. Consider a 2D square lattice with lattice constant a and a basis of
three atoms per primitive cell.
(a) What are the total number of phonon branches for each of the following:
longitudinal acoustic, transverse acoustic, longitudinal optic, and transverse
optic?
(b) Assume that the acoustic modes are represented by the Debye model. What is
the maximum acoustic phonon frequency in terms of the lattice constant a and
the velocity of sound vs ?
(c) Find the contribution to the low-temperature heat capacity (within a dimension-
less coupling constant) from the longitudinal acoustic phonon branch.
Express your answer in terms of T, TD , and N, where T is the temperature, TD is the
Debye temperature, and N is the number of cells in the crystal.
I.32. Zero-point energy correction. Find a simple expression for the zero-point energy
correction to the total energy of a solid in the following two cases.
(a) Assume a three-dimensional Debye model with one atom per cell where all
acoustic modes are degenerate.
(b) Assume a three-dimensional Einstein model with one atom per cell and all three
modes degenerate.
I.33. Charged harmonic oscillator. Consider a charged harmonic oscillator in a constant
electric field F = (−e)E. The Hamiltonian is
p2 γ
H= + x2 + eFx. (I.6)
2m 2
1 λ2
H = ω(A† A + ) − , (I.7)
2 ω
√
where) ω = * γ /m, and give the expressions for λ and A.
(b) Find A, A† .
(c) Give an expression for the position operator X(t) and give a physical explanation
for its form.
I.34. Mean position of particles in a one-dimensional potential. A particle is bound in
a one-dimensional potential, V(x) = 12 mω2 x2 − Cx3 (for small x). Describe how the
mean position of the particle changes for the different eigenstates when C is small.
(Hint: use perturbation theory and show that anharmonic terms will cause expansion.)
PART II
ELECTRON INTERACTIONS,
DYNAMICS, AND RESPONSES
5 Electron dynamics in crystals
In the preceding chapters, we discussed the electron states within a single-particle picture
in a perfect crystal. This is in some sense equivalent to considering electrons in free space
in quantum mechanics; the interaction of the electron with the crystalline potential renor-
malizes the dispersion relation of the electrons. To understand many phenomena, such as
electrical transport in a solid, we would need to describe the behavior of electrons in crys-
tals when some non-periodic perturbation, either from an external applied field or a defect
structure, is added to the crystalline periodic field. In this chapter, we introduce in a more
rigorous fashion the concepts of group velocity v, effective mass m∗ , and the equation of
motion for the dynamics of the electron.
Denoting H0 as the perfect crystal Hamiltonian, the motion of an electron in an addi-
tional potential U(r) within a one-electron picture is given by the usual time-dependent
Schrödinger equation
∂ψ(r, t)
[H0 + U(r)] ψ(r, t) = ih̄ . (5.1)
∂t
The static perturbation U(r) may arise from external applied electric or magnetic fields, or
come from imperfections in the crystal, such as structural defects and impurities. In many
important and common situations, U(r) is slowly varying in space (compared to the atomic
separation) but becomes very large as r → ∞. For example, a uniform applied electric
field E gives rise to U(r) = −eE · r, and a uniform magnetic field B produces a vector
potential A = 12 B × r in the symmetric gauge. Potential perturbations caused by charge
impurities, dislocations, surfaces, etc. are also often static and slowly varying compared to
the interatomic distance, but they can be very large. These common situations cannot be
treated straightforwardly by standard perturbation methods since, particularly for uniform
applied fields, the perturbing potential diverges at large distance r in an extended crystal.
We develop here a formulation that will allow us to determine the behavior of crystal
electrons in an external potential, which is slowly varying in both space and time. The
framework is to consider that the problem of an electron in a perfect crystal is solved, and
that we transform the dynamics of an electron under the influence of an external potential
to an electron in free space with renormalized properties.
The basic assumption or approximation is that the electron moves entirely within the
states of a single band. This approach is named the effective Hamiltonian theory or approx-
imation. This is a reasonable approximation for the case where a potential that is nearly
101
102 5 Electron dynamics in crystals
Wavevector k
Figure 5.1 Dynamics of a crystalline electron in a slowly varying field. A slowly varying potential moves an electron in a state
with wavevector k to nearby states in the same band.
constant in space and smoothly varying in time moves the electron adiabatically between
states of nearly the same energy without causing transitions to states in another band, as
illustrated in Fig. 5.1.
Electron dynamics is best discussed in terms of localized wavefunctions because we will
be thinking in terms of wavepackets. One such set of functions is the so-called Wannier
functions.1
Let φnk (r) be Bloch functions of wavevector k in the band n; the Wannier functions are,
in the simplest form, defined by
1 −ik·Rj
wn (r − Rj ) ≡ √ e φnk (r), (5.2)
N k
where N is the number of unit cells and Rj is a lattice vector. This expression is for a
single non-degenerate band. It can be generalized to a band complex. Moreover, there is
considerable freedom in defining Wannier functions. Other possible definitions of Wannier
functions include adding a well-chosen phase factor eiθk in each term of the sum on the
right-hand side of Eq. (5.2) to maximize its localization in real space.2
Wannier functions defined by Eq. (5.2) form a complete set and have the following
useful properties:
1 ik·Rj
(i) φnk (r) = √ e wn (r − Rj ),
N j
(5.3)
(ii) wn (Ri )|wn (Rj ) = δn n δij ,
(iii) wn (r − Rj ) is localized about Rj .
1 G. H. Wannier, “The structure of electronic excitation levels in insulating crystals,” Phys. Rev. 52(1937), 191.
2 N. Marzari and D. Vanderbilt, “Maximally localized Wannier functions for composite energy bands,” Phys.
Rev. B 56(1997), 12847.
103 5.2 Electron dynamics in the effective Hamiltonian approach
n (r – Rj)
Rj
Figure 5.2 Schematic diagram of a Wannier function localized around the lattice site Rj .
The Wannier functions are therefore an orthonormal set of localized functions in real space
that span the same functional space as the Bloch states. We have used the Dirac bra–
ket notation above. A real-space picture of a Wannier function, given by Eq. (5.2) and
computed with planewave Bloch functions, is illustrated in Fig. 5.2.
We now consider the motion of an electron described by Eq. (5.1) in the effective Hamil-
tonian approximation. It is most easily discussed in the notation of the state vector |ψ(t)
for the electronic state. Mathematically, |ψ(t) is just a vector in Hilbert space. It can be
expanded in any complete set of “basis” functions that spans the Hilbert space. In partic-
ular, the Bloch states |φnk , the position states |r, or the Wannier functions |wn (R) may
all be used as separate sets of such basis functions.
For the three sets of basis functions of interest here, completeness of each set of functions
gives
1= |φnk φnk |
nk
= |r r| (5.4)
r
=
wn (Rj ) wn (Rj )
.
nRj
in expanding the state vector in the basis set of |r, called the r-representation. Here we
have denoted the expansion coefficients ψ(r, t) ≡ r|ψ(t) |r. This is just the real-space
wavefunction ψ(r, t) in elementary quantum mechanics as given in Eq. (5.1).
The state vector can equivalently be expanded in the basis set of Bloch states |φnk , and
in the Bloch representation we have
|ψ(t) = |φnk φnk |ψ(t)
nk
(5.6)
≡ ψn (k, t)|φnk ,
nk
where we have denoted the expansion coefficients φnk |ψ(t) |r as ψn (k, t). Similarly, in
an expansion in Wannier functions, we have
|ψ(t) = |wn (Rj )wn (Rj )|ψ(t)
nRj
(5.7)
≡ ψn (Rj , t)|wn (Rj ),
nRj
with ψn (Rj , t) ≡ wn (Rj )|ψ(t) |r denoting the expansion coefficients in the Wannier
representation.
We would now like to derive a more simplified equation of motion for the state vector
|ψ(t), starting with
∂|ψ
H|ψ = ih̄ , (5.8)
∂t
where H = H0 + U(r) is the Hamiltonian operator. Since U(r) is slowly varying in r and t,
we assume that the motion of the electron is well approximated by the Bloch states within
a single band. Let us call this the nth band. Then Eq. (5.7) in the effective Hamiltonian
approach is taken to be well approximated by
|ψ(t) = ψn (Rj , t)|wn (Rj ). (5.9)
j
The sum is now restricted to Wannier functions of a single band of index n. Combining
Eqs. (5.5) and (5.9), we see that
ψ(r, t) = ψn (Rj , t)wn (r − Rj ). (5.10)
j
The time evolution of the wavefunction ψ(r, t) is then equivalently given by the time evo-
lution of ψn (Rj , t). In this form, we see that the electron wavepacket (as a solution to
Eq. (5.1)) can be viewed as a set of Wannier functions convoluted by an envelope function
that depends only on the lattice site index Rj . The electron function ψ(r, t) may have a
lot of internal structure, but as far as the dynamics of the wavepacket are concerned, we
105 5.2 Electron dynamics in the effective Hamiltonian approach
only need to know the dynamics of the envelope function (or the wavefunction in the Wan-
nier representation), i.e. the time evolution of ψn (Rj , t) with respect to the crystal lattice
positions. We therefore need to develop quantum operators in the Wannier representation.
We start with the usual position rop and
) momentum
* pop operators that have the properties
rop |r = r |r , pop |p = p |p, and pop , rop = −ih̄, and define the operators Rop and
kop using the relations
kop |φn,k ≡ k|φn,k
and
Rop |wn (Rj ) ≡ Rj |wn (Rj ). (5.11)
Then, it is straightforward to show that Rop and kop form a conjugate pair of dynamical
variables with
[kop , Rop ] = −i, (5.12)
∂
Rop = i , (5.13)
∂k
and
∂
kop = −i . (5.14)
∂Rj
We now consider the action of a slowly varying potential U(r) on a Wannier function.
Since the potential is position dependent, it is well-defined as an operator acting on a func-
tion that is dependent on r. We may express the Wannier function in the r-representation
and obtain
Uop |w(Rj ) = Uop |rr|w(Rj ) = U(r)|rr|w(Rj ). (5.15)
r r
Since r|w(Rj ) is a very localized function around the point Rj , and U(r) is very slowly
varying in the scale of the spacing of the R’s, we may write
Uop
w(Rj ) ∼= U(Rj ) |rr|w(Rj ) = U(Rj )|w(Rj ). (5.16)
r
Equation (5.16) gives a great simplification to Eq. (5.8). Using the single-band
approximation, Eqs. (5.13) and (5.14),
∂
H0 |ψ ∼
= En (k)φn,k |ψ|φn,k = En −i w(Rj )|ψ |w(Rj ). (5.17)
∂Rj
k Rj
Equation (5.8) can now be rewritten in the Wannier representation, within the slowly
varying potential regime, as
∂
∂
En −i + U(R) ψ(R, t)|w(R) = ih̄ ψ(R, t)|w(R), (5.18)
∂R ∂t
R R
106 5 Electron dynamics in crystals
where, for simplicity, we have dropped the index j on the lattice positions. We then have
an effective equation for the wavepacket envelope function
∂ ∂
En −i + U(R) ψ(R, t) = ih̄ ψ(R, t). (5.19)
∂R ∂t
Effectively, we have the following Hamiltonian for considering the dynamics of the crystal
electrons:
∂
Heff = En −i + U(R) (R-representation) (5.20)
∂R
or
∂
Heff = En (k) + U i (Bloch representation). (5.21)
∂k
Within this formalism, by making use of the slowly varying potential approximation,
we have eliminated the electron coordinates r from the problem, and hence all the details
of the electron’s interactions with the crystal. The effects of the crystal on the behavior
of the electron are incorporated completely in the band structure term in Eq. (5.20). The
electron behaves as if it moves only in the potential U(R) but with a renormalized kinetic
energy operator in Eq. (5.20) or Eq. (5.21). Making a connection to the usual electron
wavefunction, we have
ψ(r, t) = ψ(R, t)wn (r − R) (5.22)
R
if Eq. (5.21) is used. The two forms are completely equivalent; which one to use is a
matter of computational or conceptual convenience. The effective Hamiltonian approach
is expected to be a good approximation if the spatial and temporal variations of U are such
that (a) (∇U) · d is much less than the bandwidth where d is an interatomic distance, and
(b) h̄ω is much less than the interband energies where ω is the time variation frequency of
U. Under these conditions, a state initially composed of Bloch states within one band will
not mix with states of other bands under the influence of the potential U.
The robustness of the approach is illustrated schematically in Fig. 5.3 using a constant
applied electric field as an example. The quality U(r) is the difference between the ex-
act external potential and the piece-wise constant potential taken at the lattice sites where
the Wannier functions are located. We see that U is small (∼eEd) and oscillatory. The
largeness (diverging at large distances) of the perturbation is entirely captured in U(R).
Moreover, since U is periodic, its effects may in fact be included by redefining the
Bloch and Wannier functions as solutions to a new band structure problem, i.e. a periodic
crystalline Hamiltonian with U included. Of course, the effects of U in the Wannier
functions in Fig. 5.3 are real physical effects that were left out of the theory, e.g. those
related to polarization.
107 5.3 Shallow impurity states in semiconductors
Wannier function
n
r
i
n
r
i+1
n
r
Ri–1 Ri Ri+1
(b)
U (r) = –e ε · r
average potentials
External and cell
Ũ (r) = –e ε · R(r)
Ri–1 Ri Ri+1
(c)
r
Ri–1 Ri Ri+1
Figure 5.3 Schematic of Wannier functions (a) at different lattice sites Rj and (b) the potential due to a uniform electric field and
(c) the difference between U(r) and Ũ(r). Here R̃(r) is a piece-wise constant function that takes the value of the
closest value of Rj at a given position r.
A good example of an application of the effective Hamiltonian approach is its use in de-
termining the shallow impurity states in semiconductors. Let us consider the simple case
of a substitutional impurity (with an extra valence electron from the impurity atom) in a
semiconductor that has an isotropic and parabolic conduction band. This would be, for
example, the case of Si in GaAs replacing one of the Ga atoms. The extra electron would
be attracted to the Si site and in a state that is expected to be comprised of states from the
conduction band only (see Fig. 5.4).
We use the single-band approximation and write the impurity wavefunction in terms of
the conduction band states:
ψ(r, t) = ψ(k, t)φck (r)
k
(5.24)
= ψ(R, t)wc (r − Rj ),
R
108 5 Electron dynamics in crystals
Conduction
band
Donor state
Eg
Valence
band
Wavevector k
Figure 5.4 Diagram of a shallow donor state in a direct gap semiconductor with a bandgap Eg separating the valence band from
the conduction band.
where φck (r) and wc (r − R) are the conduction band Bloch functions and Wannier func-
tions, respectively. The conduction band dispersion relation measured from the top of the
valence band, assuming an isotropic and parabolic form, is given by
h̄2 k2
Ec (k) = Eg + , (5.25)
2m
where Eg is the bandgap and m∗ is a parameter describing the band curvature. Application
of Eq. (5.20) leads to the following effective Schrödinger equation for the wavefunction in
the Wannier or R-representation, assuming a Coulomb interaction (screened by the dielec-
tric constant of the material) between the additional electron and the positively charged
impurity core site,
2
h̄2 2 e
− ∇ + − ψ(R) = (E − Eg )ψ(R). (5.26)
2m R
This leads to the hydrogenic model for shallow impurity states in the bandgap near the
conduction band of semiconductors given in elementary treatments. Table 5.1 compares
the calculated results using the effective Hamiltonian (or effective mass) model with
experimental impurity levels.
We now proceed to analyze the motion of crystal electrons in an external field using the
formalism developed in Section 5.2. The equation of motion for the operators Rop and kop
may be obtained using the usual relation from quantum mechanics for operator A,
dA ∂A 1
= + [A, H]. (5.27)
dt ∂t ih̄
109 5.4 Motion in external fields
Table 5.1 Binding energies Eb of the ground (1s) state of shallow donor
impurities in some common semiconductors, comparing experimental data
with results from the effective Hamiltonian or effective mass theory. The
symbol XY denotes that the element Y in the host material is substitutionally
replaced by an impurity atom X. (After Cardona and Yu, 2010.)3
Semiconductor Theory Experiment
Eb (meV) Eb (meV)
Using the effective Hamiltonian in Eqs. (5.20) and (5.21), we obtain for the operator R the
following equation of motion:
dR
ih̄ = [R, Heff ]
dt
∂
= R, En −i + U(R)
∂R
∂
= R, En −i
∂R
∂
= i , En (k)
∂k
∂ ∂
=i En (k) − En (k)i
∂k ∂k
∂
=i En (k) . (5.28)
∂k
3 P. Y. Yu and M. Cardona, Fundamentals of Semiconductors: Physics and Materials Properties, 4th ed. (Berlin:
Springer, 2010).
110 5 Electron dynamics in crystals
Note that v does not depend on the form of the perturbation U. It is a general property of
the electronic band structure.
For k, again using Heff , we obtain the equation of motion
dk
ih̄ = [k, Heff ]
dt
∂
= k, U i
∂k
∂
= −i , U(R)
∂R
= −i (∇R U(R)) (5.30)
or
dk
= −∇ R U(R).
h̄ (5.31)
dt
Since −∇U(R) is the force F on the electron, we have
dk
h̄ = F. (5.32)
dt
For the special case of a constant force F, Eq. (5.32) can be integrated immediately to
t
k(t) = k(t = 0) + F . (5.33)
h̄
With Eqs. (5.29) and (5.31), we can follow the motion of an electron starting at t = 0 in
any individual Bloch state with wavevector k0 , or of a wavepacket centered about a given
k0 . For example, consider a crystal with the band structure given in Fig. 5.5 in a uniform
electric field E = εx̂. The wavevector of a state with initial wavevector at k0 = 0 at t = 0
will change according to Eq. (5.33) to
|e|εt
k(t) = − x̂. (5.34)
h̄
∂En
The group velocity is given by h̄vx = ∂kx (see Fig. 5.6), which leads to a time dependence
in the position of the particle as
t t dt ∂Endt ∂En
1 t
x(t) = vx (t)dt = = dkx
0 0 h̄ ∂kx
0 dkx ∂kx
h̄
(5.35)
1 1 2
= [En kx (t)) − En (kx (t = 0) ].
(−|e|ε)
111 5.4 Motion in external fields
C′ C
B′ B
–G/2 0 G/2 G
Wavevector kx
Figure 5.5 Schematic of a band structure in a repeated zone scheme along the kx direction, where G is the smallest reciprocal
lattice vector along this direction.
B′ A B
kx
–G/2 0 G/2 G
Figure 5.6 Group velocity vx as a function of kx for the lower band in Fig. 5.5.
W
0 t Δ xmax =
eε
Figure 5.7 Schematic of x(t) of a crystalline electron in a uniform field which shows Bloch oscillation where W is the band width.
Physically, what these equations tell us is that the wavepacket’s center of mass wavevec-
tor in k-space and its center of mass in real space are changing (due to the interaction with
the field) as a function of time. When the center of mass wavevector reaches the Brillouin
zone boundary, the wavepacket is Bragg-reflected, leading to an oscillatory behavior in
r(t) as illustrated in Fig. 5.7. As a consequence, in this idealized situation, there is no net
current for a perfect crystal. However, for real materials, impurities and zero-point motions
of various excitations will introduce scattering that would give rise to a net current even at
temperature T = 0.
The oscillations in the motion of the electrons given by Eq. (5.35) are called Bloch
oscillations. Observation of these oscillations requires that the electron travel many periods
112 5 Electron dynamics in crystals
Band 2
eε
Eg
Band 1
W Eg
Δ xmax = d=
eε eε
Figure 5.8 Semiclassical picture of Bloch oscillations with E an applied uniform electric field and W the bandwidth.
W
xmax = , (5.36)
|eε|
where W is the bandwidth. Under normal conditions, with an electric field of the order of
V/mm and W ∼ few eV, xmax ≈ a few mm, making the Bloch oscillations unobservable
for usual samples since the mean free path of the electrons is significantly shorter. Bloch
oscillations, however, have been observed in man-made periodic systems, such as semi-
conductor superlattices and atoms in optical lattices, where the large unit cell dimensions
reduce W to small values.4
The Bloch oscillation phenomenon may also be understood in a semiclassical picture
as illustrated in Fig. 5.8. In a slowly varying potential, the electron in a state at the band
extrema effectively runs into a triangular potential barrier with width d that is determined
by the bandgap to the next band and by the electric field strength. The electron wave is
reflected at the barriers, giving rise to the oscillations. There is, however, a finite probabil-
ity for tunneling to the next band if |E| is large enough (i.e. if the barrier is thin enough),
leading to phenomena such as Zener tunneling (see Fig. 5.9). Interband transitions are,
however, not included in the present effective Hamiltonian formalism, since it is assumed
that the wavefunction is always made up of states from a given band or band complexes.
4 See, for example, M. Raizen, C. Salomon, and Q. Niu, “New light on quantum transport,” Phys. Today
50(1997), 30.
113 5.5 Effective mass tensor
Energy
Wavevector k
Figure 5.9 Schematic of Zener tunneling with an electron tunneling from the lower band to the upper band under the influence
of a strong external field.
In the example in Figs. 5.5 and 5.6, under the influence of a force F, the wavevector k(t)
of an electron changes as a function of time. The velocity of the electron v goes to zero and
then changes sign as k crosses the Brillouin zone edge under a constant force. The electron
appears as if it has a mass that is a function of k and becomes negative near the top of the
band. This leads us to define the useful concept of the effective mass tensor in electron
dynamics (discussed previously in Chapter 3).
As in classical mechanics, mass is defined by F = ma or, more generally,
1
Fβ = aα . (5.37)
m∗ αβ
β
dvα 1 d ∂E 1 ∂ 2 E dkβ
= = , (5.38)
dt h̄ dt ∂kα h̄ ∂kα ∂kβ dt
β
1 ∂ 2E
aα = Fβ . (5.39)
β
h̄2 ∂kα ∂kβ
The form and properties of the effective mass tensor play an important role in determining
the transport properties of materials, especially those of carriers in semiconductors.
In semiconductors, one is mainly dealing with carriers (either electrons or holes) near
band extrema, which typically occur near a symmetry point k0 in the Brillouin zone, e.g.
114 5 Electron dynamics in crystals
at k0 equal to the point for GaAs. There is a way, as discussed in Chapter 3, that allows
one to determine the nearby band structure and the effective mass tensor given that we
know the energy states at a given point k0 . The method makes use of the knowledge of
the electron band structure at a given point of k-space to find the states nearby, and the
effective mass tensor, in the so-called k · p perturbation theory, resulting in, for the nth
band,
1 1 ∂ 2 En 1 2 pαn n pnn
β
= 2 = δαβ +
m∗ n,αβ h̄ ∂kα ∂kβ m m En (k0 ) − En (k0 )
n
αβ
1
= δαβ + fnn . (5.41)
m n
We have used the notation pαn n = n k0 |pα |nk0 and defined the quantity f through the
second line of Eq. (5.41).
The equation of motion given in Eq. (5.29) is incomplete when we consider the effects of
geometric or Berry phase5 on the motion of crystal electrons. In the effective Hamiltonian
approximation, an electron in a Bloch state |nk with wavefunction
is constrained to move within a single band n. In the presence of a slowly varying static
potential U(r), the electron will adiabatically drift in k-space along a path C (see Fig. 5.10).
Denoting the electron to initially be in the state |ψ(t = 0) = |nk(0), then at some later
time t, the state of the electron will evolve into the state |ψ(t), which within the effective
Hamiltonian approximation is |nk(t) apart from a phase factor according to the adiabatic
theorem of quantum mechanics.6 (Here |nk(t) is given by Eq. (5.42).)
The phase factor that relates |ψ(t) to |nk(t) may be determined by the following
analysis. The cell periodic part of a Bloch wavefunction |unk satisfies the equation
(p + h̄k)2
+ V(r) unk (r) = En (k)unk (r), (5.43)
2m
or
H(k) |unk = En (k) |unk . (5.44)
5 M. V. Berry, “Quantal phase factors accompanying adiabatic changes,” Proc. R. Soc. Lond. A 392(1984), 45.
6 See, for example, Chapter 10 in D. J. Griffiths, Introduction to Quantum Mechanics, 2nd edn. (Upper Saddle
River, NJ: Prentice Hall, 2005).
115 5.6 Equations of motion, Berry phase, and Berry curvature
2π
b Ωn (k)
ky γn
C unk
0 kx 2π
a
Figure 5.10 Berry curvature in k-space. Here |unk is the periodic part of a Bloch state and n (k) is the Berry curvature of band n.
C denotes a closed loop over which an electron in the nth band picks up a Berry phase of γn .
We shall use the solutions |unk in Eq. (5.44) as an instantaneous orthonormal basis to
expand |ψ(t). Equation (5.44) allows, however, an arbitrary k-dependent phase factor of
|unk . To remove this arbitrariness, we make a phase choice, or a gauge, that requires the
phase of |unk to be smooth and single-valued along the path C in the parameter or k-space.
Treating k(t) as a parameter along the path C, the periodic part of the state of the electron
(denoted by φ(r, t) ≡ e−ik·r ψ(r, t)) evolves according to the time-dependent Schrödinger
equation
d
H k(t) |φ(t) = ih̄ |φ(t) . (5.45)
dt
Adiabatically, a system prepared in |nk(0) will evolve with H k(t) and therefore be in
the state |nk(t) at t.
So, |φ(t) can be expressed as
−i t 1 2
|φ(t) = exp dt En k(t ) exp iγn (t) |unk (t) . (5.46)
h̄ 0
The first exponential factor on the right-hand side of Eq. (5.46) is the familiar dynamical
phase factor. The second exponential contains a phase which is in general non-integrable
and cannot be written as a function of k.
We determine γn (t) by requiring that |φ(t) in Eq. (5.46) satisfy the time-dependent
Schrödinger equation. Substituting Eq. (5.46) into Eq. (5.45) leads to the following
expression for the time rate of change in γn (t):
d
dk
γn (t) = i unk(t)
∇ k unk(t) · , (5.47)
dt dt
where
un(k+dkk̂i ) − |unk
∇ k unk(t) = k̂i . (5.48)
dk
i
116 5 Electron dynamics in crystals
Integrating Eq. (5.47) along the path C from some point k0 with respect to time gives
k
γn (k) = i unk | ∇ k unk · dk . (5.49)
k0
The normalization of |unk dictates that unk | ∇ k unk = − ∇ k unk | unk , and therefore
γn (k) is real.
Let us define the Berry vector potential An (k) (or the Berry connection) as the quantity
We note that A has dimensions of length and is always real, as argued above. The elec-
tron’s ket vector |φ(t), during its adiabatic evolution along the path C in parameter space,
acquires an extra phase γn (k), in addition to the dynamical phase with
k
γn (k) = dk · An (k ). (5.51)
C
The form of Eq. (5.51) reveals that we may regard this phase as an Aharonov–Bohm phase
caused by an effective “magnetic field”
which lives in parameter space k. This field n (k) is called the Berry curvature.
The Berry vector potential An (k) is gauge dependent. If we make a gauge transformation
is gauge invariant and is physically observable. The phase γn in Eq. (5.55) is known as the
Berry phase or geometric phase. Although not gauge invariant, γn (k) in Eq. (5.51) is also
often loosely referred to as the Berry phase.
The1 Berry curvature
2 n (k) defined in Eq. (5.52) is also a gauge-invariant quantity, since
∇ k × ∇ k μ(k) = 0 for any smooth function μ(k). Using Stoke’s theorem, we may rewrite
the Berry phase for a closed loop in the form
5
γn = dk · An (k) = ds · n (k), (5.56)
C S
117 5.6 Equations of motion, Berry phase, and Berry curvature
where S is any surface enclosed by the closed loop C. Unlike the Berry vector potential,
since it is gauge invariant the Berry curvature is observable. The physical origin of the
Berry curvature arises from the adiabatic approximation employed, which is equivalent to
a projection operation restricting the dynamics of the system to the nth band. The Berry
curvature can be viewed as arising from the residual interaction with states in the bands
that are not considered in the effective Hamiltonian.
The Berry curvature n (k), an intrinsic local property in k-space of the band structure
of a crystal, modifies the motion of the electron. Let us consider the effective Hamiltonian
with external potential U(r) in the Bloch representation
where R = i∇ k is the Wannier coordinate indexing the lattice sites and En (k) is the
unperturbed band structure. We may perform a gauge transformation to remove the
effective Aharonov–Bohm phase or Berry phase in Eq. (5.46), at the expense of adding the
Berry vector potential An (k) to i∇ k in the argument of U in Eq. (5.57). The transformed
Hamiltonian H is then equal to
1 2
H = En (k) + U i∇ k + An (k) . (5.58)
x = i∇ k + An (k) (5.59)
and is no longer simply R. As an operator, the components of x do not commute with each
other. Their commutation relations are given by
)
xi , xj ] = i ijk nk (k), (5.60)
dk
h̄ = −i[k, H ] = −∇ R U(R), (5.61)
dt
which is unchanged from Eq. (5.31). However, the equation for the velocity is now given
by
h̄v = −i[x, H ] = ∇ k En (k) + ∇ x U × n (k). (5.62)
By comparison with Eq. (5.29), there is an extra term involving the Berry curvature. The
new term is commonly known as the Luttinger anomalous velocity.7 For example, for a
7 M. C. Chang and Q. Niu, “Berry phase, hyperorbits, and the Hofstadter spectrum,” Phys. Rev. Lett. 75(1995),
1348.
118 5 Electron dynamics in crystals
1 2
crystal in a static electric field E with U(x) = −eE · x and a weak magnetic field B, the
semiclassical equations of motion become
dk e
h̄ = eE + v × B (5.63)
dt c
and
h̄v = ∇ k En (k) − eE × n (k). (5.64)
The importance of the contribution of the Berry curvature term to the dynamics of crystal
electrons depends on the systems of interest. In many cases, n (k) may be neglected.
However, n (k) is nonzero in a wide range of materials, particularly those with broken
time-reversal and inversion symmetries, and has given rise to insightful understanding of a
number of interesting phenomena, including the anomalous Hall effect, the quantum Hall
effect, electric polarization, orbital magnetization, topological insulators, etc.
Let us analyze the general form of the Berry curvature n (k) by examining the velocity
expression in Eq. (5.64). The velocity expression should be invariant under time-reversal
or spatial inversion operations if the crystal has these symmetries. Under time reversal,
v → −v, k → −k, and E → E; under spatial inversion, v → −v, k → −k, and
E → −E. Thus, if the crystal has time-reversal symmetry, the invariance form shown in
Eq. (5.64) dictates that
n (−k) = −n (k), (5.65)
and, if the crystal has spatial inversion symmetry,
This leads to the conclusion that, for crystals (such as Al) with both time-reversal and
spatial inversion symmetry, the Berry curvature is identically zero throughout the Brillouin
zone, with the possible exception of nonzero values at singular points in k-space due to
band degeneracy. In such a case, one may use the simpler velocity expression given in Eq.
(5.29).
For systems with either broken time-reversal or broken spatial inversion symmetry, the
proper equation for the velocity should include the Berry curvature term. In ferromagnetic
materials, where time-reversal symmetry is broken, such as fcc Fe under a uniform E-
field, the Berry curvature term in Eq. (5.64) gives rise to a current that is perpendicular
to E and constitutes an intrinsic mechanism for the so-called anomalous Hall effect (in
which a strong additional Hall resistivity results from the magnetization) in these materials.
Presently, this research area is dynamic and evolving.8
8 See, for example, the review article by D. Xiao, M. C. Chang, and Q. Niu, “Berry phase effects on electronic
properties,” Rev. Mod. Phys. 82(2010), 1959.
Many-electron interactions: the homogeneous
6 interacting electron gas and beyond
Up until now, we have treated the electronic Hamiltonian of solids as an effective Hamil-
tonian of independent electrons moving in a mean (or average) field due to the ions and
other electrons. This is a very useful conceptual approach and can be quantitative for cer-
tain physical properties if the appropriate field or self-consistent potential is used, as we
shall see in the next chapter on density functional theory.
The proper choice of the mean-field potential V(r) appearing in the one-electron
Schrödinger equation
3 4
h̄2 2
− ∇ + V(r) ψ(r) = εψ(r) (6.1)
2m
is a complex and subtle issue. Underlying this problem is the question of how best to
approximate the effects of electron–electron interactions in the exact Hamiltonian H within
the Born–Oppenheimer approximation for physical quantities of interest:
3 4
p2i 1 e2
H= + Vion (ri − R) + , (6.2)
2m 2 |ri − rj |
i R i=j
where Vion (ri − R) is the ionic potential seen by the ith electron. Very often, an effec-
tive one-electron Schrödinger-like equation of the form Eq. (6.1) is obtained by seeking
the ground-state energy and wavefunction of an interacting many-electron system using a
variational principle approach.
In the Hartree–Fock (HF) treatment, for example, we assume that the ground-state many-
electron wavefunction is of the form of a single Slater determinant:
HF (r1 σ1 , r2 σ2 , ..., rN σN )
φ1 (r1 σ1 ) φ1 (r2 σ2 ) φ1 (r3 σ3 ) ... φ1 (rN σN )
1
φ2 (r1 σ1 ) φ2 (r2 σ2 ) φ2 (r3 σ3 ) ... φ2 (rN σN )
=√
.. .. .. ..
. (6.3)
N!
. . . .
φ (r σ ) φ (r σ ) φN (r3 σ3 ) ... φN (rN σN )
N 1 1 N 2 2
Here N is the number of electrons in the system, and ri and σi denote the spatial and spin
coordinates of the ith electron, respectively. φλ are single-particle orbitals to be determined
in the calculation. Since an electron has spin 12 , σi may be up or down. Applying the
variational principle, i.e. demanding the variation of the ground-state energy with respect
119
120 6 Many-electron interactions: the homogeneous interacting electron gas and beyond
δHF |H|HF
=0 (6.4)
δφλ
h̄2 2
− ∇ φi (r, σ ) + Vion (r)φi (r, σ ) + VH (r)φi (r, σ )
2m
(6.5)
+ Vx (r, r , σ , σ )φi (r , σ )dr = εi φi (r, σ ),
σ
where
e2
occ
2
VH (r) = dr φ j (r , σ ) , (6.6)
|r − r |
σ j
and
occ
e2
Vx (r, r , σ , σ ) = − φ ∗ (r , σ )φj (r, σ ). (6.7)
|r − r | j
j
An important model system that illustrates the role of many-electron interactions in solids
without having to deal with the complexity of the atomic structure is that of the uni-
form interacting electron gas or the jellium model. In this model, we consider a system
of interacting electrons in a uniform positive-charged background so that the total system
is neutral. The jellium model is then an idealized metal for which we can study the effects
of electron interactions on the total energy, electron (or quasiparticle) excitations, and other
properties as a function of electron density.
In the zero-order approximation, we neglect electron–electron interactions; the jellium
model reduces to that of a non-interacting charged Fermi gas, i.e. the Sommerfeld model.
In the non-interacting limit, the Hamiltonian can be written as
p2
i
H0 = + VPB = Hi0 , (6.9)
2m
i i
where VPB is the potential due to the positive background. For an N electron paramagnetic
system, we may take the ground-state wavefunction to be in the form
6
N
0 = φi (ri , σi ), (6.10)
i=1
with
1
φi (r, σ ) = φkλ
0
(r) = √ eik·r χλ , (6.11)
where is the volume of the sample and χλ is a spin function with λ = ↑ or ↓. The value
of k is within the Fermi sphere, i.e. k ≤ kF , with kF the Fermi wavevector defined through
N=2 1, (6.12)
k≤kF
where the factor 2 is due to spin, and the orbital energy is independent of λ given by
h̄2 k2
ε0 (k) = .
2m
The total energy of the system is
E0 = 2 ε0 (k). (6.13)
k≤kF
The jellium model has only one variable physical parameter, the density of the electrons
n=
N
or kF = (3π 2 n)1/3 in the paramagnetic phase. In the study of the interacting electron
122 6 Many-electron interactions: the homogeneous interacting electron gas and beyond
gas, it is often useful to express physical quantities in terms of another related parameter
rs (the typical interelectron distance in units of the Bohr radius a0 ):
1 1
3 3 1 3 −1
rs = a0 . (6.14)
4π n
The dimensionless parameter rs is known as the electron gas parameter and may be viewed
as the effective radius of the volume occupied by the electron in units of the Bohr radius.
In the literature, the density dependence of physical quantitites is traditionally expressed
in terms of rs . For example, the Fermi wavevector kF = 1.92/rs (a−1 0 ), the Fermi energy
3/2
EF = h̄ kF /2m = 3.68/rs (Ry), and the plasma frequency h̄ωp = 3.46/rs (Ry). Here
0 2 2 2
Ry denotes the energy unit of one Rydberg. In the Sommerfeld model, the average kinetic
energy per electron, EKE , is given by
2 0 3 2.21
EKE = εk = E0F = 2 (Ry). (6.15)
N 5 rs
k<kF
The average potential energy per electron EPE in the interacting system can be estimated
as
e2 2
EPE ≈ = (Ry). (6.16)
rs a0 rs
Thus, rs ∼ EPE /EKE is a measure of the ratio of the average potential energy to the
kinetic energy in the electron gas. A system with rs < 1 corresponds to one with high
density and dominant kinetic energy; a system with rs > 1 corresponds to one with low
density and dominant potential energy.
Since rs is a measure of the ratio of the Coulomb energy to the kinetic energy of
the electron gas, we may use it as a small parameter (if rs is small) to do perturbation
theory to determine the effects of Coulomb interaction on the properties of the elec-
tron gas. However, real materials typically have values of rs between 1 and 5 so that
EKE ∼ EPE , making the system strongly interacting, and making it difficult to ap-
ply traditional perturbation theory. Figure 6.1 illustrates the range of rs . The importance
of treating many-electron effects going beyond Hartree–Fock will generally depend on rs
and the dimensionality of the system under consideration. For example, in three dimen-
sions and for rs > 100, electron correlation is so strong that the electrons crystallize into a
crystalline arrangement forming an electronic lattice (the Wigner crystal).1 As we will see
below, one can determine the exact correlation energy in the limit of extremely large rs .
1 E. Wigner, “On the interaction of electrons in metals,” Phys. Rev. 46(1934), 1002.
123 6.2 Hartree–Fock treatment of the interacting electron gas
KE PE KE PE
1 10 100
rs
Real materials 3D Wigner e– crystal
Figure 6.1 Schematic of the range of the electron gas parameter rs on log scale showing the ranges where the kinetic energy (KE)
dominates over the potential energy (PE) and vice versa. Typical metals in three dimensions have rs between 1 and 6.
At extremely high rs , the electrons form a Wigner lattice.
To gain some insight into electron interaction effects, we solve the uniform interacting
electron gas problem analytically within the Hartree–Fock approach. We shall assume that
the single-particle orbitals in Eq. (6.3) are planewaves times a spin function, i.e. the orbitals
are specified by the quantum numbers k and λ with
1
φiHF = φkλ = √ eik·r χλ , (6.17)
1 2 1 2
where χλ is the spin function 10 for λ = ↑ and 01 for λ = ↓, and show that φkλ are the so-
lutions of the self-consistent Hartree–Fock equations (6.5)–(6.7). We note that planewave
solutions result in a uniform electron charge density. This results in Vion (r) + VH = 0,
since the electron charge density cancels the positively charged background density. Only
the exchange term Vx given by Eq. (6.7) survives among the potential energy terms in
Eq. (6.5).
Equation (6.5) reduces to
h̄2
− ∇ 2 φi (r, σ ) + Vx (r, r , σ , σ )φi (r , σ )dr = εi φi (r, σ ). (6.18)
2m σ
h̄2 k2 1 e2
φkλ (r) − 3/2 dr |
e−ik ·r eik ·r eik·r δλ λ χλ = ελ (k)φkλ (r).
2m
|r − r
k <kF λ
(6.19)
The second term on the left-hand side, the exchange term, when summed over λ may
be rearranged (assuming equal occupation of k ↑ and k ↓ orbitals, i.e. a paramagnetic
124 6 Many-electron interactions: the homogeneous interacting electron gas and beyond
We have thus shown that the planewaves φkλ are solutions to the Hartree–Fock equation
in the case of the uniform electron gas,
with eigenvalues (we shall drop the subscript λ from now on since we are considering a
paramagnetic system)
h̄2 k2 dk 4π e2
εHF (k) = − , (6.23)
2m k <kF (2π ) |k − k |2
3
Figure 6.2 illustrates the function F(x) and the Hartree–Fock energy eigenvalues for a
paramagnetic electron gas. If one invokes Koopmans’ theorem, then the eigenvalues are
a rough estimate of the electron (or quasiparticle) excitation energy of the system. The
eigenvalues describe the energy for exciting a single-particle-like excitation within the
Hartree–Fock approximation, including only the exchange interaction. The single-particle
energy can be written as a sum of two terms,
(a)
1 Infinite
slope at
x=1
F(x)
0.5
x
0 1 2
(b)
Free
electrons eHF (k)
1 1.5
Energy (arb. units)
k / kF
Infinite slope
at k = kF
2e2
eHF (k = 0) kF
π
Figure 6.2 (a) The function F(x) given in Eq (6.25) with infinite slope at x = 1. (b) Electron energy dispersion relation of a free
electron gas and an interacting electron gas in the Hartree–Fock (HF) approximation in the paramagnetic phase.
with
2e2 k
xHF (k) = − kF F . (6.27)
π kF
The quantity x given in Eq. (6.27) is called the self energy of the excited particle. Owing
to this exchange self energy, the single-particle states are lowered in energy and the occu-
pied bandwidth WHF = εHF (kF ) − εHF (0) increases from the independent-particle value of
h̄2 k2F /2m to h̄2 k2F /2m + e2 kF /π .
From the functional form of F(x), we note that F(1) = 1/2 and its derivative is
d 1 1 + x2
1 + x
F(x) = ln
−2 , (6.28)
dx 4x x 1 − x
126 6 Many-electron interactions: the homogeneous interacting electron gas and beyond
1 ∂εHF (k)
vHF (k) = → ∞ as k → kF . (6.29)
h̄ ∂k
The effective mass in the vicinity of kF can be calculated using the relation
m ∂εHF
= ,
m∗HF ∂ε0 (k)
which gives
m m ∂ HF
=1+ 2 x (k). (6.30)
m∗HF h̄ k ∂k
Thus, the effective mass m∗HF → 0 as k → kF . Both behaviors are unphysical. Therefore
single-particle excitation calculations based on the Hartree–Fock approximation typically
give poor results for metallic systems. Some other consequences of the singular behavior
of F(x) at x = 1 include a vanishing of the electron density of states D(ε) at ε = εF and an
electronic specific heat that behaves like Ce (T) ∼ ln|T|
T
instead of a linear T dependence
at low temperature T.
The above unphysical behavior can be traced back to the divergence of the Fourier
e2 2
transform 4πk2
of the Coulomb interaction er , as k goes to zero. The divergence reflects
2
the long-range nature of the Coulomb interaction. Physically, in a metallic system, er will
be screened by the other electrons, especially at large values of r. A screened Coulomb
potential would eliminate the divergence as k → 0. For example, within Thomas–Fermi
screening, the Coulomb potential between two electrons is given by (see Chapter 8)
e
V(r) = e−Ks r , (6.31)
r
2
which would yield a non-divergent V(k) = (k4π2 +K
e
2 ) , where Ks is the Thomas–Fermi screen-
s
ing wavevector. In general, a treatment going beyond the Hartree–Fock approximation is
required to obtain accurate excited-state electronic properties.
The above discussion shows that a Slater determinant of planewave single-particle or-
bitals indeed solves the Hartree–Fock equations. However, this solution need not be the
lowest energy or the ground-state solution. Overhauser showed that, at low densities (i.e.
high rs ), there are solutions in the form of charge or spin density waves that are lower in
energy than the uniform density planewave solution.
The Hartree–Fock approximation, being a variational approach for the ground-state wave-
function, is expected to be more accurate for the ground-state energy as opposed to
individual single-particle excitation energies. The ground-state energy per particle for an
127 6.3 Ground-state energy: Hartree–Fock and beyond
where EKE is the kinetic energy per particle, which for our planewave solutions is the same
as the non-interacting kinetic energy, and EX is the exchange energy per particle with (for
the paramagnetic case)
1 2 HF
EX = x (k). (6.34)
2 N
k<kF
Since the Hartree–Fock wavefunction is a Slater determinant of planewaves, the kinetic en-
ergy is the same as the Sommerfeld model (i.e. the non-interacting case) given by Eq. (6.15)
for an electron gas of density given by rs . The exchange energy may be evaluated using
Eqs. (6.27) and (6.34):
1 HF 1 2e2 k
EX = x (k) = dk − kF F
N N (2π )3 π kF
k<kF
(6.35)
1 e2 k4F 1 2 3 e2 kF
=− x F(x)dx = − .
n π3 0 4 π
The ground-state energy per particle, as expressed in terms of the electron gas parameter
rs , is equal to
2.21 0.916
GS =
EHF − (Ry). (6.36)
r2s rs
For rs = 2, which is typical of metals like Al, we see that the kinetic energy is nearly the
same as the exchange energy, showing that the system is not weakly interacting in this
sense.
The Hartree–Fock method yields two terms for the expression of the electron excita-
tion energies and ground-state energy per particle, given by Eq. (6.26) and Eq. (6.36),
respectively. We expect that a more accurate theoretical treatment of the many-body wave-
function, going beyond Hartree–Fock, would result in more terms in both expressions. In
the limit of small rs (i.e. high density), the kinetic energy dominates, and we may expect
the ground-state energy per particle (and other quantities) to be expandable in a series in
powers of rs . It is reasonable to expect EGS to be of the form
2.21 0.916
EGS = − + A + B ln rs + Crs + . . . (Ry) (6.37)
r2s rs
Using the Rayleigh–Schrödinger perturbation theory for the uniform interacting electron
gas, Gell-Mann and Brueckner showed in 19572 that indeed EGS may be expanded in the
above form with A = –0.096 and B = 0.0622. The portion of the energy going beyond
Hartree–Fock is traditionally called the correlation energy. Thus, the expression for EGS is
usually divided into
EGS = EHF
GS + EC = EKE + EX + EC . (6.38)
in units of Rydberg.
Figure 6.3 depicts the calculated correlation energy of the interacting electron gas as
a function of rs , using the results of many-body perturbation theory given by Eq. (6.39).
Since the result is derived from perturbation theory, it is accurate only in the limit of rs < 1.
Indeed, a better wavefunction (i.e. a better description of electron–electron interaction ef-
fects) can only lower the ground-state energy from that of Hartree–Fock. The correlation
energy should always be negative. The positive values of EC in Eq. (6.39) and Fig. 6.3 at
large rs are a sign of the breakdown of perturbation theory. For real metals with 2 < rs < 6,
the expression for EC in Eq. (6.39) is inaccurate.
For intermediate values of rs , there are several formulas for EC available in the literature.
One of the simplest and yet highly accurate expressions for EC is the Wigner interpolation
formula. Wigner was able to arrive at a formula that is good for the range of rs of real mate-
rials by using simple and sound physical reasoning. He considered the interacting electron
gas at very low density, i.e. rs 1. In this limit, the potential energy dominates over
the kinetic energy and the electrons condense into an electron crystal in three dimensions
at rs ≥ 100. (Such Wigner crystals have been observed in dilute electron systems, such
as electrons on a liquid He surface.) To calculate the ground-state energy of the electron
crystal, let us divide the crystal into Wigner–Seitz cells. For a close-packed structure, the
Wigner–Seitz cell may be replaced by a sphere of volume
4π 3
VWS = r . (6.40)
3 s
The energy per electron is the energy of the electron in one of the Wigner–Seitz spheres,
which consists of a kinetic energy term EKE and a potential energy term EPE . EKE comes
−3/2
from the zero-point motion of the electron and scales like rs , which may be neglected
in the rs → ∞ limit. The potential energy (which by definition is the exchange plus
correlation energy) comes from the Coulomb interaction of the charges involved and has
two contributions:
2 M. Gell-Mann and K. A. Brueckner, “Correlation energy of an electron gas at high density,” Phys. Rev.
106(1957), 364.
129 6.3 Ground-state energy: Hartree–Fock and beyond
Perturbation
theory
Figure 6.3 Correlation energy of the interacting electron gas as a function of rs from different calculations.
Here Ee−p is the interaction energy between the electron at the center of the Wigner–Seitz
sphere and the positive neutralizing background density np = 3/(4π r3s ), and it is given by
e2 3 3
Ee−p = dr − 3
= − (Ry). (6.42)
WS r 4π r s r s
The second term on the right-hand side of Eq. (6.41), Ep−p , is the self energy of the positive
neutralizing background charge, given by
6
Ep−p = (Ry). (6.43)
5rs
1.8
EPE = − (Ry). (6.44)
rs
Using the fact that EPE = EX + EC and the Hartree–Fock result of EX = −0.916/rs , we
then deduce that in the rs → ∞ limit, the correlation energy is
0.88
EC (rs 1) = − (Ry). (6.45)
rs
0.88
EC (rs ) = − (Ry). (6.46)
rs + 7.8
The Wigner interpolation formula has been proved to be accurate and useful in many
applications, and although many other correlation energy formulas have been derived over
the years, the Wigner interpolation formula remains among the best. Very accurate values
of EC for the electron gas over a wide range of rs of interest (rs ≈ 0−100) have now
3 There was an error in the original 1934 Wigner paper. The correct formula is given here.
130 6 Many-electron interactions: the homogeneous interacting electron gas and beyond
been obtained numerically using modern computation techniques such as quantum Monte
Carlo simulations.4 Figure 6.3 shows a comparison of the Wigner formula with results
from perturbation theory and quantum Monte Carlo simulations of the interacting electron
gas.
Before we continue the discussion of the exchange-correlation energy Exc and use it in
applications, let us understand more of its physical origin. For this purpose, we introduce
the concept of a pair-correlation function in a quantum many-body system.
We consider an N-body fermion wavefunction (x1 , x2 , x3 , . . . , xN ), where xi ≡ (ri , σi ),
and is normalized such that
|(x1 , x2 , . . . , xN )|2 dN x = 1. (6.47)
4 D. M. Ceperley and B. J. Alder, “Ground state of the electron gas by a stochastic method,” Phys. Rev. Lett.
45(1980), 566.
131 6.4 Electron density and pair-correlation functions
Here spin denotes integration over the spin degree of freedom of particles 1 and 2. We
may denote
9 :
ρ2 (r, r ) = |(rσ , r σ , x3 , . . . , xN )|2 dx3 . . . dxN , (6.51)
spin
N(N − 1) ρ2 (r, r )
g(r, r ) = . (6.52)
2 ρ(r)ρ(r )
A number of general properties of the pair-correlation function can be seen from the
above expressions. For example,
N(N − 1)
g(r, r )ρ(r)ρ(r )drdr = (6.53)
2
and
The sum rule (Eq. (6.53)) arises because the integrand is the distribution of particle pairs;
integrating it over all volume gives the total number of pairs in the system, which is
N(N − 1)/2. The reciprocity relation (Eq. (6.54)) comes from the symmetry of the first
line in Eq. (6.50).
Physically, the density of electrons at r would be uncorrelated to the occupation of an
electron at r in a material if r were very far away from r. Thus we have the relation
g(r, r → ∞) = 1. (6.55)
Also, we may decompose the pair-correlation function into parallel and antiparallel spin
components,
by separating the integration over spin coordinates in Eq. (6.50) into σ = σ and σ = σ
partial sums, i.e.
N(N − 1)
g↑↑ (r, r ) = |(rσ , r σ , x3 , . . . , xN )|2 δσ σ dx3 . . . dxN (6.57)
2ρ(r)ρ(r )
σ ,σ
and
N(N − 1)
g↑↓ (r, r ) = |(rσ , r σ , x3 , . . . , xN )|2 (1 − δσ σ )dx3 . . . dxN . (6.58)
2ρ(r)ρ(r )
σ ,σ
which is a consequence of the Pauli exclusion principle since the wavefunction in Eq.
(6.57) vanishes as r → r for σ = σ . However, the value of g↑↓ at r = r will depend on
the details of electron correlations in a specific system. Another useful sum rule is
) *
dr g(r, r ) − 1 ρ(r ) = −1. (6.60)
(since
g(r, r )ρ(r ) by definition is the density at r given an electron at r), and that
ρ(r )dr = N. Equation (6.60) is then the statement that (N − 1) − N = −1.
For a homogeneous system like an electron gas, the pair-correlation function g(r, r )
depends only on the distance between r and r and not on the positions separately:
g(r, r ) = g(
r − r
) = g↑↑ (
r − r
) + g↑↓ (
r − r
). (6.62)
We now evaluate these correlation functions for the interacting electron gas within the
Hartree–Fock approximation to gain some insight into their behavior. For simplicity, we
work with a paramagnetic electron gas, i.e. the number of spin-up electrons is equal to
the number of spin-down electrons. (We note that it can be shown that at sufficiently low
density, the electron gas goes into a magnetic phase.) The electron density given by Eq.
(6.49) for a homogeneous system is then a constant ρ(r) = ρ0 = N
and the pair-correlation
function is
N(N − 1)
g(|r − r |) = ρ2 (|r − r |)
2ρ02
(6.63)
N(N − 1)
2
= 2
(rσ , r σ , x3 , . . . , xN )
dx3 . . . dxN .
2ρ0 spin
The sum over (i, j) in Eq. (6.64) is carried out over all occupied orbitals.
For the homogeneous electron gas,
eik·r
φi (rσ ) = √ χsi (σ ), (6.65)
133 6.5 g(r, r ) of the interacting electron gas
(We note here that si is the quantum number denoting the spin state and σ is the coordinate
in spin space.) Thus,
9
2 :
1
φki si (rσ ) φki si (r σ )
g(
r − r
) =
φk s (rσ ) φk s (r σ )
. (6.67)
2ρ02 k s ,k s j j j j spin
i i j j
and obtain within Hartree–Fock (using N↑ = N↓ = N2 and noting that there are four
combinations for (si , sj ): (↑, ↑), (↓, ↓), (↑, ↓), and (↓, ↑))
2
1
occ
1 2 1 N 1
gHF
↑↓ r, r = 2
2 = 2
2 = , (6.69)
N N 2 2
ki ,kj
= 1 − f(r − r )2 , (6.70)
2
where
2 ik·r
occ
3 3
f(r) = e = [sin (rkF ) − (rkF ) cos (rkF )] = j1 (rkF ), (6.71)
N
k
(rkF ) 3 rkF
and j1 (rkF ) is the first Bessel function. We see explicitly that gHF
↑↑ depends only on |r − r |.
Since f(r) in Eq. (6.71) goes to 1 as r → 0, we have
↑↑ (|r|) → 0.
lim gHF (6.72)
r→0
1
↑↑ (|r|) →
lim gHF . (6.73)
r→∞ 2
Combining Eqs. (6.69) and (6.70) yields
134 6 Many-electron interactions: the homogeneous interacting electron gas and beyond
(a)
1
2
0 r kF
(b)
1
2
0 r kF
(c)
1
2
0 r kF
Figure 6.4 Pair-correlation function of a homogeneous electron gas in the paramagnetic phase in the Hartree–Fock
↑↓ (r), (b) g↑↑ (r), (c) g (r) = g↑↓ (r)+g↑↑ (r).
approximation: (a) gHF HF HF HF HF
2
1 3
gHF (r) = gHF
↑↑ (r) + g↑↓ (r) = 1 −
HF
j1 (rkF ) , (6.74)
2 rkF
with
1
gHF (r → 0) =
2
and
gHF (r → ∞) = 1,
in Fig. 6.4. We see that they have all the general properties expected of the pair-correlation
1
functions of a many-body fermionic system. The extent over which gHF ↑↑ deviates from 2
is dictated by k1F . On the other hand, gHF 1
↑↓ is a constant of value 2 because the Hartree–
Fock Slater determinant ground-state wavefunction only has exchange correlation between
electrons of like spin; there is no correlation between electrons of unlike spin. Thus, for
the paramagnetic phase, the chance of finding a spin-down electron at any value of r
is independent of whether there is a spin-up electron at r and is equal to 12 within this
approximation.
135 6.6 The exchange-correlation hole
From the definition of the pair-correlation function in Section 6.4, the electron density at
position r as seen by an electron at position r is given by
1 2
ρr (r ) = g r, r ρ(r ). (6.75)
The difference between the electron density seen by an electron at r, ρr (r ), and that of the
electron density ρ(r ) measured, for example, in an X-ray scattering experiment, is given
by
) * ) 1 2 *
ρr (r ) = ρr (r ) − ρ(r ) = g r, r − 1 ρ(r ). (6.76)
Integrating over volume and making use of the sum rule (Eq. (6.60)) yields
) 1 2 *
ρr (r )dr = g r, r − 1 ρ(r )dr = −1. (6.77)
Equation (6.77) tells us that there is a depletion of electron density near the position of
an electron, and that the integrated missing density is exactly equal to one electron. This
is physically expected, since the Pauli exclusion principle, as well as Coulomb repulsion,
reduces the density of electrons near a specific site already occupied by an electron. Con-
servation of particle number requires that the missing density be equal to one, since there
is already one electron at site r. This depletion of (or hole in) electron density given by
Eq. (6.76) is called the exchange-correlation hole.
The exchange-correlation hole moves with the electron and, in general, changes size
and shape as the electron travels through a system of inhomogeneous electron density.
However, Eq. (6.77) always holds; that is, the integrated depletion is that of a missing
electron. Thus, the electron density ρr (r ) seen by a given electron at r is always different
from that of ρ(r ). For example, in an interacting homogeneous electron gas, although
ρ(r ) is a constant, ρr (r ) is highly inhomogeneous. The structure of the inhomogeneity
is dictated by the size and shape of the pair-correlation function. Figure 6.5 illustrates
schematically the difference between the two densities, ρ(r ) and ρr (r ).
The evaluation of the pair-correlation function and the exchange-correlation hole for
a real material is non-trivial since it requires knowing the many-electron ground-state
wavefunction. In particular, a reasonable description of g↑↓ (r, r ) means incorporation of
136 6 Many-electron interactions: the homogeneous interacting electron gas and beyond
r r
r r
Figure 6.5 Schematic of ρ(r ) and ρr (r ). (a) charge density ρ(r ) of a system as measured by an external probe. (b) charge
density ρr (r ) of the same system as seen by an electron at r.
correlation effects going beyond Hartree–Fock. One approach is to use variational quantum
Monte Carlo simulations to determine the ground-state wavefunction. Figure 6.6 shows
some results from such calculations for diamond. Different panels depict g↑↑ (r, r ) and
g↑↓ (r, r ) as functions of r on the (110) plane, with r fixed either on the bond center or the
interstitial region of the diamond crystal. The pictures show a suppression from the value
of 12 for g↑↑ , as r → r , reaching zero at r = r , showing a deep exchange hole. On the
other hand, the suppression of g↑↓ is not nearly as pronounced because diamond is only a
moderately correlated electron system. Both g↑↑ (r, r ) and g↑↓ (r, r ) show anisotropy and
strong variations as a function of r, owing to the electron density inhomogeneity of the
material.
With the exchange-correlation hole, the potential energy for the same electron at r is
ρr (r )dr g(r, r )ρ(r )dr
ε(r) = e2 = e2 . (6.79)
|r − r | |r − r |
137 6.7 The exchange-correlation energy
(a) 0.45
0.40
0.35
0.25
0.10
0.05
(b)
0.45
0.55 0.55
0.40
(c)
0.45
0.45
0.40
0.35
0.30
0.05
0.10
0.15
0.20
0.25
(d)
0.45
0.40
0.35
0.30
Figure 6.6 The pair-correlation function in diamond for (a) parallel spin and (b) opposite spin with an electron at the bond
center, and for (c) parallel spin and (d) opposite spin with an electron at the tetrahedral interstitial site. The atoms and
bonds are schematically represented by the lines and dots in the [111] direction. (After Fahy, Wang and Louie, 1990.)5
5 S. Fahy, X.W. Wang, and S.G. Louie, “Pair-correlation function and single-particle occupation numbers in
The energy difference, which is the exchange-correlation energy for the electron at r, is
) * ) 1 2 *
ρr (r ) − ρ(r ) g r, r − 1 ρ(r )
εxc (r) = ε(r) − ε0 (r) = e2
dr = e2
dr .
|r − r | |r − r |
(6.80)
The exchange-correlation energy εxc (r) of an electron at r is then just the interaction
energy of the electron with its exchange-correlation hole, which can be considered pos-
itively charged with an integrated charge of −e. The quantity εxc (r) is also called the
exchange-correlation energy density of the system. The exchange-correlation energy of
the total system is
Exc = εxc (r)ρ(r)dr. (6.81)
We may apply this concept to the case of the interacting electron gas within the Hartree–
Fock approximation. Since the system is homogeneous, εxc (r) is the same everywhere and
is the exchange-correlation energy per particle. Using Eqs. (6.80) and (6.74),
∞ j2 (x)
dr ) HF * 6e2 kF 3 e2 kF
εxc
HF
= e ρ0
2
g (r) − 1 = − 1
dx = − , (6.82)
r π 0 x 4 π
which is the result we obtained before in Eq. (6.35). More accurate theory will yield a
better description of the shape of the exchange-correlation hole as well as a change in the
electron kinetic energy, and will therefore give a more accurate value for εxc (r). However,
the basic concept is the same as developed here.
As seen in Fig. 6.6, the shape and size of the exchange-correlation hole change as an
electron moves to different positions in a material. A particularly interesting situation is
the case in which the electron leaves the crystal and stays near the surface. The exchange-
correlation hole by necessity stays behind in the crystal. The Coulomb interaction between
the positively charged exchange-correlation hole and the electron that leaves the crystal
results in an attraction that scales like the inverse of the distance from the surface, giving
rise to the classical image force. Similarly, the attractive van der Waals interaction between
neutral objects may be framed in terms of exchange-correlation effects.
The exchange-correlation energy plays an important role in determining conceptually
and quantitatively many other properties of solids. We now give a simple example of how
the cohesive energy of the alkali metals can be estimated. The cohesive energy per atom
(which is the same as the cohesive energy per electron for the alkali metal case if we
consider only one active electron per alkali metal atom) is defined by
where Eatom , using the pseudopotential concept (discussed in Chapter 7), is the binding
energy of the outer active electrons in the atom for the isolated case, and EGS is the energy
per atom in the solid state. The cohesive energy Ecoh is positive for a stable solid.
139 6.7 The exchange-correlation energy
Table 6.1 Calculated cohesive energies of the alkali metals within the
interacting electron gas model for electron–electron interaction using
different approximations – Hartree, Hartree–Fock (HF), and Hartree–Fock
plus Wigner correlations (HF + Wigner). All energies are in kcal/mol.
HF+Wigner
EHartree
coh EHF
coh Ecoh expt.6
For the alkali metal elements, Eatom is just the negative of the ionization potential and
for the solid
We may use the interacting homogeneous electron gas as a model for evaluating Eq. (6.84).
The kinetic energy is
2.21 me
EKE = 2 , (6.85)
rs mb
where mb is the conduction band effective mass of the alkali metal, and Ee-ion is the bottom
of the conduction band. Both mb and Ee-ion may be obtained by a band structure calculation.
The electron–electron interaction term Ee-e may be approximated by using the results for
the homogeneous electron gas, since the charge density of the conduction electrons of a
simple metal is expected to be nearly uniform.
1.2 0.916 0.88
Ee-e = EHartree + Ex + Ecor = + − 2 + − (Ry). (6.86)
rs rs 7.8 + rs
For the correlation energy in the above expression, we have used the Wigner interpolation
eq
formula. If a crystal is formed, EGS is a minimum at the equilibrium lattice constant rs ,
and hence Ecoh may be determined by setting
∂EGS (rs )
= 0. (6.87)
∂rs
rs =reqs
In the case of very light elements, there is a contribution Ezero-point to the cohesive energy
due to the zero-point motion of the atoms confined to the crystal, which can be evaluated
6 Experimental values from H. Brooks, “Cohesive energy of alkali metals,” Phys. Rev. 91(1953), 1027.
140 6 Many-electron interactions: the homogeneous interacting electron gas and beyond
via knowledge of the phonon spectrum. Table 6.1 shows the results for the calculated Ecoh
of alkali metals with this simple model using Hartree, Hartree–Fock, and Hartree–Fock +
Wigner correlation models. We see that the alkali metals are unbounded within the Hartree-
only approximation. Both exchange and correlation effects are needed to give a reasonable
description of the cohesive properties of this simplest class of metals.
7 Density functional theory (DFT)
In this chapter, we turn to first-principles calculations of the properties and behavior of real
materials. Tremendous progress has been made in the past decades in employing ab initio
or first-principles theories and numerical computations to explain and predict the proper-
ties, and even the existence, of condensed matter systems. These first-principles studies
have been particularly successful in investigating materials and reduced-dimensional
systems with a moderate amount of electron correlation.
As discussed in previous chapters, the structure and properties of condensed matter are
basically dictated by the outer valence electrons of its constituent atoms. The mutual in-
teractions of these electrons and their interactions with the ions determine the electronic
structure of the system, which in turn determines many of the properties of the mate-
rial. Understanding material properties from first principles involves solving an interacting
quantum many-body problem with the Hamiltonian (see Chapter 2)
For simplicity, we have omitted spin and relativistic effects in our discussion here. Gen-
erally, exact solutions to this problem are impractical, and often undesirable, since the
many-electron wavefunctions would be so complicated that it would be difficult to ob-
tain a physical understanding from the exact solutions. In this chapter, we discuss one
particular formalism, density functional theory, for solving this many-electron problem
in an ab initio fashion for properties related to the electronic ground state of the system.
Nowadays, density functional theory is arguably the most popular approach for the first-
principles study of the ground-state properties of molecules and solids, with applications
to many disciplines ranging from physics to chemistry to the biological and engineering
sciences. This formalism provides a means to treat electron–electron interactions and trans-
forms the many-electron problem to one of a self-consistent-field one-particle problem for
ground-state properties, which is exact in principle.
For an interacting many-electron system, it is useful to distinguish the ground-state prop-
erties from the electron excited-state or spectroscopic properties. At low temperature,
properties of a crystal such as the structure, cohesive energies, structural and vibrational
properties, and structural phase stability are ground-state properties because they are de-
termined collectively by all the electrons in the ground state. Many of these properties can
be obtained by knowing the ground-state total energy of the system as a function of the
atomic coordinates. However, electronic excitation properties such as those measured in
photoemission, transport, and tunneling experiments involve creating an excited particle
141
142 7 Density functional theory (DFT)
II. Ev [ρ] is a minimum for the correct physical density, with ρ(r) satisfying the constraint
N = ρ(r)dr.
III. ρ(r) and hence Ev [ρ] can in principle be exactly obtained from the solution of an
associated one-electron problem with an effective potential veff (r).
We demonstrate statement I by considering the simple case of a paramagnetic system
with a non-degenerate ground state. This is achieved using the technique of proof by con-
tradiction. Proofs of the more general cases may be found in the literature. Consider a
1 P. Hohenberg and W. Kohn, “Inhomogeneous electron gas,” Phys. Rev. 136(1964), B864.
143 7.1 The ground state and density functional formalism
collection of N electrons enclosed in a large box moving under the influence of an external
potential v(r) and their mutual Coulomb repulsion. The Hamiltonian for the electrons is
p2 1 e2
i
H= + v(ri ) +
2m 2 |ri − rj | (7.2)
i i i=j
≡ T + V + U,
where we have defined T, V, and U as the operators for the total kinetic energy, the interac-
tion energy with the external potential, and the electron–electron interaction energy of the
system, respectively. In the ground state |, the electron density is
ρ(r) = |ρ̂(r)| , (7.3)
with
ρ̂(r) = δ(r − ri ), (7.4)
i
H | = E | (7.5)
and
H
= E
. (7.6)
By the minimal property of the ground state for the primed system,
E = |H | < |H | , (7.7)
but
|H | = |H + V − V| = E + [v (r) − v(r)]ρ(r)dr, (7.8)
Adding Eqs. (7.9) and (7.10), and making use of the assumption that ρ (r) = ρ(r), we
have the contradictory result that
E + E < E + E . (7.11)
This contradiction thus proves that v(r) is a unique functional of ρ(r); that is, given a
physical electron density, there is a unique mapping to an external potential. Now we may
define a universal functional that is valid for any number of particles and any external
potential,
This is because knowing ρ implies knowing v, which in turn means knowing the ex-
act ground state [ρ]. The ground-state energy of the system is then a functional of the
electron density and is given by
Ev [ρ] = v(r)ρ(r)dr + F[ρ]. (7.13)
Statement II may also be proven using Eq. (7.13) and a variational principle of quantum
mechanics: the energy of the system is a minimum at the ground-state wavefunction |
relative to an arbitrary variation of |, keeping the number of particles constant. If F[ρ]
were a known and sufficiently simple functional of ρ(r), the problem of determining the
ground-state energy and density would be quite straightforward. However, F[ρ] is not
explicitly known.
Statement III of the Kohn–Sham formulation2 of DFT states that there exists a set of one-
body equations that exactly determine ρ(r). In fact, as shown below, these equations are
formally similar to other self-consistent field equations such as the Hartree equations, but
they yield the exact density in principle. These equations are derived using the following
theoretical construction developed by Kohn and Sham. In the expression for the ground-
state energy (Eq. (7.13)), F[ρ] contains both the electron kinetic energy and the electron–
electron interaction energy, which in general are very difficult to evaluate and for which
there are no known explicit expressions in terms of the density ρ(r). However, we do
know the contributions to F[ρ] in some limiting cases. If exchange and correlation effects
are neglected, the interaction energy is just the classical Hartree term
e2 ρ(r)ρ(r )
EHartree = drdr . (7.14)
2 |r − r |
2 W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. Rev.
140(1965), A1133.