Phonons: Theory

To understand them we start by looking at the Taylor expansion of the total energy ($E$) around the set of equilibrium positions of the nuclei ($\{\mathbf{R}^0\}$)

$E(\{\mathbf{R}\})= E(\{\mathbf{R}^0\})+ \sum_{I\alpha} \frac{\partial E(\{\mathbf{R^0}\})}{\partial R_{I\alpha}} (R_{I\alpha}-R^0_{I\alpha})+ \sum_{I\alpha J\beta} \frac{\partial E(\{\mathbf{R}^0\})}{\partial R_{I\alpha} \partial R_{J\beta}} (R_{I\alpha}-R^0_{I\alpha}) (R_{J\beta}-R^0_{J\beta})+ \mathcal{O}(\mathbf{R}^3),$

where $\{\mathbf{R}\}$ the positions of the nuclei. The first term in the expansion corresponds to the forces

$F_{I\alpha} (\{\mathbf{R}^0\}) = - \left. \frac{\partial E(\{\mathbf{R}\})}{\partial R_{I\alpha}} \right|_{\mathbf{R} =\mathbf{R^0}}$,

and the second to the second-order force-constants

$C_{I\alpha J\beta} (\{\mathbf{R}^0\}) = \left. \frac{\partial E(\{\mathbf{R}\})}{\partial R_{I\alpha} \partial R_{J\beta}} \right|_{\mathbf{R} =\mathbf{R^0}} = - \left. \frac{\partial F_{I\alpha}(\{\mathbf{R}\})}{\partial R_{J\beta}} \right|_{\mathbf{R} =\mathbf{R^0}}.$

We can define a variable that corresponds to the displacement of the atoms with respect to the equilibrium position $R_{I\alpha} = R^0_{I\alpha}+u_{I\alpha}$ which leads to

$E(\{\mathbf{R}\})= E(\{\mathbf{R}^0\})+ \sum_{I\alpha} -F_{I\alpha} (\{\mathbf{R}^0\}) u_{I\alpha}+ \sum_{I\alpha J\beta} C_{I\alpha J\beta} (\{\mathbf{R}^0\}) u_{I\alpha} u_{J\beta} + \mathcal{O}(\mathbf{R}^3)$

If we take the system to be in equilibrium, the forces on the atoms are zero and then the Hamiltonian of the system is

$H = \frac{1}{2} \sum_{I\alpha} M_I \ddot{u}^2_{I\alpha} + \frac{1}{2} \sum_{I\alpha J\beta} C_{I\alpha J\beta} u_{I\alpha} u_{J\beta},$

and the equation of motion

$M_I \ddot{u}^2_{I\alpha} = - C_{I\alpha J\beta} u_{J\beta}$

Using the following ansatz

$\mathbf{u}^\mu_{I\alpha}(\mathbf{q},t) = \frac{1}{2} \frac{1}{\sqrt{M_I}} \left\{ A^\mu(\mathbf{q}) \xi^{\mu }_{I\alpha}(\mathbf{q}) e^{ i [\mathbf{q} \cdot \mathbf{\mathbf{R}}_I -\omega_\mu(\mathbf{q})t]}+ A^{\mu*}(\mathbf{q}) \xi^{\mu*}_{I\alpha}(\mathbf{q}) e^{-i [\mathbf{q} \cdot \mathbf{\mathbf{R}}_I -\omega_\mu(\mathbf{q})t]} \right\}$

where $\xi^{\mu }_{I\alpha}(\mathbf{q})$ are the phonon mode eigenvectors. Replacing we obtain the following eigenvalue problem

$\sum_{J\beta} D_{I\alpha J\beta} (\mathbf{q}) \xi^{\mu }_{J\beta}(\mathbf{q}) = \omega^\mu(\mathbf{q})^2 \xi^{\mu }_{I\alpha}(\mathbf{q})$

with

$D_{I\alpha J\beta} (\mathbf{q}) = \frac{1}{\sqrt{M_I M_J}} C_{I\alpha J\beta} e^{i\mathbf{q} \cdot (\mathbf{R}_J-\mathbf{R}_I)}$

the dynamical matrix in the harmonic approximation. Now by solving the eigenvalue problem above we can obtain the phonon modes $\xi^{\mu }_{I\alpha}(\mathbf{q})$ and frequencies $\omega^\mu(\mathbf{q})^2$ at any arbitrary q point.

Finite differences

The second-order force constants are computed using finite differences of the forces when each ion is displaced in each independent direction. This is done by creating systems with finite ionic displacement of atom $a$ in direction $i$ with magnitude $\lambda$, computing the orbitals $\psi^{u^a_i}_{\lambda}$ and the forces for these systems. The second-order force constants are then computed using

$C_{I\alpha J\beta}= \frac{\partial^2E}{\partial u_{I\alpha} \partial u_{J\beta}}= -\frac{\partial F_{I\alpha}}{\partial u_{J\beta}} \approx -\frac{ \left( \mathbf{F}[\{\psi^{u_{J\beta}}_{\lambda/2}\}]- \mathbf{F}[\{\psi^{u_{J\beta}}_{-\lambda/2}\}] \right)_{I\alpha}}{\lambda}, \quad {I=1,..,N_\text{atoms}} \quad {J=1,..,N_\text{atoms}} \quad {\alpha=x,y,z} \quad {\beta=x,y,z}$

where $u^a_i$ corresponds to the displacement of atom $a$ in the cartesian direction $i$ and $\mathbf{F}[\psi]$ retrieves the set of forces acting on all the ions given the $\psi_{n\mathbf{k}}$ orbitals.

Similarly, the internal strain tensor is

$\Xi_{I\alpha l}=\frac{\partial^2 E}{\partial u_{I\alpha} \partial \eta_l}= \frac{\partial \sigma_l}{\partial u_{I\alpha}} \approx \frac{ \left( \mathbf{\sigma}[\{\psi^{u_{I\alpha}}_{\lambda/2}\}]- \mathbf{\sigma}[\{\psi^{u_{I\alpha}}_{-\lambda/2}\}] \right)_l }{\lambda} ,\qquad {l=xx, yy, zz, xy, yz, zx} \quad {\alpha=x,y,z} \quad {J=1,..,N_\text{atoms}}$

where $\mathbf{\sigma}[\psi_{n\mathbf{k}}]$ computes the strain tensor given the $\psi_{n\mathbf{k}}$ orbitals.