FourPhonon: An extension module to ShengBTE for computing four-phonon scattering rates and thermal conductivity

FourPhonon is built within ShengBTE framework and its execution is fully compatible with ShengBTE. Before explaining FourPhonon, we briefly review the basic ShengBTE workflow [12]. It uses second-order and third-order force constants as inputs, retrieves phonon properties from harmonic calculations, estimates the available three-phonon scattering phase space, and then calculates the three-phonon scattering rates and lattice thermal conductivity. Besides the $2^{nd}$- and $3^{rd}$-order force constants files, users need to provide another input file, named CONTROL, which contains all the user-specified settings and parameters, including crystal structural information, temperature, $\mathbf{q}$-mesh, broadening factor, etc. On the basis of this workflow, FourPhonon requires an additional input file with fourth-order force constants. Flag is also added to the CONTROL file to enable FourPhonon utilities, namely four_phonon. When combined with other inherent flags or settings, users can perform various tasks, gaining insight into four-phonon scatterings. We now illustrate the modified workflow in more details. A workflow chart is presented in Fig. 1.

Due to the great computational cost of four-phonon scattering in most materials, one may need to estimate first how significant four-phonon scattering is compared to three-phonon scattering in a certain system. One can use the available three-phonon and four-phonon phase spaces as a criterion. To calculate these values, users can turn on the flags onlyharmonic and four_phonon, which will count the number of four-phonon scattering processes so that one can have a quick estimation of the computational demand when going forward.

Another particular problem is one may need to determine whether to solve the BTE exactly with an iterative scheme beyond single mode relaxation time approximations (SMRTA, or RTA) to get converged thermal conductivity. Some theoretical works have indicated that this is closely related to the relative strength of phonon normal/Umklapp processes [25, 18]. Since the iterative scheme would require huge computational resources (both processors and memory), users can first perform calculations in SMRTA scheme by turning on four_phonon flag and disable convergence flag. SMRTA calculations usually take much less time to finish than iterative scheme and would output scattering rates from normal/Umklapp scattering events separately. If Umklapp processes are dominant, then the accuracy in thermal conductivity may be acceptable without further utilizing iterative scheme. If normal processes are relatively strong, one can apply the iterative scheme to BTE up to three-phonon scattering by turning on convergence flag, with four-phonon scattering rates treated at the RTA level. Presently, this package is not able to involve the four-phonon scattering phase space in the iteration scheme due to the exceptionally high memory demands. Further optimization could be made to include this functionality in the future.

The linearized phonon Boltzmann transport equation with three-phonon scattering has been extensively studied in the literature, and the formalism on four-phonon scattering has also been well documented in several preceding papers [13, 14, 18]. Here, we only introduce a unified formalism implemented in our program and make all the corresponding notations consistent with the one used in ShengBTE [12]. For a certain mode $\lambda$ , the linearized BTE is expressed as:

where $\mathbf{F_{\lambda}}$ is the linear coefficient in nonequilibrium phonon distribution function and $\tau^{0}_{\lambda}$ is the relaxation time for mode $\lambda$ under SMRTA. $\mathbf{\Delta_{\lambda}}$ works for iterative scheme and is a quantity showing the phonon population deviation from the SMRTA scheme. With the inclusion of four-phonon scattering, $\mathbf{\Delta_{\lambda}}$ and $\tau^{0}_{\lambda}$ (0 represents zeroth-order in iterations) are computed as:

where $N$ is the total grid of $\vec{q}$ points. $\xi_{\lambda\lambda^{\prime}}=\omega_{\lambda^{\prime}}/\omega_{\lambda^{\prime\prime}}$. The superscripts $(\pm)$ or $(\pm\pm)$ represent the three-phonon (3ph) and four-phonon (4ph) processes, namely $\vec{q}^{\prime\prime}=\vec{q}\pm\vec{q^{\prime}}+\mathbf{Q}$ and $\vec{q}^{\prime\prime\prime}=\vec{q}\pm\vec{q^{\prime}}\pm\vec{q^{\prime\prime}}+\mathbf{Q}$, respectively. $\mathbf{Q}$ is a reciprocal lattice vector with $\mathbf{Q}=0$ implying normal process. Note that Eq.(2) is the full expression for three- and four-phonon iterative scheme, whereas the current version of FourPhonon does not support such functionality and excludes terms in the right bracket in Eq.(2). $\Gamma$ denotes the scattering rates for different processes and are computed by scattering probability matrices [26, 13]:

