International
Tables for
Crystallography
Volume D
Physical properties of crystals
Edited by A. Authier

International Tables for Crystallography (2006). Vol. D, ch. 2.1, pp. 266-274

## Section 2.1.2. Fundamentals of lattice dynamics in the harmonic approximation

G. Eckolda*

aInstitut für Physikalische Chemie, Universität Göttingen, Tammannstrasse 6, D-37077 Göttingen, Germany
Correspondence e-mail: geckold@gwdg.de

### 2.1.2. Fundamentals of lattice dynamics in the harmonic approximation

| top | pdf |

#### 2.1.2.1. Hamiltonian and equations of motion

| top | pdf |

In order to reduce the complexity of lattice dynamical considerations, we describe the crystal's periodicity by the smallest unit needed to generate the whole (infinite) lattice by translation, i.e. the primitive cell. Each individual primitive cell may be characterized by a running index l and a vector pointing to its origin. Let there be N atoms per cell, the equilibrium positions of which are given by being the vector of the κth atom with respect to the origin of the primitive cell (see Fig. 2.1.2.1).

 Figure 2.1.2.1 | top | pdf |Definition of position vectors.

The set of vectors describes the structure of the perfect lattice. At a particular time t, however, the κth atom within the lth primitive cell, denoted by (κl), may be found at a position which differs slightly from the equilibrium position, the time-dependent displacement being The potential energy V of the whole crystal depends on the position vectors of all atoms, and is minimal if all atoms occupy their equilibrium positions. For small displacements, it can be expanded in a Taylor series with respect to :where denotes the Cartesian coordinate of in direction α. In the harmonic approximation, third and higher-order terms are neglected. In order to simplify the formulae, we now drop the time argument, keeping in mind that we are always dealing with dynamical displacements. The expansion coefficients in equation (2.1.2.4) are the partial derivatives of the potential energy with respect to the atomic displacements taken at the equilibrium positions:Using the matrix notationand dropping the constant , equation (2.1.2.4) readsThe product is just the force acting upon atom if the atom is displaced by (Fig. 2.1.2.2). Hence, the matrix may be regarded as a force constant matrix and its elements as force constants. These parameters may be calculated with the help of specific interaction models such as pair potentials, tensor-force models or more complicated many-body interactions.

 Figure 2.1.2.2 | top | pdf |Definition of the force acting on atom () when atom () is displaced by .

The Hamiltonian of the perfect harmonic crystal can now be written in the formif and are the momentum and the mass of atom (κl), respectively. Consequently, the equation of motion for a particular atom (κl) is given by Solutions of this set of coupled differential equations are of the form which are plane waves with wavevector q and polarization vector . If a finite crystal is considered or if periodic boundary conditions are applied, the wavevector is restricted to a sequence of discrete and equidistant values which are, however, very close to each other. Thus, for practical work q can be treated as a continuous variable. The polarization vectors are, in general, different for every atom κ. Moreover, they depend on q and for each wavevector there are 3N different modes of vibration characterized not only by different 's but also by different vibrational frequencies ω. Hence, equation (2.1.2.10) can be written more specifically aswhere the running index labels the different fundamental vibrations or phonons. The upper index distinguishes two waves with identical frequencies which are travelling in opposite directions.

#### 2.1.2.2. Stability conditions

| top | pdf |

Not all of the elements of the force-constant matrix are independent. From its definition, equation (2.1.2.5), it is clear that the force-constant matrix is symmetric: Moreover, there are general stability conditions arising from the fact that a crystal as a whole is in mechanical equilibrium: If a macroscopic crystal is rigidly translated by a vector , no interatomic interactions are affected and, hence, the force acting on any particular atom must vanish: and, consequently, This relation is known as the condition of translational invariance.

In a similar way, it is argued that no interatomic interactions are affected when the crystal is rigidly rotated by infinitesimal amounts about arbitrary axes. This condition of rotational invariance leads to the following restrictions for the force constants: for all and .

