# 10. Formalism for program anphon¶

## 10.1. Dynamical matrix¶

The dynamical matrix is given by

(1)$D_{\mu\nu}(\kappa\kappa^{\prime};\boldsymbol{q}) = \frac{1}{\sqrt{M_{\kappa}M_{\kappa^{\prime}}}} \sum_{\ell^{\prime}}\Phi_{\mu\nu}(\ell\kappa;\ell^{\prime}\kappa^{\prime})\exp{\left[i\boldsymbol{q}\cdot(\boldsymbol{r}(\ell^{\prime})-\boldsymbol{r}(\ell))\right]},$

where $$M_{\kappa}$$ is the atomic mass of atom $$\kappa$$. By diagonalizing the dynamical matrix, one can obtain $$m$$ ($$=3N_{\kappa}$$) eigenvalues $$\omega_{\boldsymbol{q}j}^{2}$$ ($$j = 1, 2, \dots, m$$) and corresponding eigenvectors $$\boldsymbol{e}_{\boldsymbol{q}j}$$ for each $$\boldsymbol{q}$$ point. Here, $$\boldsymbol{e}_{\boldsymbol{q}j}$$ is a column vector consisting of atomic polarization $$e_{\mu}(\kappa;\boldsymbol{q}j)$$. Let $$D(\boldsymbol{q})$$ denote a matrix form of equation (1), the eigenvalues may be written as

(2)$\omega_{\boldsymbol{q}j}^{2} = (\boldsymbol{e}_{\boldsymbol{q}j}^{*})^{\mathrm{T}} D(\boldsymbol{q})\boldsymbol{e}_{\boldsymbol{q}j}.$

Next, we introduce $$m\times m$$ matrices $$\Lambda$$ and $$W$$ which are defined as $$\Lambda(\boldsymbol{q}) = \mathrm{diag} (\omega_{\boldsymbol{q}1}^{2},\dots,\omega_{\boldsymbol{q}m}^{2})$$ and $$W(\boldsymbol{q}) = (\boldsymbol{e}_{\boldsymbol{q}1},\dots,\boldsymbol{e}_{\boldsymbol{q}m})$$, respectively. Then, equation (2) can be denoted as

$\Lambda(\boldsymbol{q}) = W^{\dagger}(\boldsymbol{q})D(\boldsymbol{q})W(\boldsymbol{q}).$

When one needs to capture the LO-TO splitting near the zone-center by the supercell approach, it is necessary to add the non-analytic part of the dynamical matrix defined by

$D_{\mu\nu}^{\mathrm{NA}}(\kappa\kappa^{\prime};\boldsymbol{q}) = \frac{1}{\sqrt{M_{\kappa}M_{\kappa^{\prime}}}} \frac{4\pi e^{2}}{\Omega} \frac{(Z_{\kappa}^{*}\boldsymbol{q})_{\mu}(Z_{\kappa^{\prime}}^{*}\boldsymbol{q})_{\nu}}{\boldsymbol{q}\cdot\epsilon^{\infty}\boldsymbol{q}},$

where $$\Omega$$ is the volume of the primitive cell, $$Z_{\kappa}^{*}$$ is the Born effective charge tensor of atom $$\kappa$$, and $$\epsilon^{\infty}$$ is the dielectric constant tensor, respectively. In program anphon, either the Parlinski’s way [1] or the mixed-space approach [2] can be used. In the Parlinski’s approach (NONANALYTIC = 1), the total dynamical matrix is given by

$D(\boldsymbol{q}) + D^{\textrm{NA}}(\boldsymbol{q})\exp{(-q^{2}/\sigma^{2})},$

where $$\sigma$$ is a damping factor. $$\sigma$$ must be chosen carefully so that the non-analytic contribution becomes negligible at Brillouin zone boundaries. In the mixed-space approach (NONANALYTIC = 2), the total dynamical matrix is given by

$D(\boldsymbol{q}) + D^{\textrm{NA}}(\boldsymbol{q})\frac{1}{N}\sum_{\ell^{\prime}}\exp{\left[i\boldsymbol{q}\cdot(\boldsymbol{r}(\ell^{\prime})-\boldsymbol{r}(\ell))\right]}.$

