Phonons#

Phonons: Harmonic vibrations in solids#

To determine the vibrations in a solid, we approximate the potential energy surface for the nuclei by performing a Taylor expansion of the total energy $\mathcal V({\bf R})$ around the equilibrium positions:

The linear term vanishes, since no forces $\vec{F} = - \nabla \mathcal V$ are acting on the system in equilibrium $\vec{R}^0$. Assessing the Hessian $\Phi_{IJ} = \frac{\partial^2 \mathcal V}{\partial \vec{R}_I\partial \vec{R}_J}$ involves some additional complications: In contrast to the forces $\vec{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 its response to nuclear displacements. One can either use Density Functional Perturbation Theory (DFPT) [Baroni2001] to compute the response or one can circumvent this problem by performing the second order derivative numerically by finite differences

as we will do in this tutorial. The definition in Eq.$~\eqref{eq:FinDiff}$ is helpful to realize that the Hessian describes a coupling between different atoms, i.e., how the force acting on an atom $\vec{R}_J$ changes if we displace atom $\vec{R}_I$, as you have already learned in tutorial 1. However, an additional complexity arises in the case of periodic boundary conditions, since beside the atoms in the unit cell $\vec{R}_J$ we also need to account for the periodic images $\vec{R}_{J'}$. Accordingly, the Hessian is in principle a matrix of infinite size. In non-ionic crystals, however, the interaction between two atoms$~I$ and $J$ quickly decays with their distance$~\vec{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:

In reciprocal space [AshcroftMermin], this dynamical matrix determines the equation of motion for such a periodic array of harmonic atoms for each reciprocal vector$~\vec{q}$:

The dynamical matrix has dimension $3N_\t{A} \times 3N_\t{A}$, where $N_\t{A}$ is the number of atoms in the primitive unit cell. Equation$~\eqref{eq:eigenproblem}$ thus constitutes an eigenvalue problem with $3N_\t{A}$ solutions at each $\vec 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_\t{A} - 3)$ branches are denoted as optical branches.

The eigenvalues$~\omega_s^2(\vec{q})$ and eigenvectors$~\vec e_s(\vec{q})$ of the dynamical matrix$~D(\vec{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 (\vec{q})$.

From the set of eigenvalues $\{ \omega_s (\vec{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 + \t d \omega]$:

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$~\vec{q}$. This is much easier to handle both in numerical and in analytical models. For instance, we can compute the associated thermodynamic potential¹, i.e., the (harmonic) Helmholtz free energy²

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

For a comprehensive introduction to the field of lattice dynamics and its foundations, we refer to [BornHuang].

To compute the quantities introduced above, we will use FHI-vibes, which uses the package phonopy [Togo2015] as a backend to compute vibrational properties via the finite-displacements method as outlined above. Please note that phonopy makes extensive use of symmetry analysis [Parlinski1997], which allows to reduce numerical noise and to speed up the calculations considerably. Our system of choice will be fcc-diamond Silicon (you can run the tutorial as well with LJ-Argon).

Recap on solid state physics#

Brillouin Zone#

As you know, the periodicity of the atomic positions in a lattice is reflected by the fact that it is sufficient to look at $\bf q$ values within the first Brillouin zone of the reciprocal lattice: A wave vector $\bf q$ that lies outside the first Brillouin zone corresponds to a wave whose wavelength is shorter than the distance between periodic images of the same atom. It can thus be represented equally well by a wave with longer wavelength, i.e. a smaller wave vector ${\bf q}'$ taken from within the first Brillouin zone.

High Symmetry Points#

The Brillouin zone of our Silicon fcc diamond structure is displayed below

The labelled points correspond to $\bf q$ values of high symmetry. This means that there are symmetry operations in the point group of the lattice that leave this point invariant (up to a reciprocal space vector) [p. 218 in Dresselhaus].

You can list the high symmetry points of the lattice of your geometry with vibes by running

vibes info geometry geometry.in -v

...
Special k points:
G: [0. 0. 0.]
K: [0.375 0.375 0.75 ]
L: [0.5 0.5 0.5]
U: [0.625 0.25  0.625]
W: [0.5  0.25 0.75]
X: [0.5 0.  0.5]

Please note that the list of points is given in fractional coordinates as coefficients of the reciprocal lattice. For the meaning of the Symbols $\Gamma$ (G), $X$, etc., you can take a look at the Wikipedia article.

Bandstructure#

If we connect two or more $\bf q$ points from the Brillouin zone, solve the eigenvalue problem for any $\bf q$ point in between, and plot the obtained dispersions $\omega (q)$ versus $q$, we obtain the so-called phonon bandstructure. The bandstructure is typically computed for a path in the Brillouin zone that connects several or all of the high symmetry points.