Skip to content

Theory: Harmonic vibrations in solids

To determine the vibrations in a solid, we approximate the potential energy surface of the nuclei by performing a Taylor expansion of the total energy \(E\) around the equilibrium positions \(\mathbf{R}^0\):

\[\begin{align} E \left(\{\mathbf{R}^0 + \Delta \mathbf{R}\}\right) \approx E\left(\{\mathbf{R}^0\}\right) + \cancel{ \sum\limits_{I} \left.\frac{\partial E}{\partial \mathbf{R}_I}\right\vert_{\mathbf{R}^0} \Delta\mathbf{R}_{I} } + \frac{1}{2} \sum\limits_{I,J} \left.\frac{\partial^2 E}{\partial \mathbf{R}_I\partial \mathbf{R}_J}\right\vert_{\mathbf{R}^0} \Delta\mathbf{R}_{I}\Delta\mathbf{R}_{J} + \mathcal{O}(\Delta\mathbf{R}^3) \qquad (1) \end{align}\]

The linear term vanishes, since no forces \(\mathbf{F} = - \nabla E\) are acting on the system in equilibrium \(\mathbf{R}^0\). Assessing the Hessian \(\Phi_{IJ} = \frac{\partial^2 E}{\partial \mathbf{R}_I\partial \mathbf{R}_J}\) involves some additional complications: In contrast to the forces \(\mathbf{F}\), which only depend on the density, the Hessian \(\Phi_{IJ}\) also depends on its derivative with respect to the nuclear coordinates, i.e., on density response to nuclear displacements. One can either use Density Functional Perturbation Theory (DFPT) 2 to compute the response or one can circumvent this problem by performing the second order derivative numerically by finite differences

\[ \Phi_{IJ} = \left.\frac{\partial^2 E}{\partial \mathbf{R}_I\partial \mathbf{R}_J}\right\vert_{\mathbf{R}^0} = - \left.\frac{\partial }{\partial \mathbf{R}_I} \mathbf{F}_J\right\vert_{\mathbf{R}^0} \approx - \frac{ \mathbf{F}_J(\mathbf{R}_I^0 + \varepsilon \,\mathbf{d}_I)}{\varepsilon} \quad ,\qquad (2) \]