In mechanical equilibrium, there must not be any strains within the crystal. The conditions of an unstrained crystal are also known as Huang conditions and may be formulated as for all .

All these stability conditions are independent of the particular crystal structure. There are other restrictions that are due to the symmetry of the atomic arrangement. They will be considered in detail in Section 2.1.3.

#### 2.1.2.3. The dynamical matrix

| top | pdf |

If the ansatz (2.1.2.10a) is inserted into the equation of motion (2.1.2.9), the following eigenvalue equation is obtained: The summation over all primitive cells on the right-hand side of equation (2.1.2.15) yields the Fourier-transformed force-constant matrix which is independent of l for infinite crystals. contains all interactions of type atoms with type atoms. Using this notation, equation (2.1.2.15) reduces toIf for a given vibration characterized by we combine the three-dimensional polarization vectors of all atoms within a primitive cell to a 3N-dimensional polarization vector ,and simultaneously the matrices to a matrix F(q)equation (2.1.2.17) can be written in matrix notation and takes the simple formwhere the diagonal matrix contains the masses of all atoms. The matrix is called the dynamical matrix. It contains all the information about the dynamical behaviour of the crystal and can be calculated on the basis of specific models for interatomic interactions. In analogy to the matrices , we introduce the submatrices of the dynamical matrix: Owing to the symmetry of the force-constant matrix, the dynamical matrix is Hermitian:1or more specifically Obviously, the squares of the vibrational frequency and the polarization vectors are eigenvalues and corresponding eigenvectors of the dynamical matrix. As a direct consequence of equation (2.1.2.20), the eigenvalues are real quantities and the following relations hold:Moreover, the eigenvectors are mutually orthogonal and can be chosen to be normalized.

#### 2.1.2.4. Eigenvalues and phonon dispersion, acoustic modes

| top | pdf |

The wavevector dependence of the vibrational frequencies is called phonon dispersion. For each wavevector q there are 3N fundamental frequencies yielding 3N phonon branches when is plotted versus q. In most cases, the phonon dispersion is displayed for wavevectors along high-symmetry directions. These dispersion curves are, however, only special projections of the dispersion hypersurface in the four-dimensional q–ω space. As a simple example, the phonon dispersion of b.c.c. hafnium is displayed in Fig. 2.1.2.3. The wavevectors are restricted to the first Brillouin zone (see Section 2.1.3.1) and the phonon dispersion for different directions of the wavevector are combined in one single diagram making use of the fact that different high-symmetry directions meet at the Brillouin-zone boundary. Note that in Fig. 2.1.2.3, the moduli of the wavevectors are scaled by the Brillouin-zone boundary values and represented by the reduced coordinates ξ. Owing to the simple b.c.c. structure of hafnium with one atom per primitive cell, there are only three phonon branches. Moreover, for all wavevectors along the directions [00ξ] and [ξξξ], two exhibit the same frequencies – they are said to be degenerate. Hence in the corresponding parts of Fig. 2.1.2.3 only two branches can be distinguished.

 Figure 2.1.2.3 | top | pdf |Phonon dispersion of b.c.c. hafnium for wavevectors along the main symmetry directions of the cubic structure. The symbols represent experimental data obtained by inelastic neutron scattering and the full lines are the results of the model. From Trampenau et al. (1991). Copyright (1991) by the American Physical Society.

Whereas in this simple example the different branches can be separated quite easily, this is no longer true for more complicated crystal structures. For illustration, the phonon dispersion of the high-Tc superconductor Nd2CuO4 is shown in Fig. 2.1.2.4 for the main symmetry directions of the tetragonal structure (space group , seven atoms per primitive cell). Note that in many publications on lattice dynamics the frequency is used rather than the angular frequency ω.

 Figure 2.1.2.4 | top | pdf |Phonon dispersion of Nd2CuO4 along the main symmetry directions of the tetragonal structure. The symbols represent experimental data obtained by inelastic neutron scattering and the full lines are drawn to guide the eye. Reprinted from Pintschovius et al. (1991), copyright (1991), with permission from Elsevier.