The second term vanishes at commensurate $$\boldsymbol{q}$$ points other than $$\Gamma$$ point ($$\boldsymbol{q} = 0$$).

To include the non-analytic term, one needs to set NONANALYTIC > 0 and give appropriate BORNINFO and NA_SIGMA tags.

## 10.2. Group velocity¶

The group velocity of phonon mode $$\boldsymbol{q}j$$ is given by

$\boldsymbol{v}_{\boldsymbol{q}j} = \frac{\partial \omega_{\boldsymbol{q}j}}{\partial \boldsymbol{q}}.$

To evaluate the group velocity numerically, we employ a central difference where $$\boldsymbol{v}$$ may approximately be given by

$\boldsymbol{v}_{\boldsymbol{q}j} \approx \frac{\omega_{\boldsymbol{q}+\Delta\boldsymbol{q}j} - \omega_{\boldsymbol{q}-\Delta\boldsymbol{q}j}}{2\Delta\boldsymbol{q}}.$

If one needs to save the group velocities, please turn on the PRINTVEL-tag.

## 10.3. Thermodynamics functions¶

The specific heat at constant volume $$C_{\mathrm{v}}$$, the internal energy $$U$$, the vibrational entropy $$S$$, and the Helmholtz free energy $$F$$ of individual harmonic oscillator are given as follows:

\begin{align*} U &= \frac{1}{N_{q}}\sum_{\boldsymbol{q},j} \hbar\omega_{\boldsymbol{q}j} \left[\frac{1}{e^{\hbar\omega_{\boldsymbol{q}j}/kT} - 1} + \frac{1}{2}\right], \\ C_{\mathrm{v}} &= \frac{k}{N_{q}}\sum_{\boldsymbol{q},j} \left(\frac{\hbar\omega_{\boldsymbol{q}j}}{2kT}\right)^{2} \mathrm{cosech}^{2}\left(\frac{\hbar\omega_{\boldsymbol{q}j}}{2kT}\right),\\ S &= \frac{k}{N_{q}}\sum_{\boldsymbol{q},j} \left[\frac{\hbar\omega_{\boldsymbol{q}j}}{kT} \frac{1}{e^{\hbar\omega_{\boldsymbol{q}j}/kT} - 1} - \log{\left( 1 - e^{-\hbar\omega_{\boldsymbol{q}j}/kT}\right)}\right], \\ F &= \frac{1}{N_{q}}\sum_{\boldsymbol{q},j}\left[ \frac{\hbar\omega_{\boldsymbol{q}j}}{2} + kT\log{\left( 1 - e^{-\hbar\omega_{\boldsymbol{q}j}/kT}\right)} \right]. \end{align*}

Here, $$k$$ is the Boltzmann constant. These quantities will be saved in the PREFIX.thermo file.

## 10.4. Mean square displacement¶

The mean square displacement tensor of atom $$\kappa$$ is given by

\begin{align} \left< u_{\mu}(\kappa)u_{\nu}(\kappa) \right> & = \frac{\hbar}{2M_{\kappa}N_{q}}\sum_{\boldsymbol{q},j} \frac{1}{2\omega_{\boldsymbol{q}j}}\left(e_{\mu}(\kappa;\boldsymbol{q}j)e_{\nu}^{*}(\kappa;\boldsymbol{q}j)+ e_{\mu}^{*}(\kappa;\boldsymbol{q}j)e_{\nu}(\kappa;\boldsymbol{q}j)\right) \notag \\ & \hspace{25mm} \times \coth{\left(\frac{\hbar\omega_{\boldsymbol{q}j}}{2kT}\right)}. \end{align}

When PRINTMSD is turned on, the code print the diagonal part of the mean square displacement tensor

$\left< u_{\mu}^{2}(\kappa)\right> = \frac{\hbar}{M_{\kappa}N_{q}}\sum_{\boldsymbol{q},j}\frac{1}{\omega_{\boldsymbol{q}j}} |e_{\mu}(\kappa;\boldsymbol{q}j)|^{2} \left(n_{\boldsymbol{q}j}+\frac{1}{2}\right),$

