Theory: Harmonic vibrations#

Info

This section functions as a brief introduction the basics of thermal transport with phonons. Complete phonon tutorials are available in Phonons with FHI-vibes.

A Taylor expansion of the potential energy surface \(E(\{\boldsymbol{R}\})\) around equilibrium \(\boldsymbol{R}^0\) is usually made for describing the lattice vibrations in a solid:

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

The linear term vanishes, since no forces \(\boldsymbol{F} = - \nabla E\) are acting on the system in equilibrium \(\boldsymbol{R}^0\). For the study of lattice dynamics, the harmonic approximation is usually made, which means preserving the Taylor series up to second-order term, the Hessian \(\Phi_{IJ} = \frac{\partial^2 E}{\partial \boldsymbol{R}_I\partial \boldsymbol{R}_J}\). Higher-order residual terms \(\mathcal{O}(\Delta\boldsymbol{R}^3)\) are neglected in the harmonic approximation, but play an important role in thermal transport, which will be introduced in later tutorials.

Accessing the Hessian \(\Phi_{IJ} = \frac{\partial^2 E}{\partial \boldsymbol{R}_I\partial \boldsymbol{R}_J}\) involves some additional computations of the derivative of forces \(\boldsymbol{F}\), with respect to the nuclear displacements. One can either use Density Functional Perturbation Theory (DFPT) ¹ to compute the response or one can circumvent this problem by the finite displacements method

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

as we will do in this tutorial. As shown in Eq. (2), the Hessian describes how the force acting on an atom \(\boldsymbol{R}_J\) changes if we displace an atom \(\boldsymbol{R}_I\). In the case of periodic boundary conditions, as well as the atoms in the unit cell \(\boldsymbol{R}_J\), we also need to account for the periodic images \(\boldsymbol{R}_{J} + \boldsymbol{R}_{L}\), where \(\boldsymbol{R}_L\) is the Bravais lattice vector. 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 \(\boldsymbol{R}_{IJ}\), therefore 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 in the mass-scaled Fourier transform of the Hessian:

\[D_{IJ}^{\alpha\beta}(\boldsymbol{q}) = \sum_{L} \frac{1}{\sqrt{M_I M_J}} \Phi_{IJ}^{\alpha \beta} \exp(i\boldsymbol{q}(\boldsymbol{R}_I - \boldsymbol{R}_J - \boldsymbol{R}_L))\]

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

\[ D(\boldsymbol{q}) \, \boldsymbol{e}_s (\boldsymbol{q}) = \omega_s^2(\boldsymbol{q}) \, \boldsymbol{e}_s (\boldsymbol{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 \(\boldsymbol{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(\boldsymbol{q})\) and eigenvectors \(\boldsymbol{e}_s(\boldsymbol{q})\) of the dynamical matrix \(D(\boldsymbol{q})\) completely describe the dynamics of the system (in the harmonic approximation), which is nothing more than a superposition of harmonic oscillators, one for each mode, i.e., for each eigenvalue \(\omega_s (\boldsymbol{q})\).

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