where Eq.(4) is for three-phonon processes and Eq.(5) for the four-phonon processes, with $n^{0}_{\lambda}$ being the phonon Bose-Einstein distribution at equilibrium, $\omega_{\lambda}$ being the phonon frequency for a certain mode $\lambda$. Conservation of energy is enforced by the Dirac delta function $\delta$. In Eq.(4) and (5), the matrix elements $V$ are given by the Fourier transformation of force constants, or transition probability matrices:

where $i,j,k,l$ denote the atomic indices and $\alpha,\beta,\gamma,\theta$ denote the Cartesian dimensions $x,y,orz$. $\Phi_{ijk}^{\alpha\beta\gamma}$ and $\Phi_{ijkl}^{\alpha\beta\gamma\theta}$ are the third-order and fourth-order force constants, respectively. $e_{\alpha}^{\lambda}(i)$ is the eigenvector component for a mode. $\mathbf{r}_{j}$ is the position vector of the unit cell where $j$th atom lies, and $M_{j}$ is its mass.

The remaining formulae including isotope scatterings and the expression of thermal conductivity are the same as ShengBTE, so are not elaborated in this paper.

The potential energy of crystals can be written as a Taylor expansion with respect to the atomic displacement based on perturbation theory [26]

where $r_{i}^{\alpha}$ is the displacement of atom $i$ along the $\alpha(x,y,z)$ direction, $E_{0}$ is the potential energy when $r=0$, and the coefficients $\Phi$ in the $n$-th order term are the corresponding $n$-th order interatomic force constants (IFCs). Given by the translational symmetry, fourth-order IFCs depend only on the relative coordinates of atom quartette $(i,j,k,l)$ rather than the absolute coordinates of the four atoms. The IFCs are usually determined by either a finite-difference scheme on the atomic forces [12] or fitting the atomic displacement-force relations [27].

In Fourthorder.py we employ a real-space finite-difference method to calculate the fourth-order IFCs:

where $h$ is a small displacement from the equilibrium position, and $F_{l}^{\theta}$ is the $\theta$ component of the force on the $l$-th atom. The default value of h is 0.04 times the Bohr radius, which can be changed by editing the header of the Fourthorder.py script. Each element of $\Phi_{ijkl}^{\alpha\beta\gamma\theta}$ requires eight DFT calculations with different supercell configurations, and there are 81$n^{4}N^{3}$ elements a supercell with $N$ unit cells with $n$ basis atoms per unit cell. A typical fourth-order IFCs calculation requires supercells containing more than 100 atoms, meaning that tens of thousands of single DFT simulation would be needed. Since directly performing these calculations is impractical for most computing infrastructures, it is crucial to make use of the symmetry analysis to reduce the computational cost.

We begin by populating the force constant matrices with transposition symmetries, and in the case of fourth-order IFCs this permutation of indices can have 24 sets of equality constraints

Then, we apply a space group symmetry operation $\sum_{\alpha}{\bf T}^{\alpha^{\prime}\alpha}{\bf R}_{i}^{\alpha}+{\bf b}^{\alpha^{\prime}}={\bf R}_{{\bf T}_{b(i)}}^{\alpha^{\prime}}$ to the tensor, where $\bf{T}$ and $\bf{b}$ represent the point-group and translation operators respectively, and ${\bf{T}}_{b(i)}$ specifies the atom to which the $i$th atom is mapped under the corresponding operation. The fourth-order IFC tensor must satisfy the following relation:

Each symmetry operation enables to map the four atom indices $i,j,k,l$ into themselves or to a different set. The first class of operations enforces $m$ constraints on the set of 81 fourth-order IFCs related to an atomic quartette, and each constraint can be described by a linear equation. Using Gaussian elimination, the $m$ linear equations can be transformed into the following form:

where asterisks in the columns stand for arbitrary numbers, corresponding to the independent elements among the set of 81 IFCs. The remaining force constants are linear combinations of them. Likewise with thirdorder.py , Fourthorder.py relies on the spglib library for detecting all kinds of crystal symmetries.