where $$n_{\boldsymbol{q}j} = 1/(e^{\hbar\omega_{\boldsymbol{q}j}/kT}-1)$$ is the Bose-Einstein distribution function.

## 10.5. Phonon DOS¶

When KPMODE = 2, the program anphon saves the (one) phonon density of states (DOS) to the file PREFIX.dos. The one-phonon DOS is given by

$\mathrm{DOS}(\omega) = \frac{1}{N_{q}}\sum_{\boldsymbol{q},j}\delta(\omega - \omega_{\boldsymbol{q}j}).$

If PDOS = 1 is given, the program also prints the atom-projected phonon DOS which is given by

$\mathrm{PDOS}(\kappa;\omega) = \frac{1}{N_{q}}\sum_{\boldsymbol{q},j}|\boldsymbol{e}(\kappa;\boldsymbol{q}j)|^{2}\delta(\omega - \omega_{\boldsymbol{q}j}).$

In addition, TDOS-tag is available to compute the two-phonon DOS defined by

$\mathrm{DOS2}(\omega;\boldsymbol{q};\pm) = \frac{1}{N_{q}}\sum_{\boldsymbol{q}_{1},\boldsymbol{q}_{2}, j_{1}, j_{2}} \delta(\omega\pm\omega_{\boldsymbol{q}_{1}j_{1}}-\omega_{\boldsymbol{q}_{2}j_{2}})\delta_{\boldsymbol{q}\pm\boldsymbol{q}_{1},\boldsymbol{q}_{2}+\boldsymbol{G}},$

where $$\boldsymbol{G}$$ is a reciprocal lattice vector. The sign $$\pm$$ correspond to absorption and emission processes, respectively. Please note that the computation of the two-phonon DOS can be expensive especially when $$N_{q}$$ or $$N_{\kappa}$$ is large.

## 10.6. (Atomic) participation ratio¶

Participation ratio (PR) and atomic participation ratio (APR) defined in the following may be useful to analyze the localized nature of the phonon mode $$\boldsymbol{q}j$$.

• Participation ratio (PR)
$PR_{\boldsymbol{q}j} = \left(\sum_{\kappa}^{N_{\kappa}} \frac{|\boldsymbol{e}(\kappa;\boldsymbol{q}j)|^{2}}{M_{\kappa}}\right)^{2} \Bigg/ N_{\kappa} \sum_{\kappa}^{N_{\kappa}} \frac{|\boldsymbol{e}(\kappa;\boldsymbol{q}j)|^{4}}{M_{\kappa}^{2}}$
• Atomic participation ratio (APR)
$APR_{\boldsymbol{q}j,\kappa} = \frac{|\boldsymbol{e}(\kappa;\boldsymbol{q}j)|^{2}}{M_{\kappa}} \Bigg/ \left( N_{\kappa} \sum_{\kappa}^{N_{\kappa}} \frac{|\boldsymbol{e}(\kappa;\boldsymbol{q}j)|^{4}}{M_{\kappa}^{2}} \right)^{1/2}$

For an extended eigenmode, the PR value is of order 1, whereas for a localized eigenmodes PR is of order $$1/N_{\kappa}$$ [3]. APR is an atomic decomposition of PR that satisfies $$PR_{\boldsymbol{q}j} = \sum_{\kappa} (APR_{\boldsymbol{q}j,\kappa})^{2}$$. To print the PR and APR, please set MODE = phonons and PRINTPR = 1 in the &analysis entry field.

## 10.7. Scattering phase space¶

When KPMODE = 2 and SPS = 1, the three-phonon scattering phase space $$P_{3}$$ is calculated and saved to the file PREFIX.sps. $$P_{3}$$ is defined as

$P_{3}(\boldsymbol{q}j) = \frac{1}{3m^{3}} (2P_{3}^{(+)}(\boldsymbol{q}j) + P_{3}^{(-)}(\boldsymbol{q}j)),$

where $$m$$ is the number of phonon branches and