The 21 phonon branches of Nd2CuO4 with their more complicated dispersion reflect the details of the interatomic interactions between all atoms of the structure. The phonon frequencies ν cover a range from 0 to 18 THz. In crystals with strongly bonded molecular groups, like SiO4 tetrahedra in quartz or SO4 tetrahedra in sulfates, for example, the highest frequencies are found near 35 THz and correspond to bond-stretching vibrations. Soft materials like organic molecular crystals, on the other hand, exhibit a large number of phonon branches within a rather small frequency range which cannot easily be separated. Deuterated naphthalene (C10D8) is a well investigated example. The low-frequency part of its phonon dispersion is shown in Fig. 2.1.2.5.

 Figure 2.1.2.5 | top | pdf |Low-frequency part of the phonon dispersion of deuterated naphthalene at 6 K. The symbols represent experimental data obtained by inelastic neutron scattering and the full lines are drawn to guide the eye. Reproduced with permission from Natkaniec et al. (1980). Copyright (1980) IOP Publishing Limited.

Whereas neutron inelastic scattering is the most powerful method for the determination of phonons at arbitrary wavevectors, long wavelength phonons may also be detected by optical spectroscopy. The determination of phonon frequencies alone is, however, not sufficient for a concise determination of dispersion branches. Rather, individual phonons have to be assigned uniquely to one of the 3N branches, and this may prove a rather hard task if N is large. Here, symmetry considerations of eigenvectors are of special importance since phonons belonging to the same branch must exhibit the same symmetry properties. Moreover, inspection of Figs. 2.1.2.3 to 2.1.2.5 shows that some of the branches cross each other and others do not. It is a general statement that crossing is only allowed for branches with different symmetries – a property which yields a classification scheme for the different phonon branches. The symmetry of fundamental vibrations of a lattice will be discussed in some detail in Section 2.1.3.

In the limit of long wavelengths, there are always three particular modes with identical polarization vectors for all atoms, which will be considered in the following. At exactly (the Γ point) or infinite wavelength, the eigenvalue equation (2.1.2.15) reduces to One immediately recognizes that there are special solutions withi.e. the (mass-weighted) eigenvectors of all atoms are identical. There are three orthogonal eigenvectors of this kind and the displacement pattern of such phonons corresponds to rigid translations of the whole lattice along the three orthogonal coordinates in direct space. These motions do not affect any interatomic interaction. Hence, there is no change in potential energy and the condition of translational invariance (cf. Section 2.1.2.2) guarantees that the frequencies of these modes are zero:The phonon branches that lead to zero frequency at the Γ point () are called acoustic, whereas all other branches are called optic. The dispersion of acoustic branches in the vicinity of the Γ point can be investigated by expanding the phase factor in equation (2.1.2.15) in powers of q. Using (2.1.2.28) one obtainsNeglecting higher-order terms, summing up both sides of equation (2.1.2.30) over κ and multiplying by yieldsM being the total mass of all atoms within the primitive cell (). The first term on the right-hand side is zero according to equation (2.1.2.29). The second term vanishes due to the symmetry property of the force-constant matrices, equation (2.1.2.23). Hence (2.1.2.31) is simplified toThe right-hand side no longer depends on the moduli of the wavevector and displacement but only on their orientation with respect to the crystal lattice. Consequently, acoustic dispersion curves always leave the Γ point as a straight line with a constant slope ().