The fourth-order IFC tensor also obey the acoustic sum rules (ASRs)

Considering transposition symmetries mentioned above, Eq.(13) must also be valid if the summation is performed over $i$, $j$, or $k$.

Whether these sums are exactly zero has a significant effect on the calculated phonon scattering rates at low frequencies near the ${\bf\Gamma}$ point. Numerical uncertainties usually lead to small violations in ASRs, leading to significant deviations in calculated scattering rates. To tackle this problem, we adopt the same method with the thirdorder.py [12], namely adding a small compensation to each independent non-zero IFC, where the compensation is chosen by minimizing the sum of squares of these compensations by introducing a Lagrange multiplier.

It should also be noted that the fourth-order IFCs convergence with respect to the cutoff interatomic distance has to be checked manually, by examining either the force constants directly or the calculated thermal conductivity values. The computational cost increases significantly with increasing cutoff interatomic distance. Fortunately, it has been found [14, 28, 29, 21, 30] that, for many three-dimensional crystals, including up to second nearest neighbors in the fourth-order IFCs can give a satisfactorily converged values of scattering rates and lattice thermal conductivity.

The delta function $\delta$ for energy conservation is approximated by using the adaptive Gaussian broadening method [31, 32]. The expression for four-phonon processes is similar to the three-phonon case. Energy spacing level is described by taking the derivative of energy difference respect to the third phonon involved $\frac{\partial\Delta\omega}{\partial\mathbf{q}^{\prime\prime}}$. This adaptive broadening then needs the group velocity of a certain mode $\mathbf{v}_{\lambda}$ and the spacing of sampling $\mathbf{q}$ points $|\Delta\vec{q}^{\prime\prime}|$. The formula is expressed as:

while the parameter $\sigma$ for different processes is calculated as:

The detailed derivations for this formula can be found in B. Similar to ShengBTE, users of this program can always adopt a smaller scalebroad to reduce the computational cost by including fewer four-phonon processes and it can often give reasonable results. Our implementation also addresses a technical problem. When $\sigma$ evaluated is very small, the calculated scattering rates may be some abnormal high values. This numerical error is corrected by setting the $\delta$ function as 1.

Being the most commonly used semiconductor, silicon has been extensively studied as a benchmark material for theoretical predictions. We use silicon to illustrate the common practice for calculating thermal conductivity with four-phonon scattering.

Force constants are all calculated by first principles using VASP package [33]. Using a $6\times 6\times 6$ supercell, we have considered up to the fifth-nearest neighboring atoms in the third-order IFCs and second-nearest neighboring atoms in the fourth-order IFCs. Calculation of these IFCs is nearly a standard process and is well documented in some helpful tutorials [11]. In the calculation of silicon, we have turned on the iterative flag for three-phonon scattering.

Regarding the convergence with respect to the $\vec{q}$-points, our tests on silicon show that (see Fig. 2) the inclusion of four-phonon scattering would only require a $17\times 17\times 17$ $\vec{q}$-points grid to produce converged thermal conductivity values. In contrast, the three-phonon scattering typically requires a $25\times 25\times 25$ $\vec{q}$-points grid to converge. This is beneficial given the fact that the calculation of four-phonon scattering is very expensive.

Details on phonon-phonon interactions in Si at room temperature are shown in Fig. 3. Consistent with the literature [14], normal three-phonon processes do not overwhelm Umklapp processes. And in the four-phonon scattering, Umklapp process is dominant. Thus, applying only RTA scheme to Si is acceptable in calculating thermal conductivity.

Four-phonon effect is known to be important in boron arsenide (BAs), a zinc-blende structure material promising for thermal management applications. When determining the intrinsic thermal transport process, four-phonon scattering is even stronger than three-phonon scatterings. It was found that the higher-order anharmonicity would result in a substantial reduction of thermal conductivity [14], and subsequent experimental studies have validated the predicted thermal conductivity values. Therefore, we use BAs as an example material to further illustrate FourPhonon workflow in calculating thermal conductivity.