$P_{3}^{(\pm)}(\boldsymbol{q}j) = \frac{1}{N_{q}}\sum_{\boldsymbol{q}_{1},\boldsymbol{q}_{2}, j_{1}, j_{2}}\delta(\omega_{\boldsymbol{q}j}\pm\omega_{\boldsymbol{q}_{1}j_{1}}-\omega_{\boldsymbol{q}_{2}j_{2}})\delta_{\boldsymbol{q}\pm\boldsymbol{q}_{1},\boldsymbol{q}_{2}+\boldsymbol{G}}.$

anphon also print the total scattering phase space

$P_{3} = \frac{1}{N_{q}}\sum_{\boldsymbol{q}j} P_{3}(\boldsymbol{q}j).$

When SPS = 2, the three-phonon scattering phase space with the occupation factor $$W_{3}^{(\pm)}$$ will be calculated and saved to the file PREFIX.sps_Bose. $$W_{3}^{(\pm)}$$ is defined as

$\begin{split}W_{3}^{(\pm)}(\boldsymbol{q}j) = \frac{1}{N_{q}}{\sum_{\boldsymbol{q}_{1},\boldsymbol{q}_{2}, j_{1}, j_{2}}} \left\{ \begin{array}{c} n_{2} - n_{1} \\ n_{1} + n_{2} + 1 \end{array} \right\} \delta(\omega_{\boldsymbol{q}j}\pm\omega_{\boldsymbol{q}_{1}j_{1}}-\omega_{\boldsymbol{q}_{2}j_{2}})\delta_{\boldsymbol{q}\pm\boldsymbol{q}_{1},\boldsymbol{q}_{2}+\boldsymbol{G}}.\end{split}$

Here, $$n_{1}=n(\omega_{\boldsymbol{q}_{1}j_{1}})$$ and $$n_{2}=n(\omega_{\boldsymbol{q}_{2}j_{2}})$$ where $$n(\omega) = \frac{1}{e^{\hbar\omega/k_B T}-1}$$ is the Bose-Einstein distribution function. Since $$n(\omega)$$ is temperature dependent, $$W_{3}^{(\pm)}$$ is also temperature dependent. The file PREFIX.sps_Bose contains $$W_{3}^{(\pm)}$$ for all phonon modes at various temperatures specified with TMIN, TMAX, and DT tags.

## 10.8. Grüneisen parameter¶

The mode Grüneisen parameter, defined as $$\gamma_{\boldsymbol{q}j} = - \frac{\partial \log{\omega_{\boldsymbol{q}j}}}{\partial \log{V}}$$, is calculated by

$\gamma_{\boldsymbol{q}j}= -\frac{(\boldsymbol{e}_{\boldsymbol{q}j}^{*})^{\mathrm{T}} \delta D(\boldsymbol{q})\boldsymbol{e}_{\boldsymbol{q}j}}{6\omega_{\boldsymbol{q}j}},$

where $$\delta D(\boldsymbol{q})$$ is a change in the dynamical matrix due to a volume change $$\delta V$$, which is given by

\begin{align} \delta D_{\mu\nu}(\kappa\kappa^{\prime};\boldsymbol{q}) &= \frac{1}{\sqrt{M_{\kappa}M_{\kappa^{\prime}}}} \sum_{\ell^{\prime}}\delta\Phi_{\mu\nu}(\ell\kappa;\ell^{\prime}\kappa^{\prime})\exp{\left[i\boldsymbol{q}\cdot(\boldsymbol{r}(\ell^{\prime})-\boldsymbol{r}(\ell))\right]},\\ \delta\Phi_{\mu\nu}(\ell\kappa;\ell^{\prime}\kappa^{\prime}) &= \sum_{\ell^{\prime\prime},\kappa^{\prime\prime},\lambda}\Phi_{\mu\nu\lambda}(\ell\kappa;\ell^{\prime}\kappa^{\prime};\ell^{\prime\prime}\kappa^{\prime\prime})r_{\lambda}(\ell^{\prime\prime}\kappa^{\prime\prime}). \end{align}

Please set GRUNEISEN = 1 and give an appropriate FCSXML file containing cubic IFCs to print Grüneisen parameters.

## 10.9. Anharmonic self-energy¶

The anharmonic self-energy due to cubic anharmonicity to the lowest order is given by