The displacement pattern of these long-wavelength modes corresponds to a continuous deformation of a rigid body. Hence, acoustic phonons near the Γ point can be regarded as sound waves and the slope of the dispersion curve is given by the corresponding sound velocity, Sound velocities, on the other hand, can be calculated from macroscopic elastic constants using the theory of macroscopic elasticity (cf. Chapter 1.3 ). Thus we are able to correlate macroscopic and microscopic dynamic properties of crystals. Using the generalized Hooke's law, the equation of motion for the dynamic deformation of a macroscopic body may be written as being the macroscopic density, the ith Cartesian component of the deformation and () the symmetric tensor of elastic stiffnesses, which is discussed in detail in Chapter 1.3 . The solution of this differential equation using plane waves,leads to the following relation: If we define the components of the propagation tensor byequation (2.1.2.36) may be written as the eigenvector equationFor a given propagation direction as defined by the Cartesian components of q, the eigenvectors of the corresponding propagation tensor yield the polarization of three mutually orthogonal deformation waves. Its eigenvalues are related to the respective sound velocities . If the tensor of elastic stiffnesses is known, the elements of and, hence, the velocity of elastic (sound) waves can be calculated for arbitrary propagation directions (see Section 1.3.4 ). These data, in turn, allow the prediction of the slopes of acoustic phonon dispersion curves near .

#### 2.1.2.5. Eigenvectors and normal coordinates

| top | pdf |

The plane-wave solutions (2.1.2.10) of the equations of motion form a complete set of orthogonal functions if q is restricted to the first Brillouin zone. Hence, the actual displacement of an atom (κl) can be represented by a linear combination of the : Since this displacement is an observable quantity, it must correspond to a real vector, not a complex one. Hence, the coefficients obey the relationand equation (2.1.2.39) reduces towhere If the displacement vectors are combined to form a 3N-dimensional vector in analogy to the formation of the eigenvector from the individual polarization vectors [equation (2.1.2.18)] we obtainThus, the atomic displacement is a linear combination of the eigenvectors of the dynamical matrix. The coefficients are called normal coordinates. They reflect the relative weight and amplitude of a particular vibrational mode () which is temperature-dependent and may be determined by statistical methods.

In terms of these normal coordinates, the Hamiltonian of the lattice (2.1.2.8) is reduced to a sum of independent harmonic oscillators. These are called phonons and may be regarded as quantum quasiparticles.( is the number of primitive cells within the crystal.)

#### 2.1.2.6. Amplitudes of lattice vibrations

| top | pdf |

Lattice vibrations that are characterized by both the frequencies and the normal coordinates are elementary excitations of the harmonic lattice. As long as anharmonic effects are neglected, there are no interactions between the individual phonons. The respective amplitudes depend on the excitation level and can be determined by quantum statistical methods. The energy levels of a lattice vibration () are those of a single harmonic oscillator: as illustrated in Fig. 2.1.2.6. The levels are equidistant and the respective occupation probabilities are given by Boltzmann statistics: In the quasiparticle description, this quantity is just the probability that at a temperature T there are n excited phonons of frequency . Moreover, in thermal equilibrium the average number of phonons is given by the Bose factor:and the corresponding contribution of these phonons to the lattice energy is The mean-square amplitude of the normal oscillator coordinate is obtained asAt high temperatures (), the phonon number, the corresponding energy and the amplitude approach the classical values of respectively. Note that occupation number, energy and amplitude merely depend on the frequency of the particular lattice vibration. The form of the corresponding eigenvector is irrelevant.

 Figure 2.1.2.6 | top | pdf |Energy levels of a quantum-mechanical harmonic oscillator.

#### 2.1.2.7. Density of states and the lattice heat capacity

| top | pdf |

