Formalism for program anphon

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.

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.

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.

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

\[\begin{split}\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),\end{split}\]

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

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}\) [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.

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}} 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}},\]

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

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.

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

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.

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.


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