# 6.3. ANPHON: Theoretical background¶

## 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}}{V} \frac{(Z_{\kappa}^{*}\boldsymbol{q})_{\mu}(Z_{\kappa^{\prime}}^{*}\boldsymbol{q})_{\nu}}{\boldsymbol{q}\cdot\epsilon^{\infty}\boldsymbol{q}},$

where $$V$$ 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  or the mixed-space approach  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.

## 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.

## Mode effective charge¶

The mode effective charge is defined as

$Z^{*}_{j,\alpha} = \sum_{\kappa\beta}Z^{*}_{\kappa,\alpha\beta}\frac{e_{\beta}(\kappa;\boldsymbol{0}j)}{\sqrt{M_{\kappa}}}$

where $$Z^{*}_{\kappa,\alpha\beta}$$ is the atomic Born effective charge. To compute the mode effective charges of the zone-center modes, please set ZMODE = 1 in the &analysis field and supply BORNINFO.

## 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 are saved in the PREFIX.thermo file.

When the self-consistent phonon mode (MODE = SCPH) is selected, the anharmonic free-energy defined by the following equation will be calculated and saved in the PREFIX.scph_thermo file:

\begin{align*} F^{\mathrm{SCP}} &= \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] \\ & - \frac{1}{4N_{q}}\sum_{\boldsymbol{q},j}\left[ \Omega_{\boldsymbol{q}j}^{2} - (C_{\boldsymbol{q}}^{\dagger}\Lambda_{\boldsymbol{q}}^{(\mathrm{HA})}C_{\boldsymbol{q}})_{jj} \right] \times \frac{\hbar [1 + 2n_{\boldsymbol{q}j} ]}{2\Omega_{\boldsymbol{q}j}}. \end{align*}

Details of the derivation of the above expression can be found in Ref. .

## Mean square displacement¶

The displacement-displacement correlation function is given by

\begin{align} \left< u_{\mu}(0\kappa)u_{\nu}(\ell'\kappa') \right> & = \frac{1}{\sqrt{M_{\kappa}M_{\kappa'}}N_{q}}\sum_{\boldsymbol{q}j} \frac{\hbar (2n_{\boldsymbol{q}j}+1)}{2\omega_{\boldsymbol{q}j}} \mathrm{Re}\bigg[e_{\mu}(\kappa;\boldsymbol{q}j)e_{\nu}^{*}(\kappa';\boldsymbol{q}j) e^{-i\boldsymbol{q}\cdot\boldsymbol{r}(\ell')} \bigg], \end{align}

where $$n_{\boldsymbol{q}j} = 1/(e^{\hbar\omega_{\boldsymbol{q}j}/kT}-1)$$ is the Bose-Einstein distribution function. When UCORR = 1, the code prints the above correlation function in PREFIX.ucorr. The vector $$\boldsymbol{r}(\ell')$$ is the lattice translation vector to the $$\ell'$$ th unit cell, which can be specified by the SHIFT_UCORR tag. When PRINTMSD is turned on, the code print the mean square displacements, which are the diagonal components of the correlation function:

$\left< u_{\mu}^{2}(0\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).$

## 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.

## (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}$$ . 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.

## 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.

## 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}^{2}},$

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.

## 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}j}\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.

## Isotope scattering¶

The effect of isotope scatterings can be considered by the mass perturbation approach proposed by S. Tamura  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.

## Lattice thermal conductivity (Peierls term)¶

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}{V 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 $$V$$ is the unit cell volume, $$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 is saved in 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$$.

## 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}{V 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}{V 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.

## Coherent component of lattice thermal conductivity¶

The coherent components of lattice thermal conductivity (see Ref. ), which are associated with the band off-diagonal components of the harmonic heat-flux operator, is calculated as

$\begin{split}\kappa_{\mathrm{c}}^{\mu\nu}(T) = \frac{1}{V N_{q}} \sum_{\substack{\boldsymbol{q},jj'\\ j\neq j'}}\frac{c_{\boldsymbol{q}j}\omega_{\boldsymbol{q}j'} + c_{\boldsymbol{q}j'}\omega_{\boldsymbol{q}j}}{\omega_{\boldsymbol{q}j}+ \omega_{\boldsymbol{q}j'}} v_{\boldsymbol{q}jj'}^{\mu}v_{\boldsymbol{q}j'j}^{\nu} \frac{\Gamma_{\boldsymbol{q}j}+\Gamma_{\boldsymbol{q}j'}}{(\omega_{\boldsymbol{q}j}-\omega_{\boldsymbol{q}j'})^{2}+(\Gamma_{\boldsymbol{q}j}+\Gamma_{\boldsymbol{q}j'})^2},\end{split}$

where $$c_{\boldsymbol{q}j} = \hbar\omega_{\boldsymbol{q}j}\partial n_{\boldsymbol{q}j}/\partial T$$ and $$\Gamma_{\boldsymbol{q}j}$$ is the total phonon linewidth (half width) of phonon mode $$\boldsymbol{q}j$$.

$$\boldsymbol{v}_{\boldsymbol{q}jj'}$$ is a band off-diagonal generalization of the group velocity . When KAPPA_COHERENT = 1 | 2 the coherent component is calculated and saved in PREFIX.kl_coherent. When KAPPA_COHERENT = 2, all components of the coherent term before summation are saved in PREFIX.kc_elem.

## Delta function¶

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  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.

## Self-consistent phonon (SCPH) calculation¶

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

(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_dfc2.

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})$$.

  K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 81, 3298 (1998).
  Y. Wang et al., J. Phys.: Condens. Matter 22, 202201 (2010).
  J. Hafner and M. Krajci, J. Phys.: Condens. Matter 5, 2489 (1993).
  S. -I. Tamura, Phys. Rev. B 27, 858 (1983).
  P. E. Blöchl, O. Jepsen, and O. K. Andersen, Phys. Rev. B 49, 1450555 (1994).
  T. Tadano and S. Tsuneyuki, Phys. Rev. B 92, 054301 (2015).
  Y. Oba, T. Tadano, R. Akashi, and S. Tsuneyuki, Phys. Rev. Materials 3, 033601 (2019).
  M. Simoncelli, N. Marzari, and F. Mauri, Nat. Phys. 15, 809 (2019).
  P. B. Allen and J. L. Feldman, Phys. Rev. B 48, 12581 (1993).