(3)$\begin{split}\Sigma_{\boldsymbol{q}j}(i\omega_m) &= \frac{1}{2\hbar^{2}}\sum_{\boldsymbol{q}_{1},\boldsymbol{q}_{2}}\sum_{j_{1},j_{2}} |V^{(3)}_{-\boldsymbol{q}j,\boldsymbol{q}_{1}j_{1},\boldsymbol{q}_{2}j_{2}}|^{2} \notag \\ & \times \left[ \frac{n_{1}+n_{2} + 1}{i\omega_{m} + \omega_{1} + \omega_{2}} - \frac{n_{1}+n_{2} + 1}{i\omega_{m} - \omega_{1} - \omega_{2}} + \frac{n_{1}-n_{2}}{i\omega_{m} - \omega_{1} + \omega_{2}} - \frac{n_{1}-n_{2}}{i\omega_{m} + \omega_{1} - \omega_{2}} \right],\end{split}$

where $$i\omega_{m}$$ is the Matsubara frequency. In equation (3), we simply denoted $$\omega_{\boldsymbol{q}_{i}j_{i}}$$ as $$\omega_{i}$$. The matrix element $$V^{(3)}$$ is given by

$\begin{split}V^{(3)}_{\boldsymbol{q}j,\boldsymbol{q}^{\prime}j^{\prime},\boldsymbol{q}^{\prime\prime}j^{\prime\prime}} & = \left( \frac{\hbar}{2N_{q}}\right)^{\frac{3}{2}} \frac{1}{\sqrt{\omega_{\boldsymbol{q}n}\omega_{\boldsymbol{q}^{\prime}j^{\prime}}\omega_{\boldsymbol{q}^{\prime\prime}j^{\prime\prime}}}} \sum_{\ell,\ell^{\prime},\ell^{\prime\prime}} \exp{\left[\mathrm{i}(\boldsymbol{q}\cdot\boldsymbol{r}(\ell)+\boldsymbol{q}^{\prime}\cdot\boldsymbol{r}(\ell^{\prime})+\boldsymbol{q}^{\prime\prime}\cdot\boldsymbol{r}(\ell^{\prime\prime}))\right]} \notag \\ & \times \sum_{\kappa,\kappa^{\prime},\kappa^{\prime\prime}} \frac{1}{\sqrt{M_{\kappa}M_{\kappa^{\prime}}M_{\kappa^{\prime\prime}}}} \sum_{\mu,\nu,\lambda} \Phi_{\mu\nu\lambda}(\ell\kappa;\ell^{\prime}\kappa^{\prime};\ell^{\prime\prime}\kappa^{\prime\prime}) e_{\mu}(\kappa;\boldsymbol{q}j)e_{\nu}(\kappa^{\prime};\boldsymbol{q}^{\prime}j^{\prime})e_{\lambda}(\kappa^{\prime\prime};\boldsymbol{q}^{\prime\prime}j^{\prime\prime}) \; ,\end{split}$

which becomes zero unless $$\boldsymbol{q}+\boldsymbol{q}^{\prime}+\boldsymbol{q}^{\prime\prime}$$ is an integral multiple of $$\boldsymbol{G}=n_{1}\boldsymbol{b}_{1}+n_{2}\boldsymbol{b}_{2}+n_{3}\boldsymbol{b}_{3}$$. Phonon linewidth $$\Gamma_{\boldsymbol{q}j}$$, which is the imaginary part of the phonon self-energy, can be obtained by the analytic continuation to the real axis ($$i\omega_{m}\to \omega + i0^{+}$$) as

(4)$\begin{split} \Gamma^{\mathrm{anh}}_{\boldsymbol{q}j}(\omega) &= \frac{\pi}{2\hbar^{2}}\sum_{\boldsymbol{q}_{1},\boldsymbol{q}_{2}}\sum_{j_{1},j_{2}} |V^{(3)}_{-\boldsymbol{q}j,\boldsymbol{q}_{1}j_{1},\boldsymbol{q}_{2}j_{2}}|^{2} \notag \\ & \times \left[ -(n_{1}+n_{2} + 1)\delta{(\omega + \omega_{1} + \omega_{2})} + (n_{1}+n_{2} + 1) \delta{(\omega - \omega_{1} - \omega_{2})} \right. \notag \\ & \left. \hspace{12mm} - (n_{1}-n_{2})\delta{(\omega - \omega_{1} + \omega_{2})} + (n_{1}-n_{2})\delta{(\omega + \omega_{1} - \omega_{2})} \right].\end{split}$