The total energy stored in the harmonic phonon system is given by the sum over all phonon states (): Related thermodynamic quantities like the internal energy or the heat capacity are determined by the frequency distribution of the lattice vibrations rather than by details of the phonon dispersion. Hence, it is useful to introduce the phonon density of states in such a way that is the number of phonons with frequencies between and . Using this quantity, the sum in (2.1.2.52) may be replaced by an integral expression:Here, is the energy at . The derivative with respect to temperature provides the lattice heat capacity at constant volume: As an example, Fig. 2.1.2.7 displays the phonon dispersion of GaAs as determined by inelastic neutron scattering along with the phonon density of states. Obviously, even in this relatively simple substance (DOS) exhibits a rather complicated multi-peak structure. Integral properties like the heat capacity are, however, not very sensitive to details of . There are two well known approximations for that are able to reproduce some prominent features of :

 Figure 2.1.2.7 | top | pdf |Phonon dispersion and density of states for GaAs. The experimental data are from Strauch & Dorner (1990); the full lines and the density of states (DOS) are results of ab initio model calculations by Giannozzi et al. (1991). From Giannozzi et al. (1991). Copyright (1991) by the American Physical Society.

(1) The Einstein model.

In the Einstein model, it is assumed that all phonons exhibit the same frequency (Einstein oscillator) and is represented by a delta function: Consequently, the heat capacity turns out to bewhere we use the abbreviation which is the Einstein temperature.

At low temperatures, this model predicts an exponential temperature dependence of the heat capacity , which does not correspond to the experimental findings in most substances. Here, the Debye model provides a significant improvement.

(2) The Debye model.

In contrast to the Einstein model, which takes only one optic mode into account, the Debye model is restricted to acoustic modes that exhibit a linear dispersion close to the Γ point (see Section 2.1.2.4). Neglecting any deviation from linear behaviour, we get the simple result that the density of states is proportional to the square of the phonon frequency. The total number of phonon states is, however, given by , which is the number of all dynamical degrees of freedom of the whole system. Consequently, the frequency spectrum is assumed to be limited to frequencies below a particular value according to This limiting frequency is called the Debye frequency and is related to an appropriate average of (longitudinal and transverse) sound velocities and exhibits large values for hard materials. Fig. 2.1.2.8 compares schematically the true phonon density of states with the Debye approximation. The density of phonon states may thus be represented by and, correspondingly, the heat capacity is yielding a temperature dependence as shown in Fig. 2.1.2.9.

 Figure 2.1.2.8 | top | pdf |Schematic representation of the true phonon density of states (solid line) along with the Debye approximation (dotted line). Note that the areas under the two curves are identical.
 Figure 2.1.2.9 | top | pdf |Temperature dependence of the heat capacity at constant volume according to the Debye model.

is the Debye temperature, which is defined as At low temperatures, the heat capacity is proportional to , in excellent agreement with most experiments:

It is not surprising that the Debye model provides a reasonable description of the low-temperature heat capacity, since in this temperature regime well below the Debye temperature, optical phonons are hardly excited and the heat capacity is dominated by the low-frequency acoustic modes which are modelled exactly. At higher temperatures it is, however, necessary to take into account the thermal excitation of (in general less dispersive) optic modes. This can be achieved either by introducing a temperature dependence of the Debye temperature or by mixing a Debye term like (2.1.2.60) and Einstein terms like (2.1.2.56).

As an example, we consider the case of GaAs, the density of states of which is shown in Fig. 2.1.2.7. Obviously there are two very pronounced peaks at high frequencies, which are due to nearly dispersionless optical phonon branches. These modes may therefore be regarded as Einstein oscillators. The remaining acoustic branches lead to the more continuous part of the spectrum at lower frequencies, which may be approximated by a Debye law.

#### 2.1.2.8. Thermal expansion, compressibility and Grüneisen parameters

| top | pdf |

So far, we have always assumed that the crystal volume is constant. As long as we are dealing with harmonic solids, the thermal excitation of phonons does not result in a mean displacement of any atom. Consequently, thermal expansion cannot be understood in the harmonic approximation. It is due to the fact that there are anharmonic contributions to the lattice energy, i.e. third- and higher-order terms in the expansion with respect to atomic displacements [equation (2.1.2.4)]. Moreover, in an anharmonic lattice phonons are no longer independent elementary excitations. Rather, different lattice vibrations interact with each other leading to temperature-dependent frequency shifts, damping etc. Quantitatively, anharmonic effects may be analysed by means of perturbation theory, which is, however, beyond the scope of the present article. Details may be found, for example, in the monograph The Physics of Phonons (Reissland, 1973).