The sets of IFCs are obtained from first-principles on a $4\times 4\times 4$ supercell. Following parameters used in literature, we firstly perform calculations on a $16\times 16\times 16$ $\vec{q}$-points grid with scalebroad=0.1 at 300 K. For illustrating purpose, we assume that at this stage users have no knowledge on how important four-phonon effect is in BAs. Instead of directly calculating scattering rates with potentially very high cost, one can first look at the weighted phase space as defined in Ref.[36]. This is done by turning on onlyharmonic and four_phonon flags in the CONTROL file. Phase space calculation is less expensive to perform since it does not involve anharmonic IFCs. The results are shown in Fig. 4. We also plot the results of Si for comparison. It is found that phase space of four-phonon scattering is quite large in BAs, especially in optical branches. In contrast in Si, four-phonon phase space is generally smaller than three-phonon phase space. This is consistent with the fact that four-phonon scattering is not significant in silicon at room temperature and only has minor effect in its thermal conductivity [14]. Thus, the calculation of phase space can provide a preliminary estimation if scattering needs to be considered for thermal transport study or not.

Figure 5 presents the phonon scattering rates in BAs at 300 K and 1000 K. The four-phonon scattering rates can be comparable to or even larger than three-phonon ones even at room temperature. The scattering rates calculation takes about 1680 CPU hours for a single temperature point.

Thermal conductivity is then obtained by exactly solving BTE in an iterative scheme involving three-phonon interactions. To ensure the accuracy of scalebroad, we also perform calculations on the same mesh with scalebroad=1. The thermal conductivity values in W/mK are 1221.9 (scalebroad=1) vs. 1266.9 (scalebroad=0.1), and the relative difference is within 4%. Considering the fact that four-phonon scattering calculations are very expensive at a larger broadening factor, we set scalebroad=0.1 in calculations at all temperature points.

After performing calculations from 200 K to 1000 K, we get the temperature-dependent thermal conductivity, as shown in Fig. 6. Consistent with results in Fig. 5, larger reductions are observed for higher temperatures when considering four-phonon scatterings. Our results are in good agreement with experimental data.

$\mathrm{LiCoO_{2}}$ is a representative material of layered lithium transition-metal oxides ($\mathrm{Li_{x}TMO_{2}}$). This family of materials store lithium ions in between the $\mathrm{TMO_{2}}$ layers, and have significant applications in rechargeable batteries due to their supreme electrical and mechanical properties. However, their thermal transport properties seem to be hindering the improvement of device performance [37, 30]. As revealed by Feng et al. [30], the theoretical upper limit of $\mathrm{LiCoO_{2}}$’s lattice thermal conductivity would be lowered 2-6 times when including higher-order lattice anharmonicity. Compared to the preceding examples, lithium transition-metal oxides have more complex structures and involve four atoms in the unit cell. In this part, we intend to calculate the thermal transport properties of $\mathrm{LiCoO_{2}}$ using our new computational tool, and decompose the scattering rates into different scattering channels (i.e., recombination, the redistribution and splitting processes).

A $3\times 3\times 3$ supercell is employed to obtain IFCs by first-principles (for detailed information see Ref. [30]) and the Brillouin zone is discretized by $10\times 10\times 10$ $\vec{q}$-points grid. For this material with complex lattice structures, we only use scalebroad=0.01 to perform four-phonon calculations, and the available four-phonon processes are already around $10^{10}$.

Results of phonon scattering rates are shown in Fig .7. Even at room temperature, four-phonon scattering rates of $\mathrm{LiCoO_{2}}$ are generally comparable to three-phonon scattering rates over the entire frequency range. Different from BAs, four-phonon scattering rates are not that significant in optical phonon branches. Decomposing the four-phonon interactions further into different scattering channels provides us with more information. One can observe that redistribution process($\lambda+\lambda^{\prime}\to\lambda^{\prime\prime}+\lambda^{\prime\prime\prime}$) contributes the most and Umklapp process dominates this specific scattering channel.

We then calculate the in-plane ($\kappa_{\parallel}$) and through-plane ($\kappa_{\perp}$) thermal conductivity and compare their values with literature. Such information is available in a file named BTE.kappa_tensor. At 300 K, our calculations yield 9.35 W/mK for $\kappa_{\parallel}$ and 1.39 W/mK for $\kappa_{\perp}$, which agree reasonably well with Ref.[30] that gives 9.7 W/mK and 1.4 W/mK, respectively.