The computation of equation (4) is the most expensive part of the thermal conductivity calculations. Therefore, we employ the crystal symmetry to reduce the number of triplet pairs $$(\boldsymbol{q}j,\boldsymbol{q}^{\prime}j^{\prime},\boldsymbol{q}^{\prime\prime}j^{\prime\prime})$$ of $$V^{(3)}$$ to calculate. To disable the reduction, please set TRISYM = 0.

## 10.10. Isotope scattering¶

The effect of isotope scatterings can be considered by the mass perturbation approach proposed by S. Tamura [4] by the ISOTOPE-tag. The corresponding phonon linewidth is given by

$\Gamma_{\boldsymbol{q}j}^{\mathrm{iso}}(\omega)= \frac{\pi}{4N_{q}} \omega_{\boldsymbol{q}j}^{2}\sum_{\boldsymbol{q}_{1},j_{1}}\delta(\omega-\omega_{\boldsymbol{q}_{1}j_{1}}) \sum_{\kappa}g_{2}(\kappa)|\boldsymbol{e}^{*}(\kappa;\boldsymbol{q}_{1}\boldsymbol{j}_{1})\cdot\boldsymbol{e}(\kappa;\boldsymbol{q}\boldsymbol{j})|^{2},$

where $$g_{2}$$ is a dimensionless factor given by

$g_{2}(\kappa)=\sum_{i}f_{i}(\kappa)\left(1 - \frac{m_{i}(\kappa)}{M_{\kappa}}\right)^{2}.$

Here, $$f_{i}$$ is the fraction of the $$i$$th isotope of an element having mass $$m_i$$, and $$M_{\kappa}=\sum_{i}f_{i}m_{i}(\kappa)$$ is the average mass, respectively. The $$g_{2}$$ values should be provided by the ISOFACT-tag. The average mass $$M_{\kappa}$$ is substituted by the value specified in the MASS-tag.

## 10.11. Lattice thermal conductivity¶

The lattice thermal conductivity tensor $$\kappa_{\mathrm{ph}}^{\mu\nu}(T)$$ is estimated within the relaxation-time approximation as

$\kappa_{\mathrm{ph}}^{\mu\nu}(T) = \frac{1}{\Omega N_{q}} \sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}(T)v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\nu}\tau_{\boldsymbol{q}j}(T),$

where $$c_{\boldsymbol{q}j} = \hbar\omega_{\boldsymbol{q}j}\partial n_{\boldsymbol{q}j}/\partial T$$ and $$\tau_{\boldsymbol{q}j}(T)$$ is the phonon lifetime. The phonon lifetime is estimated using the Matthiessen’s rule as

$\tau_{\boldsymbol{q}j}^{-1}(T) = 2 (\Gamma_{\boldsymbol{q}j}^{\mathrm{anh}}(T) + \Gamma_{\boldsymbol{q}j}^{\mathrm{iso}}).$

The lattice thermal conductivity will be written to the file PREFIX.kl.

The spectra of the lattice thermal conductivity $$\kappa_{\mathrm{ph}}^{\mu\mu}(\omega)$$ can also be calculated by setting KAPPA_SPEC = 1 in the &analysis field. $$\kappa_{\mathrm{ph}}^{\mu\mu}(\omega)$$ is defined as

$\kappa_{\mathrm{ph}}^{\mu\mu}(\omega) = \frac{1}{\Omega N_{q}}\sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\mu}\tau_{\boldsymbol{q}j} \delta(\omega-\omega_{\boldsymbol{q}j}).$

If we integrate this quantity over $$\omega$$, we then obtain the bulk thermal conductivity, namely $$\kappa_{\mathrm{ph}}^{\mu\mu} = \int_{0}^{\infty} \kappa_{\mathrm{ph}}^{\mu\mu}(\omega) \; \mathrm{d}\omega$$.

## 10.12. Cumulative thermal conductivity¶