as we will do in this tutorial. The definition in Eq. (2) is helpful to realize that the Hessian describes a coupling between different atoms, i.e., how the force acting on an atom \(\mathbf{R}_J\) changes if we displace atom \(\mathbf{R}_I\), as you have already learned in tutorial 1. An additional complexity arises in the case of periodic boundary conditions, since beside the atoms in the unit cell \(\mathbf{R}_J\) we also need to account for the periodic images \(\mathbf{R}_{J'}\). Accordingly, the Hessian is in principle a matrix of infinite size. However, in non-ionic crystals the interaction between two atoms \(I\) and \(J\) quickly decays with their distance \(\mathbf{R}_{IJ}\), so that we can compute the Hessian from finite supercells, the size convergence of which must be accurately inspected.

Once the real-space representation of the Hessian is computed, we can determine the dynamical matrix by adding up the contributions from all periodic images \(J'\) in the mass-scaled Fourier transform of the Hessian:

\[ D_{IJ}({\mathbf{q}}) = \sum\limits_{J'} \frac{\text{e}^{\text{i}\left({\mathbf{q}}\cdot{\mathbf{R}_{JJ'}}\right)}}{\sqrt{M_I M_J}} \;\Phi_{IJ'} \quad .\qquad (3) \]

In reciprocal space 3, this dynamical matrix determines the equation of motion for such a periodic array of harmonic atoms for each reciprocal vector \(\mathbf{q}\):

\[ D(\mathbf{q}) \, \boldsymbol{e}_s (\mathbf{q}) = \omega_s^2(\mathbf{q}) \, \boldsymbol{e}_s (\mathbf{q}) \; . \qquad (4) \]

The dynamical matrix has dimension \(3N_\text{A} \times 3N_\text{A}\), where \(N_\text{A}\) is the number of atoms in the primitive unit cell. Equation (4) thus constitutes an eigenvalue problem with \(3N_\text{A}\) solutions at each \(\mathbf{q}\) point. The solutions are labelled by \(s\) and are denoted as phonon branches. The lowest three branches are commonly called the acoustic branches, whereas in solids with more than one atom in the primitive unit cell, the remaining \((3N_\text{A} - 3)\) branches are denoted as optical branches.

The eigenvalues \(\omega_s^2(\mathbf{q})\) and eigenvectors \(\boldsymbol{e}_s(\mathbf{q})\) of the dynamical matrix \(D(\mathbf{q})\) completely describe the dynamics of the system (in the harmonic approximation), which is nothing else than a superposition of harmonic oscillators, one for each mode, i.e., for each eigenvalue \(\omega_s (\mathbf{q})\).

From the set of eigenvalues \(\{ \omega_s (\mathbf{q}) \}\), also denoted as spectrum, the density of states (DOS) can be obtained by calculating the number of states in an infinitesimal energy window \([\omega, \omega + \,\text{d}\omega]\):

\[ g(\omega) = \sum_s\int\frac{\,\text{d}\mathbf{q}}{(2\pi)^3}\delta(\omega - \omega(\mathbf{q})) = \sum_s\int\limits_{\omega(\mathbf{q}) = \omega}\frac{\text{d}S}{(2\pi)^3}\frac{1}{\vert\nabla\omega(\mathbf{q})\vert}~. \qquad (5) \]

The DOS is a very useful quantity, since it allows to determine any integrals (the integrand of which only depends on \(\omega\)) by a simple integration over a one-dimensional variable \(\omega\) rather than a three-dimensional variable \(\mathbf{q}\). This is much easier to handle both in numerical and in analytical models. For instance, we can compute the associated thermodynamic potential1, i.e., the (harmonic) Helmholtz free energy:

\[ F^{\mathrm{ha}}(T,V) = \int \,\text{d}\omega\; g(\omega) \left ( \frac{\hbar\omega}{2} + k_\text{B}\, T~ \ln\left(1-\text{e}^{-\frac{\hbar\omega}{k_\text{B}\,T}}\right) \right)\;. \qquad (6) \]

In turn, this allows to calculate the heat capacity at constant volume:

\[ C_V = - T \left(\frac{\partial^2 F^{\mathrm{ha}}(T,V)}{\partial T^2} \right)_V \;. \qquad (7) \]

For a comprehensive introduction to the field of lattice dynamics and its foundations, we refer to the text book literature763.

To compute the quantities introduced above, we will utilize the program package FHI-vibes 8 that uses the package phonopy 5 as a backend to compute vibrational properties. Please note that phonopy makes extensive use of symmetry analysis 4, which allows to reduce numerical noise and to speed up the calculations considerably.

  1. Given that the Bose-Einstein distribution is used for the derivation of the harmonic free energy in this case, we get the correct quantum-mechanical result including zero-point effects by this means. ↩

  2. S. Baroni, et.al., Rev. Mod. Phys. 73, 515 (2001). ↩

  3. N. W. Ashcroft, N. D. Mermin, Solid State Physics, Saunders College Publishing, New York, (1976). ↩↩

  4. K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 78, 4063 (1997). ↩

  5. A. Togo, F. Oba, and I. Tanaka, Phys. Rev. B 78, 134106 (2008). ↩

  6. M. Dove, Introduction to Lattice Dynamics, Cambridge University Press (1993). ↩

  7. M. Born, and K. Huang, Dynamical Theory of Crystal Lattices, Oxford University Press (1962). ↩

  8. Knoop, F., Purcell, T., Scheffler, M., & Carbogno, C. (2020). FHI-vibes: Ab Initio Vibrational Simulations. The Journal of Open Source Software, 5(56). ↩