Chapter 7
Chapter 7
Chapter 7
Lattice vibrations
7.1
Introduction
Up to this point in the lecture, the crystal lattice was always assumed to be completely rigid,
i.e. atomic displacements away from the positions of a perfect lattice were not considered.
For this case, we have developed a formalism to compute the electronic ground state for an
arbitrary periodic atomic configuration and therewith the total energy of our system. Although
it is obvious that the assumption of a completely rigid lattice does not make a lot of sense (it
even violates the uncertainty principle), we obtained with it a very good description for a wide
range of solid state properties, among which were the electronic structure (band structure), the
nature and strength of the chemical bonding, or low temperature cohesive properties like the
stable crystal lattice. The reason for this success is that the internal energy U of the solid (or
similarly of a molecule) can be written as
U = E static + E vib
(7.1)
where E static is the electronic ground state energy with a fixed lattice (on which we have focused so far) and E vib is the additional energy due to lattice motion. In general (and we will
substantiate this below), E static E vib and thus many properties are correctly obtained, even
when simply neglecting E vib .
There is, however, also a multitude of solid state properties, for which the consideration of
lattice mobility is crucial. Mobility is in this respect already quite a big word, because for
many of these properties it is enough to take small vibrations of the ions around their (ideal)
equilibrium position into account. And although the amplitude of these vibrations is small,
the important aspect about them is that they can describe deformations whose wavelength is
comparable to the interatomic distances which is thus outside of the validity of continuum
elasticity theory. The vibrational modes of crystalline lattices are called phonons, and most
salient examples of solid state properties for which they are of paramount importance are:
Heat capacity: In metals with a rigid perfect lattice, only the free electron like conduction
electrons can take up small amounts of energy (of the order of thermal energies). The
specific heat capacity due to these conduction electrons was computed to vary linearly
with temperature for T 0 in the first exercise. In insulators with the same rigid perfect
lattice assumption, not even these degrees of freedom are available. Correspondingly,
the only excitations would be electronic excitations across the (huge) band gap Eg
kB T . Their (vanishingly small) probability will be exp(Eg /kB T ), and the specific heat
35
36
would show a similar scaling at low temperatures accordingly. In reality, one observes,
however, that for both insulators and metals the specific heat varies predominantly with
T 3 at low temperatures, in clear contradiction to the above cited rigid lattice dependencies.
Thermal expansion: Due to the negligibly small probability of electronic excitations
over the gap, there is nothing left in a rigid lattice insulator that could account for the
pronounced thermal expansions typically observed in experiments. These expansions are,
on the other hand, key to many engineering applications of materials sciences (different
pieces simply have to fit together over varying temperatures). Likewise, for structural
phase transitions upon heating and finally melting the static lattice approximation fails
with a vengeance.
Transport properties: In a perfectly periodic potential Bloch electrons suffer no collisions (being solutions to the many-body Hamiltonian). Assuming only point defects
as responsible for finite path lengths, one can not explain quantitatively the really measured, and obviously finite electric and thermal conductance of metals and semiconductors. Scattering of electrons, propagation of sound, heat transport, or optical absorption
are all fundamental properties of real devices (like transistors, LEDs, lasers). They are
among other things responsible for the electric resistance and loss of the latter, and can
only be understood when considering lattice vibrations.
7.2
(7.2)
(7.3)
which does not contain the electronic degrees of freedom explicitly anymore. The complete
interaction with the electrons (i.e. the chemical bonding) is instead contained implicitly in the
form of the PES V BO ({RI }).
For a perfect crystal, the minima of this PES correspond to the points of a Bravais lattice (or
better: to its basis) {RI }. And it was the situation with the atoms fixed exactly at {RI }, which
we discussed exclusively in the preceding chapters in the rigid lattice model. There is, however, nothing that would confine the validity of the Born-Oppenheimer approximation to just
the minima, or at least (recalling the discussion of the first chapter) the Born-Oppenheimer approximation will be equally good or bad for the minima and for the PES vicinity around them.
37
Correspondingly, understanding the properties of a crystal with the atoms at arbitrary positions
{RI } boils in the adiabatic approximation equally down to evaluating the PES at these positions. However, as long as the solid is not yet molten, we dont even need really arbitrary {RI },
but only ones that are very close to the (pronounced) minima corresponding to the ideal lattice.
This allows us to write
RI = RI + sI
|sI | a
with
(7.4)
For such small displacements around a minimum, the PES can be approximated by a parabola
and leads us to the harmonic approximation.
2
V (x) s +
V (x) s +
V (x) s3 + . . . ,
(7.5)
V (x) = V (x0 ) +
2
3
x
2 x
3! x
x0
x0
x0
with the displacement s = x x0 . The linear term vanishes, because x0 is an equilibrium geometry. For small displacements, the cubic (and higher) terms will be comparably small, leaving
in the harmonic approximation only
1 2
V (x) V (x0 ) +
(7.6)
V (x) s2 .
2 x2
x0
Introducing the force acting on the particle
F =
V (x)
x
(7.7)
2
V (x)
x2
x0
V (x) = cs ,
x
(7.8)
BO
(R1 , R2 , . . . , RM ) = V
BO
({RI })
1 M,M 3,3
2
BO
V ({RI })
+ sI, sJ,
2 I=1, =1,
RI, RJ,
{R }
J=1 =1
(7.9)
38
where , = x, y, z run over the three Cartesian coordinates, and the linear term has vanished
again close to the equilibrium geometry. Defining the (3 3) matrix
2
BO
,
(7.10)
V ({RI })
(RI , RJ ) =
RI, RJ,
{R }
I
the Hamilton operator for the M ions of the solid (and the N electrons implicitly contained in
V BO ) can then be written as
1 M,M 3,3
p2I
BO
0
+V
({R
})
+
sI,(RI , RJ )sJ,
I
2
I=1 =1
I=1 2Mion,I
M
H =
(7.11)
J=1 =1
(7.12)
This way, the second derivative can be obtained numerically, by simply displacing an atom
slightly from its equilibrium position and monitoring the arising forces on this and all other
atoms. From this force point of view, it also becomes clear that many elements of the second
derivative matrix will vanish: a small displacement of atom I will only lead to non-negligible
forces on the atoms in its immediate vicinity. Atoms further away will not notice the displacement, and hence the second derivative between these atoms J and the initially displaced atom I
will be zero.
(7.13)
(RI , RJ ) = 0 for |RI RJ | a0 ,
where a0 is the lattice parameter. Very often, one even finds
(RI , RJ ) 0 for
(7.14)
Despite these simplifications, this numerical procedure is at this stage still unfortunately intractable, because we would have to do it for each of the M atoms in the solid (i.e. a total of 3M
times). With M 1023 , this is unfeasible, and we need to exploit further properties of solids to
reduce the problem. Obviously, just like in the case of the electronic structure problem a much
more efficient use of the periodicity of the lattice must be made. This is, however, not done
at the level of further specifying the total Hamilton operator, but in the equations of motion
following from the general form of eq. (7.11).
39
In the classical case and knowing the structure of the Hamilton function from eq. (7.11), the
total energy (as sum of kinetic and potential energy) can be written as
Mion,I .
1 M,M 3,3
BO
s
+V
({R
})
+
I
=1,
sI,(RI , RJ )sJ,
I
2 I=1,
I=1 2
M
E =
(7.15)
J=1 =1
The classical equation of motion of the nuclei in our system (within the harmonic approximation) can be written as
..
(7.16)
Mion,I sI, = (RI , RJ )sJ, ,
J,
where the right hand side represents the -component of the force vector acting on atom I of
mass Mion,I and the double dot above s represents the second time derivative. As a first step, we
will get rid of the time dependence. A general form of solution for the second-order differential
equation (eq. (7.16)) is
sI, (t) = uI, eit .
(7.17)
Inserting into eq. (7.16) leads to
Mion,I 2 uI, =
(RI , RJ )uJ,
(7.18)
J,
This is a system of coupled algebraic equations with 3M unknowns (still with M 1023 !). As
already remarked above, it is hence not directly tractable in this form, and we exploit (just like
in the electronic case), that the solid has in its equilibrium geometry a periodic lattice structure.
Correspondingly, we decompose the coordinate RI into a Bravais vector Rn (pointing to the unit
cell) and a vector R pointing to atom of the basis within the unit cell,
RI = Rn + R
(7.19)
Formally, we can therefore replace the index I over all atoms in the solid by the index pair
(n, ) running over the unit cells and different atoms in the basis. With this, the matrix of
second derivatives takes the form
(RI , RJ ) = (n, n )
(7.20)
And in this form, we can now start to exploit the periodicity of the Bravais lattice. From the
translational invariance of the lattice follows e.g. immediately that
(n, n ) = (n n )
(7.21)
From a classical point of view, {sI (t)} is simply a snap shot of the displacement field of the
atoms at time t. This suggests, that uI may be describable by means of a vibration
1
u (Rn ) = p
c eikRn
Mion,
(7.22)
(7.23)
(n n )eik(Rn Rn ) c .
2 c = p
M
M
40
Note that the part in curly brackets is only formally still dependent on the position RI (through
the index n). Since all unit cells are equivalent, RI is just an arbitrary zero for the summation
over the whole Bravais lattice. This part in curly brackets
1
D (k) = p
(n n )eik(Rn Rn )
M M n
(7.24)
is called dynamic matrix, and its definition allows us to rewrite the homogeneous linear system
of equations in a very compact way
2 c =
D(k)c
(7.25)
Compared to the original equation of motion, cf. eq. (7.18), a dramatic simplification has been
reached: due to the lattice symmetry we achieved to reduce the system of 3r M equations to 3r
equations, where r is the number of atoms in the basis of the unit cell. The price we paid for it is
that we have to solve the equations of motion for each k anew. This is comparably easy, though,
since this scales linearly with the number of k and not cubic like the matrix inversion, and has
to be done only for few k (since we will see that the k-dependent functions are quite smooth).
From eq. (7.25) one gets then as condition for the existence of solutions the eigenvalue problem
(i.e. solutions must satisfy the following:)
det D(k)
2 1 = 0 .
(7.26)
This problem has 3r eigenvalues, which are functions of the wave vector k,
= i (k)
i = 1, 2, . . . , 3r
(7.27)
and to each of these eigenvalues i (k), eq. (7.25) provides the eigenvector
(i)
c = e (k)
(i)
viz. c = e (k)
(7.28)
These eigenvectors are determined up to a constant factor, and this factor is finally obtained by
requiring all eigenvectors to be normalized (and orthogonal to each other).
In conclusion, one then determines for the displacement sn (t) of atom I = (n, ) away from
the equilibrium position at time t
1
(i)
(i)
sn (t) = p
e (k)ei(kRn i (k)t)
Mion,
(7.29)
as a special solution to the equation of motion (7.16), i.e. the one corresponding to the particular normal mode (k, (k)). With this, the system of 3M three-dimensional coupled oscillators
is transformed into 3M decoupled (but collective) oscillators, and the solutions (or so-called
normal modes) correspond to waves, spreading over the whole crystal (i.e. these waves have
only a formal resemblance to what we think of traditionally as an oscillator). Any arbitrary displacement situation of the atoms (i.e. the general solution) can now be written as an expansion
in these special solutions, i.e. the waves form the basis for the description of lattice vibrations.
41
(7.30)
i.e. the sum over all atoms breaks down to left and right nearest neighbor, and the second
derivative of the PES corresponds simply to the spring constant. We solve this problem as
42
unit cell
a
M
s(n-1)
s(n)
s(n+1)
s(n+2)
Figure 7.1: Linear chain with one atomic species (Mass M, and lattice constant a).
(k) / (f/M)1/2
2.0
1.5
1.0
0.5
0.0
/a
+/a
Figure 7.2: Dispersion relation (k) for the linear chain with one atomic species.
2 M = f 2 eika eika
(7.31)
(7.32)
(7.33)
(k) = 2
f
M
sin ka
2
i.e. as expected is a periodic function of k and symmetric with respect to k and k. The first
period lies between k = /a and k = +/a (first Brillouin zone), and the form of the dispersion
is plotted in Fig. 7.2. Equation (7.31) gives the explicit solutions for the displacement pattern
in the linear chain for any wave vector k. They describe waves propagating along the chain with
phase velocity /k and group velocity vg = /k.
s(n-1) = -s(n)
s(n)
43
s(n+4) = s(n)
Figure 7.3: Snap shot of the displacement pattern in the linear chain at k = +/a.
Let us look in more detail at the two limiting cases k 0 and k /a (i.e. the middle and
border of the Brillouin zone). For k 0 we can approximate (k) given by eq. (7.33) with its
leading term in the small k Taylor expansion (sin(x) x for x 0)
r !
f
k .
(7.34)
a
M
p
is thus linear in k, and the proportionality constant is simply the group velocity vg = a f /M,
which is in turn identical to the phase velocity and independent of frequency. For the displacement pattern we obtain in this limit
g
1
(7.35)
sn = c eik(anv t) ,
M
i.e. at these small k-values (long wavelength) the phonons spread simply as vibrations with
constant speed vg (the ordinary sound speed in this 1D crystal!). This is the same result we
would have equally obtained within linear elasticity theory, which is understandable because for
such very long wavelength displacements the atomic structure and spacing is less important (the
atoms move almost identically on a short length scale). Fundamentally, the existence of a branch
of solutions whose frequency vanishes as k vanishes results from the symmetry requiring the
energy of the crystal to remain unchanged, when all atoms are displaced by an identical amount.
We will therefore always obtain such modes, regardless of whether we consider more complex
interactions (e.g. spring constants to second nearest neighbors) or higher dimensions. Since
their dispersion is close to the center of the Brillouin zone linear in k, which is characteristic of
sound waves, these modes are commonly called acoustic branches.
One of the characteristic features of waves in discrete media, however, is that such a linear
behaviour ceases to hold at wavelengths short enough to be comparable to the interparticle
spacing. In the present case falls below vg k as k increases, the dispersion curve bends over
and becomes flat (i.e. the group velocity drops to zero). At the border of the Brillouin zone (for
k = /a) we then obtain
1
sn = cein eit .
(7.36)
M
Note that ein = 1 for even n, and ein = 1 for odd n. This means that neighboring atoms
vibrate against each other as shown in Fig. 7.3. This also explains, why in this case the highest
frequency occurs. Again, the flattening of the dispersion relation is required by pure symmetry:
From the repetition in the next Brillouin zone and with kinks forbidden, /k = 0 must result
at the Brillouin zone boundary, regardless of how much more complex and real we make our
model.
44
a
-a/4
+a/4
M1
s(1,n)
s(2,n)
s(1,n+1)
M2
s(2,n+1)
Figure 7.4: Linear chain with two atomic species (Masses M1 and M2). The lattice constant is
a and the interatomic distance is a/2.
sn
(2)
sn
This leads to
2
c1
2 c2
1
1
= c1 ei(k(n 4 )at )
M1
1
1
= c2 ei(k(n+ 4 )at )
M2
(7.39)
.
2f
2f
ka
=
c1 +
c2 cos
2
M1 M2
M1 M2
2f
2f
ka
=
c2 +
c1 cos
2
M1 M2
M1 M2
(7.40)
(7.41)
,
(7.42)
i.e. in this case the frequencies of the two atoms are still coupled to each other. The eigenvalues
follow then from the determinant
2f
2f
ka
2
cos
2
M1 M2
M 1 M2
= 0 ,
(7.43)
2 f cos ka
2
f
2
2
M1 M2
M1 M2
which yields the two solutions
2
1
ka
1
1 2
4
1
+
f
+
sin2
= f
M1 M2
M1 M2
M1 M2
2
s
4
f
1
ka
f
=
sin2
,
e
e 2 M1 M2
2
M
M
(7.44)
45
(k) / (f/M)1/2
2.0
1.5
1.0
0.5
0.0
/a
+/a
Figure 7.5: Dispersion relations (k) for the linear chain with two atomic species (the specific
values are for M1 = 3M2 = M). Note the appearance of the optical branch not present in Fig.
7.2.
1
e= 1 + 1
. We therefore obtain two branches for (k), the
with the reduced mass M
M1
M2
dispersion of which is shown in Fig. 7.5. As apparent, one branch goes to zero for k 0,
(0) = 0, as in the case of the monoatomic chain, i.e. we have again an acoustic branch.
q The
other branch, on the other hand, exhibits interestingly a (high) finite value + (0) = 2 ef for
M
k 0 and shows a much weaker dispersion over the Brillouin zone. Whereas in the acoustic
branch both atoms within the unit cell move in concert, they vibrate against one another in opposite directions in the other solution. In ionic crystals (i.e. where there are opposite charges
connected to the two species), these modes can therefore couple to electromagnetic radiation
and are responsible for much of the characteristic optical behavior of such crystals. Correspondingly, they are generally (i.e. not only in ionic crystals) called optical modes, and their
dispersion is called the optical branch.
q
2f
At the border of the Brillouin zone at k = /a we finally obtain the values + = M
and
1
q
2f
, i.e. as apparent from Fig. 7.5 a gap has opened up between the acoustic and the
= M
2
optical branch. The width of this gap depends only on the mass difference of the two atoms
(M1 vs. M2 ). For identical mass of the two species it closes, which is comprehensible, since
then our model with equal spring constants and equal masses gets identical to the previously
discussed model with one atomic species: the optical branch is then nothing but the upper half
of the acoustic branch, folded back into the half-sized Brillouin zone.
46
Figure 7.6: DFT-LDA and experimental phonon band structure (left) and corresponding density
of states (right) for fcc bulk Rh [from R. Heid, K.-P. Bohnen and W. Reichardt, Physica B 263,
432 (1999)].
dynamic matrix, and the different branches of the band structure follow from the eigenvalues
after diagonalizing this matrix. This is also, what one computes in a quantitative approach to the
phonon band structure problem, e.g. with density-functional theory. Since the dynamic matrix
derives from the second derivatives of the total energy, one exploits for its computation in the
so-called direct (supercell) approach (or frozen-phonon approach) the relation with the forces
expressed in eq. (7.12). In other words, the different atoms in the unit cell are displaced by
small amounts and the forces on all other atoms are recorded. In order to be able to compute the
dynamic matrix with this procedure at different points in the Brillouin zone, larger supercells
are employed and the (in principle equivalent) atoms within this larger periodicity displaced differently, so as to account for the different phase behavior. Obviously, this way only k-points that
are commensurate to the used supercell periodicity can be accessed, and one runs into the problem that quite large supercells need to be computed in order to get the full dispersion curves.
Approaches based on linear response (i.e. by inverting the dielectric matrix) are therefore also
frequently employed. The details of the different methods go beyond the scope of the present
lecture, but a good overview can for example be found in R.M. Martin, Electronic Structure:
Basic Theory and Practical Methods, Cambridge University Press, Cambridge (2004). Experimentally, phonon dispersion curves are primarily measured using inelastic neutron scattering,
and Ashcroft/Mermin dedicates an entire chapter to the particularities of such measurements.
Overall, the measured and DFT-LDA/GGA calculated phonon band structures agree extremely
well these days (regardless of whether computed by a direct or a linear response method), and
instead of discussing the details of either measurement or calculation we rather proceed to some
practical examples.
If one has a mono-atomic basis in the three-dimensional unit cell, we expect from the onedimensional model studies that only acoustic branches should be present in the phonon band
structure. There are, however, now 3 such branches, because in three-dimensions the orientation of the polarization vector (i.e. the eigenvector belonging to the eigenvalue (k)), cf. eqs.
(7.28) and (7.29) matters. In the linear chain, the displacement was only considered along the
chain axis and such modes are called longitudinal. If the atoms vibrate in the three-dimensional
case perpendicular to the propagation direction of the wave, two further transverse modes are
obtained, and longitudinal and transverse modes must not necessarily exhibit the same disper-
47
Figure 7.7: DFT-LDA and experimental phonon band structure (left) and corresponding density
of states (right) for Si (upper panel) and AlSb (lower panel) in the diamond lattice. Note the
opening up of the gap between acoustic and optical modes in the case of AlSb as a consequence
of the large mass difference. 8 cm1 1 meV [from P. Giannozzi et al., Phys. Rev. B 43, 7231
(1991)].
sion. Fig. 7.6 shows the phonon band structure for bulk Rh in the fcc structure, as well as
the phonon density of states (DOS) resulting from it (which is obtained analogously to the
electronic structure case by integrating over the Brillouin zone). The dispersion of the bands
particularly along the lines X and L, which run both from the Brillouin zone center to
the border, has qualitatively the form discussed already for the linear chain model and is thus
conceptually comprehensible.
Turning to a poly-atomic basis, the anticipated new feature is the emergence of optical branches
due to vibrations of the different atoms against each other. This is indeed apparent in the phonon
band structure of Si in the diamond lattice shown in Fig. 7.7, and we find for these branches
again the rough shape that also emerged out of the di-atomic linear chain model (check again
the X and L curves). Since there are in general 3r different branches for a basis of r
atoms, cf. section 7.2.4, 3r 3 optical branches result (i.e. in the present case of a diatomic
basis, there are 3 optical branches, one longitudinal one and two transverse ones). Although the
diamond lattice has two inequivalent lattice sites (and represents thus a case with a di-atomic
basis), these sites are both occupied by the same atom type in bulk Si. From the discussion of
the di-atomic linear chain, we would therefore expect a vanishing gap between the upper end
of the acoustic branches and the lower end of the optical ones. Comparing with Fig. 7.7 this is
indeed the case. Looking at the phonon band structure of AlSb also shown in this figure, a huge
gap can now, however, be discerned. This is obviously the result of the large mass difference
between Al and Sb in this compound, and shows us that we can indeed understand most of the
qualitative features of a real phonon band structure from the analysis of the two simple linear
chain models. One should also have this inverse variation of the phonon frequencies with mass,
48
cf. eq. (7.24), in mind, in order to develop a rough feeling for the order of magnitude of lattice
vibrational frequencies (or in turn the energies contained in these vibrations). For transition
metals, cf. Fig. 7.6, we are talking about a range of 20-30 meV for the optical modes, while
for lighter elements, cf. Fig. 7.7 this becomes increasingly higher. First row elements like
oxygen in solids (thus oxides) exhibit optical vibrational modes around 70-80 meV, whereas
the acoustic modes go in all cases down to zero at the center of the Brillouin zone. We are thus
with all frequencies in the range of thermal energies, which already tells us how important these
modes will be for the thermal behavior of the solid.
7.3
In the preceding section we have learnt to compute the lattice vibrational frequencies, when the
ions are treated as classical particles. In practice, this approximation is not too bad for quite a
range of materials properties. This has not only to do with that the mass of the ions is (with the
exception of H and He) quite heavy, but also with that the quantum mechanical and the classical
harmonic oscillator (on which all of the harmonic theory of lattice vibrations is based) exhibit
many similar properties. In a very rough view, one can say that if only the vibrational frequency
(k) of a mode itself is relevant, the classical treatment gives already a good answer (which is
ultimately why we could understand the experimentally measured phonon band structure with
our classical analysis). If, on the other hand, the energy connected with this mode (k) is
important, or more precisely, the way how one excites (or populates) this mode is concerned,
then the classical picture will fail. The prototypical example for this is the specific heat due to
the lattice vibrations, which is nothing else but the measure of how effectively the vibrational
modes of the lattice can take up thermal energy. In order to understand material properties of
this kind, we need a quantum theory of the harmonic crystal. Just as in the classical case, let
us first recall this theory for a one-dimensional harmonic oscillator before we really address the
problem of the full solid state Hamiltonian.
(7.46)
(7.47)
49
which is in this case nothing but the Heisenberg uncertainty principle. The resulting Hamiltonian has then still the same form as in eq. (7.45), but p and q have now to be treated as
(non-commuting) operators.
This structure can be considerably simplified by introducing the lowering operator
r
r
m
1
q+i
p ,
(7.48)
a =
2
h
2
hm
and the raising operator
a =
r
m
1
qi
p .
2
h
2
hm
(7.49)
Note that the idea behind this pair of adjoint operators leads to the concept called second quantization, where a and a are denoted as annihilation and creation operator, respectively. The
canonical commutation relation of eq. (7.47) implies then that
[a, a ] = 1 ,
(7.50)
1
= h
a a+
2
(7.51)
From this commutation relation and the form of the Hamiltonian it follows that the eigenstates
|n> (n = 0, 1, 2, . . .) of this problem fulfill the relations (see any decent textbook on quantum
mechanics for a derivation)
a |n> = (n + 1)1/2 |n + 1>
1/2
a|n> = n |n 1>
a|0> = 0 .
(7.52)
(n 6= 0)
(7.53)
(7.54)
Application of the operator a brings the system therefore into one excited state higher, and
application of the operator a one state lower. This explains the names given to these operators,
in particular when considering that each excitation can be seen as one quantum. The creation
operator creates therefore one quantum, while the annihilation operator wipes one out. With
these recursion relations between the different states, the eigenvalues (i.e. the energy levels) of
the one-dimensional harmonic oscillator are finally obtained as
QM
|n> = (n + 1/2) h
En = < n|H1D
(7.55)
Note, that E0 6= 0, i.e. even in the ground state the oscillator has some energy (so-called energy
due to zero-point vibrations or zero-point energy), whereas with each higher excited state one
quantum h
is added to the energy. With the above mentioned picture of creation and annihilation operators, one therefore commonly talks about the state En of the system as corresponding
to having excited n phonons of energy h.
50
sI,(RI , RJ )sJ,
I
2
I=1 I=1
I=1 2Mion
M
H =
(7.56)
J=1 J=1
Inserting the displacement field for the normal modes sI = s(n,) given in eq. (7.29), this can be
simplified to the form
H
vib
p2(i) (k)
2Mion
i=1 k
1 3
(i) 2
2
+ ki (k) s (k)
2 i=1
(7.57)
where we have focused now only on the vibrational part of the ionic Hamiltonian (denoted by
H vib ), i.e. we discard the contribution from the static minimum of the potential energy surface
(V BO ({RI }) = 0). With this transformation the Hamiltonian has apparently decoupled to a
sum over 3M harmonic oscillators (compare with eq. (7.45)!), and it is thus straightforward
to generalize the solution discussed in detail in the last section to this more complicated case.
We define annihilation and creation operators for a phonon (i.e. normal mode) with frequency
i (k) and wave vector k as
s
r
Mion i (k) (i)
1
ai (k) =
s (k) + i
p (k)
(7.58)
2
h
2
hMion i (k) (i)
s
r
Mion i (k) (i)
1
ai (k) =
s (k) i
p (k) .
(7.59)
2
h
2
hMion i (k) (i)
Inverting the last two expressions gives the displacement s(i) (k) and momentum p(i) (k) as a
function of annihiliation and creation operators. If this is inserted into eq. (7.57) one can show
(in a lengthy but straightforward calculus) that the Hamiltonian is brought to a form equivalent
to eq. (7.51), i.e.
3
1
vib
H
= h
i (k) ai (k)ai (k) +
.
(7.60)
2
i=1 k
This also implies that also the eigenvalues must have an equivalent form, namely
E
vib
i (k)
h
i=1 k
1
ni (k) +
2
(7.61)
The energy due to the lattice vibrations comes therefore from 3M harmonic oscillators with frequency (k). This frequency is the same as in the classic case, which is (as already remarked
in the beginning of this section) why the classical analysis permitted us already to determine
51
the phonon band structure. In the harmonic approximation the 3M modes are completely independent of each other: Each such mode can be excited independently and contributes then the
energy h(k).
While this sounds very familiar to the concept of the independent electron gas
(with its energy contribution n (k) from each state (n, k)), the fundamental difference is that in
the electronic system each state (or mode in the present language) can only be occupied with
one electron (or two in the case of non-spin polarized systems). This results from the Pauli principle and we talk about fermionic particles. In the phonon gas, on the other hand, each mode
can be excited with an arbitrary number of (indistinguishable) phonons, i.e. ni (k) is not restricted to the values 0 and 1. Phonons are therefore bosonic particles. The number of phonons
is furthermore not constant, but is dictated by the energy of the lattice and therewith through
the temperature. At T = 0 K there are no phonons (ni (k) 0), and the lattice is static so to
speak. Note, however, that even then the lattice energy is not zero, cf. eq. (7.61), but receives
a contribution from the quantum mechanical zero-point vibrations. At increasing temperatures,
more and more phonons are created, the lattice starts to vibrate and the lattice energy rises.
kEnT
B
l e
El
BT
l e
nhi (k)
kB T
lhi (k)
kB T
xn
l xl
(7.62)
i (k)
where we have defined the short hand variable x = exp(h
kB T ) in the last step. From mathematics, the value of the geometric series is known as
xn
n=0
1
1x
(7.63)
(7.64)
With this probability the mean energy of the oscillator at temperature T becomes
E i (T, (k)) =
n=0
n=1
nPn (T )
EnPn(T ) = E0 +h
=
h
i (k)
+h
i (k) (1 x)
2
nxn
(7.65)
n=0
nxn
n=0
x
(1 x)2
(7.66)
52
which finally leads to
E i (T, (k)) = h
i (k)
e
h
i (k)
kB T
1
+
2
(7.67)
as the mean energy of this oscillating mode. Comparing with the energy levels of the mode
En = h
i (k)(n + 1/2), we can identify
1
ni (k) =
e
h
i (k)
kB T
(7.68)
1
as the mean excitation (or occupation with phonons) of this particular mode at T . This is simply
the famous Bose-Einstein statistics, consistent with the above made remark that phonons behave
as bosonic particles.
Knowing the energy of one mode at T , cf. eq. (7.67), it is easy to generalize over all modes and
arrive at the total energy due to all lattice vibrations at temperature T
E vib (T ) =
i (k)
h
i=1 k
1
e
h
i (k)
kB T
1
+
2
(7.69)
This is still a formidable expression and we notice that the mean occupation of mode i with
frequency i (k) and wave vector k does not explicitly depend on k itself (only indirectly via
the dependence of the frequency on k). It is therefore convenient to reduce the summations over
modes and wave vectors to one simple frequency integral, by introducing a density of normal
modes (or phonon density of states) g()d, defined so that g() is the total number of modes
with frequencies in the infinitesimal range between and + d per volume,
1
g()d =
(2)3
i
This implies obviously, that
dk ( i (k)) .
d g() = 3M/V
(7.70)
(7.71)
i.e. the integral over all modes yields the total number of modes per volume, which is three
times the number of atoms per volume, cf. section 7.2.4.
With this definition of g(), the total energy due to all lattice vibrations at temperature T of a
solid of M atoms in a volume V takes the simple form
E vib (T ) = V
Z
0
d g()
"
#
1
+ h
,
2
1
1
e
kB T
(7.72)
i.e. the complete information about the specific vibrational properties of the solid is contained in
the phonon density of states g(). Once this quantity is known, either from experiment or DFT
(see the examples in Figs. 7.6 and 7.7), the contribution from the lattice vibrations can be evaluated immediately. Via thermodynamic relations, also the entropy due to the vibrations is then
accessible, allowing to compute the complete contribution from the phonons to thermodynamic
quantities like the free energy or Gibbs free energy.
53
(7.73)
i.e. it measures (normalized per volume) the amount of energy per temperature taken up by a
solid. If we want to specifically evaluate the contribution from the lattice vibrations, we arrive
with the help of eq. (7.72) at the general expression
cVvib (T )
#
"
Z
1
1
1 E vib
+ h
d g() h
=
=
V T V
T 0
2
kB T 1
e
#
"
Z
1
h
d g() h
T 0
e kB T 1
(7.74)
where in the last step we have dropped the temperature-independent zero-point vibrational term.
If one plugs in the phonon densities of states from experiment or DFT, the calculated cVvib agree
excellently with the measured specific heats of solids, proving that indeed the overwhelming
contribution to this quantity comes from the lattice vibrational modes. The curves have for
most elements roughly the form shown in Fig. 7.8. In particular, cV approaches a constant
value at high temperatures, while in the lowest temperature range it scales with T 3 (in metals,
this scaling is in fact T +T 3 , because there is in addition the contribution from the conduction
electrons, which you calculated in the first exercise). But, let us again try to see, if we can
understand these generic features on the basis of the theory developed in this chapter (and not
simply accept the data coming out from DFT or experiment. . .).
54
Figure 7.8: Qualitative form of the specific heat of solids as function of temperature (here
normalized to a temperature D below which quantum effects become important). At high
temperatures, cV approaches the constant value predicted by Dulong-Petit, while it drops to
zero at low temperatures.
cVvib (T ) =
d g() h
h
T 0
kB T 1
e
"
#
Z
1
d g()
h
T 0
(1 + kh
+
.
.
.)
1
T
B
Z
(kB T )
T
0
= 3kB M/V = constant ,
=
d g()
(7.75)
where in the last step we exploited eq. (7.71) for the total number of phononic states. We
therefore indeed find the specific heat to become constant at high temperatures, and the absolute
value of cVvib is entirely given by the density of the solid (M/V ). This was long known as the
Dulong-Petit law from classical statistical mechanics. That we recover it in this limit is simply
a consequence of the fact that at such high temperatures the quantized nature of the vibrations
does not show up anymore. Energy can be taken up in a quasi-continuous manner and we arrive
at the classical result that each degree of freedom yields kB T to the energy of the solid.
Intermediate temperature range (Einstein approximation)
The constant cVvib obtained in the high-temperature limit was not a surprising result, since it
coincides with what we expect from classical statistical mechanics. This classical understanding can, however, not at all grasp, why the specific heat deviates at intermediate temperatures
(k)
55
/a
+/a
Figure 7.9: Left: Illustration of the dispersion curves behind the Einstein and Debye model
(dashed curves). Einstein approximates optical branches, while Debye focuses on the acoustic
branches relevant at low temperatures. Right: Specific heat curves obtained with both models.
from the constant value predicted by the Dulong-Petit law. This lowering of cVvib can only be
understood on the basis of a quantum theory, namely the discrete nature of the excitation levels
available for lattice vibrations. The simplest model that would go beyond a classical theory
and consider this quantized nature of the oscillations is historically due to Einstein. Ignoring to
zeroth order any interaction between the vibrating atoms in the solid, it is reasonable to assume
that all atoms of one species vibrate with exactly one characteristic frequency (independent harmonic oscillators). This is also obvious, because without interaction, the dispersion curves of
the various phonon branches become simply constant. For a mono-atomic solid, the phonon
density of states breaks therefore in the Einstein model down to
g()Einstein = 3(M/V ) ( E )
(7.76)
where E is this characteristic frequency (often called Einstein frequency), with which all atoms
vibrate. Inserting this form of g() into eq. (7.74) leads to
2 x
x e
Einstein
cV
(T ) = 3kB M/V
,
(7.77)
x
(e 1)2
with the short hand variable x = h
E /(kB T ), which expresses the relation between vibrational
and thermal energy. And from eq. (7.77) we see already the fundamentally important result: While in the high-temperature limit (x 1) the term in brackets approaches unity and
we recover the Dulong-Petit law, it is exactly this term that also declines for intermediate temperatures and finally drops to zero for x 1 (i.e. T 0). The resulting specific heat curve
shown in Fig. 7.9 exhibits therefore all qualitative features of the experimental data, which
is a big success for such a simple model (and was a sensation hundred years ago!). To be
more specific, the cVEinstein (T ) curve starts to deviate significantly from the Dulong-Petit value
for x 1, i.e. for temperatures E h
E /kB . E is known as the Einstein temperature, and
for temperatures above it, all modes can easily be excited and a classical behavior is obtained.
For temperatures of the order or below E , however, the thermal energy is no longer enough
to quasi-continuously excite the corresponding phonon mode, the discrete nature of the levels
starts to show up, and one says that this mode starts to freeze-out. In other words, any phonon
56
57
model recovers therefore the correct low temperature limit. It fails on the other hand for the
intermediate temperature range. As illustrated in Fig. 7.9 the Debye model deviates in this
range substantially from the Einstein model, which was found to reproduce the experimental
data very well. The reason is obviously that at these temperatures the optical branches get
noticeably excited, the dispersion of which is badly approximated by the Debye model, cf. Fig.
7.9.
The specific heat data from a real solid can therefore be well understood from the discussion
of the three ranges just discussed. At low temperatures, only the acoustic modes contribute
noticeably to cV and the Debye model provides a good description. With increasing temperatures, optical modes start to become excited, too, and the cV curve changes gradually to the
form predicted by the Einstein model. Finally, at highest temperatures all modes contribute
quasi-continuously and the classical Dulong-Petit value is approached. Similar to the role of
E in the Einstein model, also the Debye temperature D indicates the temperature range above
which classical behavior sets in and below which modes begin to freeze-out due to the quantized nature of the lattice vibrations. Debye temperatures are normally obtained by fitting the
predicted T 3 form of the specific heat to low temperature experimental data. They are tabulated
for all elements and are typically of the order of a few hundred degree Kelvin, which is thus the
natural temperature scale for phonons.
Note that the Debye temperature plays therefore the same role in the theory of lattice vibrations
as the Fermi temperature plays in the theory of metals, separating the low-temperature region
where quantum statistics must be used from the high-temperature region where classical statistical mechanics is valid. For conduction electrons, TF is of the order of 20000 K, though, i.e.
at actual temperatures we are always well below TF in the electronic case, while for phonons
both classical and quantum regimes can be encountered. This large difference between TF and
D is also responsible, why the contribution from the conduction electrons to the specific heat
of metals is only noticeable at lowest temperatures. In this regime you had derived in the first
exercise that cVel scales linearly with (T /TF ). On the other hand, the phonon contribution scales
with (T /D )3 . If we therefore equate the two, we see that the conduction electron contribution
will show up only for temperatures
3D
TF
(7.80)
58
7.4
So far our general theory of lattice vibrations has been based on the very initial assumption that
the amplitudes of the oscillations are quite small. This has led us to the so-called harmonic approximation, where higher order terms beyond the quadratic one in the Taylor expansion around
the PES minimum are neglected. This approximation seems plausible, keeps the mathematics
tractable and can indeed explain a wide range of materials properties. We had exemplified this
in the last section for a fundamental equilibrium property like the specific heat. Similarly, we
would find that lattice vibrations are central in explaining transport properties like the electric
conductance. More precisely it is the scattering of the Bloch electrons as solutions of a perfect lattice at the irregularities introduced by a slightly vibrating lattice, which is important for
such transport properties (which is shortly denoted as scattering at phonons in the quasi-particle
picture). There are, however, also a number of materials properties, which can not be understood within a harmonic theory of lattice vibrations. The most important of these are thermal
expansion and heat transport. Lets consider them in turn:
1
=
3B
as
V
P
P
T
(7.82)
(7.83)
U
P=
V
(7.84)
59
hi,k
ni,k
1
=
3B
V
T V
T
i,k
(7.85)
As 7.85 contains terms common to CV it is common to express the thermal expansion coefficient
in terms of the specific heat:
CV
=
,
(7.86)
3B
where is the Gruneisen parameter. Clearly the Gruneisen parameter is not a simple function as
it embodies the variation of all normal modes in the crystal with changes of volume. Normally
one defines the Gruneisen parameter for individual normal modes and it is given by
i,k =
(ln i,k )
(lnV )
(7.87)
The overall Gruneisen parameter is then determined as a weighted average in which the contribution of each mode is weighted by its contribution to CV ,
i,k i,kCV
i,k
i,k
i,k CV
(7.88)
Often turns out to be 1 (e.g. NaCl = 1.6, KBr = 1.5), which indicates that will exhibit a
similar temperature dependence as CV . Specifically becomes constant at high T and at low T
approaches zero as T 3 .
60
approach and what is due to anharmonic effects on a conceptual level, the actual computation
of thermal expansion coefficients or thermal conductivities becomes quickly rather involved,
without yielding a lot of fundamental new insight into the physics of solids. We will therefore
not dwell longer on this aspect in this lecture, but refer to the literature for further reading.
Corresponding treatments can e.g. be found in a complete chapter dedicated to this subject in
Ashcroft and Mermin.