The accumulative lattice thermal conductivity $$\kappa_{\mathrm{ph,acc}}^{\mu\nu}(L)$$ is defined as

$\kappa_{\mathrm{ph,acc}}^{\mu\mu}(L) = \frac{1}{\Omega N_{q}} \sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\mu}\tau_{\boldsymbol{q}j}\Theta (L-|\boldsymbol{v}_{\boldsymbol{q}j}|\tau_{\boldsymbol{q}j}),$

where $$\Theta(x)$$ is the step function. This quantity can be calculated by using the script analyze_phonons.py with --calc cumulative flag. One can also use another definition for the accumulative thermal conductivity:

$\kappa_{\mathrm{ph,acc}}^{\mu\nu}(L) = \frac{1}{\Omega N_{q}} \sum_{\boldsymbol{q},j}c_{\boldsymbol{q}j}v_{\boldsymbol{q}j}^{\mu}v_{\boldsymbol{q}j}^{\nu}\tau_{\boldsymbol{q}j}\Theta (L-|v_{\boldsymbol{q}j}^{\mu}|\tau_{\boldsymbol{q}j}).$

In this case, the contribution to the total thermal conductivity is limited only from phonon modes whose mean-free-path along the $$\mu$$-direction is smaller than $$L$$. To calculate this, please use the --calc cumulative2 flag and specify the direction $$\mu$$ by the --direction option.

## 10.13. Delta function¶

In order to compute the phonon DOSs and the imaginary part of phonon self-energies, it is necessary to evaluate the Brillouin-zone integration containing Dirac’s delta function. For that purpose, we provide 3 options through the ISMEAR-tag.

When ISMEAR = 0, the delta function is replaced by the Lorentzian function as

$\delta(\omega) \approx \frac{1}{\pi}\frac{\epsilon^{2}}{\omega^{2}+\epsilon^{2}}.$

When ISMEAR = 1, the delta function is replaced by the Gaussian function as

$\delta(\omega) \approx \frac{1}{\sqrt{\pi}\epsilon}\exp{(-\omega^{2}/\epsilon^{2})},$

which decays faster than the Lorentzian function. For both cases, $$\epsilon$$ should be given by the EPSILON-tag, which must be chosen carefully to avoid any unscientific results. $$\epsilon$$ should be small enough to capture detailed phonon structures such as phonon DOS or energy conservation surface related to three-phonon processes, but it should be large enough to avoid unscientific oscillations. Choosing an appropriate value for $$\epsilon$$ is not a trivial task since it may depend on the phonon structure and the density of $$\boldsymbol{q}$$ points.

To avoid such issues, the program anphon employs the tetrahedron method [5] by default (ISMEAR = -1) for numerical evaluations of Brillouin zone integration containing $$\delta(\omega)$$. When the tetrahedron method is used, the EPSILON-tag is neglected. We recommend using the tetrahedron method whenever possible.

## 10.14. Self-consistent phonon (SCPH) calculation¶

The self-consistent phonon mode (MODE = SCPH) computes temperature-dependent phonon frequencies by solving the following equation self-consistently [6]:

(5)$V_{\boldsymbol{q}ij}^{[n]} = \omega_{\boldsymbol{q}i}^{2}\delta_{ij}+\frac{1}{2}\sum_{\boldsymbol{q}_{1},k,\ell}F_{\boldsymbol{q}\boldsymbol{q}_{1},ijk\ell}\mathcal{K}_{\boldsymbol{q}_{1},k\ell}^{[n-1]}.$

Here, $$\omega_{\boldsymbol{q}j}$$ is the harmonic phonon frequency and $$F_{\boldsymbol{q}\boldsymbol{q}_{1},ijk\ell} = \Phi(\boldsymbol{q}i;-\boldsymbol{q}j;\boldsymbol{q}_{1}k;-\boldsymbol{q}_{1}\ell)$$ is the reciprocal representation of fourth-order force constants. The updated phonon frequency in the $$n$$th iteration is obtained by diagonalizing the matrix $$V_{\boldsymbol{q}ij}^{[n]}$$ as