Some aspects of anharmonicity can, however, be discussed on the basis of the quasi-harmonic model. This approach makes use of the fact that the atomic interactions vary with the interatomic spacing and, hence, with the volume or, more generally, with any type of lattice deformation. The phonon frequencies will therefore depend on the deformation as well. Using the deformed lattice as a new reference frame for lattice dynamical calculations, the corresponding frequencies may be obtained again on the basis of a harmonic model with modified effective force constants. The comparison of phonons of both the original and the (arbitrarily) deformed lattice finally yields the partial derivatives of the frequencies with respect to the components of the strain tensor.

If the deformation is exclusively due to a change of temperature, the phonon frequency shifts can thus be related to the coefficients of thermal expansion. In this approximation, any intrinsic temperature dependence of phonon frequencies due to phonon interactions is neglected. Note that just those effects are, however, of particular importance if displacive phase transitions that are associated with soft phonon modes are considered.

In the quasi-harmonic approximation, we use the thermodynamic relation between the Helmholtz free energy A and the partition function Z: Z is given in terms of the energy levels of the independent harmonic oscillators: where Φ is the potential energy of the crystal if all atoms occupy their equilibrium positions. Hence, the following expression for A results: Elementary thermodynamics yields the pressure p as the partial derivative of A with respect to the volume at constant temperature: This relation may be generalized if not only volume changes are taken into account but also arbitrary deformations as described by the strain tensor : with the stress tensor . Note that the hydrostatic pressure p and the relative volume change are given by the traces of and , respectively: Taking the temperature derivative of (2.1.2.66), we obtain Using Euler's relations, the left-hand side may be written as where the tensor of the elastic stiffnesses and the tensor of thermal expansion have been introduced.

The free energy depends on the lattice deformations via the phonon frequencies. Hence, using (2.1.2.65), the right-hand side of (2.1.2.69) is evaluated as under the assumption that the phonon frequencies do not depend explicitly on the temperature. [ is the Bose factor, (2.1.2.47)].

Let us denote the contribution of a single phonon () to the heat capacity at constant volume by Then the combination of (2.1.2.69), (2.1.2.70) and (2.1.2.71) yields the result with the generalized-mode Grüneisen parametersThe set of equations (2.1.2.75) (for ) provide relations between the variation of phonon frequencies with the lattice deformations on the one hand and the tensors of elastic stiffnesses and thermal expansion on the other hand.

For cubic crystals, the tensor of the thermal expansion is diagonal,where β represents the coefficient of volume expansion. If Voigt's notation is used for the elastic stiffnesses (cf. Section 1.3.3.2.2 ), equation (2.1.2.75) reduces toIf there is an isotropic deformation, the shift of phonon frequencies may be described by an averaged-mode Grüneisen parameter: and (2.1.2.78) may be rewritten asIntroducing the mean Grüneisen parameter γ by the summation over all phonon states,we arrive at Remembering that in cubic crystals the expression represents the isothermal compressibility, we find the commonly used scalar form of equation (2.1.2.75): which relates the thermodynamic quantities thermal expansion, compressibility and heat capacity with the mean Grüneisen parameter. For most substances, γ exhibits values between 1 and 4 which are hardly temperature dependent. Hence, equation (2.1.2.84) may be regarded as an equation of state for solid systems.

Experimentally, it is almost impossible to determine the heat capacity at constant volume since the thermal expansion cannot be easily compensated. The more convenient quantity is , the heat capacity at constant pressure. There is a simple thermodynamic relation between the two quantities, and hence the following equation is obtained:

### References

Reissland, J. A. (1973). The physics of phonons. London: Wiley. (ISBN 0–471–71585–9.)