$\Lambda^{[n]}_{\boldsymbol{q}} = C^{[n]\dagger}_{\boldsymbol{q}}V^{[n]}_{\boldsymbol{q}}C^{[n]}_{\boldsymbol{q}},$

where $$\omega_{\boldsymbol{q}j}^{[n]} = (\Lambda^{[n]}_{\boldsymbol{q}jj})^{\frac{1}{2}}$$ and $$C^{[n]}_{\boldsymbol{q}}$$ is the unitary matrix that transforms the harmonic phonon eigenvectors into anharmonic ones as $$W^{[n]}_{\boldsymbol{q}} = W_{\boldsymbol{q}}C^{[n]}_{\boldsymbol{q}}$$. The matrix $$\mathcal{K}$$ in Eq. (5) is defined as

$\begin{split}\mathcal{K}_{\boldsymbol{q},ij}^{[n]} &= \alpha K_{\boldsymbol{q},ij}^{[n]} + (1-\alpha) K_{\boldsymbol{q},ij}^{[n-1]}, \\ K_{\boldsymbol{q},ij}^{[n]} &= \sum_{k} C_{\boldsymbol{q},ki}^{[n]} C_{\boldsymbol{q},kj}^{[n]*} \frac{\hbar\big[1+2n(\omega_{\boldsymbol{q}_{1}k}^{[n]})\big]}{2\omega_{\boldsymbol{q}_{1}k}^{[n]}}.\end{split}$

$$\alpha$$ is the mixing parameter, which can be changed via the MIXALPHA tag.

The SCPH equation is solved on the irreducible $$\boldsymbol{q}$$ grid defined by the KMESH_INTERPOLATE tag. The $$\boldsymbol{q}_{1}$$ grid in Eq. (5), given by the KMESH_SCPH tag, can be finer than the $$\boldsymbol{q}$$ grid. After the SCPH iteration converges, the code computes the anharmonic correction to the harmonic force constant $$\Delta D(\boldsymbol{r}(\ell))$$ as follows:

$\begin{split}&\Delta D(\boldsymbol{r}(\ell)) = \frac{1}{N_{q}}\sum_{\boldsymbol{q}} \Delta D(\boldsymbol{q}) e^{-i\boldsymbol{q}\cdot\boldsymbol{r}(\ell)}, \\ &\Delta D(\boldsymbol{q}) = D_{\mathrm{SCPH}}(\boldsymbol{q}) - D_{\mathrm{Harmonic}}(\boldsymbol{q}), \\ &D_{\mathrm{SCPH}}(\boldsymbol{q}) = W_{\boldsymbol{q}}C_{\boldsymbol{q}}^{[n]}\Lambda_{\boldsymbol{q}}^{[n]}C_{\boldsymbol{q}}^{[n]\dagger}W_{\boldsymbol{q}}^{\dagger}.\end{split}$

$$\Delta D(\boldsymbol{r}(\ell))$$ is saved in PREFIX.scph_fc2_correction.

The most computationally expensive part is the calculation of matrix elements of $$F_{\boldsymbol{q}\boldsymbol{q}_{1},ijk\ell}$$. When SELF_OFFDIAG = 0 (default), the code only computes the elements of $$F_{\boldsymbol{q}\boldsymbol{q}_{1},iikk}$$. Therefore, the computational complexity is $$\mathcal{O}(N_{q}^{\mathrm{irred.}}N_{q_{1}}m^{2})$$. When SELF_OFFDIAG = 1, the off-diagonal elements are also calculated, and the computational complexity is $$\mathcal{O}(N_{q}^{\mathrm{irred.}}N_{q_{1}}m^{4})$$.

 [1] K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 81, 3298 (1998).
 [2] Y. Wang et al., J. Phys.: Condens. Matter 22, 202201 (2010).
 [3] J. Hafner and M. Krajci, J. Phys.: Condens. Matter 5, 2489 (1993).
 [4] S. -I. Tamura, Phys. Rev. B 27, 858 (1983).
 [5] P. E. Blöchl, O. Jepsen, and O. K. Andersen, Phys. Rev. B 49, 1450555 (1994).
 [6] T. Tadano and S. Tsuneyuki, Phys. Rev. B 92, 054301 (2015).