The role of high-order anharmonicity and off-diagonal terms in thermal conductivity: a case study of multi-phase CsPbBr₃

Phonons are crucial for understanding thermal transport in semiconductors and dielectrics. In the traditional approach of first-principles calculations, important quantities such as the mode-Grüneisen parameter, thermal expansion, phonon group velocity, three-phonon lifetime, and linewidth can be obtained by using anharmonic lattice dynamics under the quasi-harmonic approximation through the linearized phonon Boltzmann transport equation Togo et al. (2015); Carrete et al. (2014); van Roekeghem et al. (2016).

However, the traditional approach, in which some important factors have been overlooked, is facing several challenges:

1) the higher-order interatomic interactions like the fourth-order have been ignored for a long time like in BAs in which the four-phonon scattering is responsible for around 40% suppression of $\kappa_{L}$ compared to that with only three-phonon interactions;

2) interatomic anharmonicity increases as temperature is raised. However, the perturbation approach is difficult to deal with highly anharmonic systems, such as cubic ABX₃ perovskites with imaginary frequencies of harmonic phonons Tadano and Tsuneyuki (2015); Zhao et al. (2021). Therefore, we should illustrate the significance of anharmonic phonon renormalization, by utilizing the self-consistent phonon (SCPH) with temperature-dependent frequencies in several systems Xia et al. (2020); Kang et al. (2019);

3) the off-diagonal terms in the heat flux operator representing the heat transfer through the tunneling of wave-like coherent phonons could provide a potential method to bridge the gap between the traditional Peierls-Boltzmann transport equation and the experimental measurement in ultralow thermal conductivity materials and glass-like materials Simoncelli et al. (2019). Items 1) and 2) are related but independent from 3).

CsPbBr₃ is a classical chalcogenide material and a promising candidate for thermoelectric application, garnering significant attention in recent years. However, current studies have mainly focused on the thermal transport properties of individual phase Simoncelli et al. (2019); Tadano and Saidi (2022), often considering only three-phonon scattering and calculating the lattice thermal conductivity with zero temperature phonon dispersion Hu et al. (2021); Lee et al. (2017).

The systematic investigation of all three phases of CsPbBr₃, considering both renormalization effects and four-phonon scattering, as well as the contribution of off-diagonal terms, remains largely unexplored. Previous studies have primarily focused on phase transitions, vibrational mechanisms, and dielectric properties Lanigan-Atkins et al. (2021); Svirskas et al. (2020). To the best of our knowledge, no systematic studies have been conducted to date that encompass all three phases while considering renormalization, four-phonon scattering, and off-diagonal contributions.

In this study, we conduct a systematic investigation of the influence of quartic anharmonicity on the lattice dynamics and thermal transport properties of the three distinct phases of CsPbBr₃. We employ recent advancements in first-principles simulations, incorporating: (i) efficient construction of high-order interatomic force constants (IFCs) from the CSLD method Zhou et al. (2014, 2019a, 2019b); (ii) rigorous calculations of temperature-dependent phonons through SCPH theory and higher-order multi-phonon scattering rates Tadano and Tsuneyuki (2015); Feng and Ruan (2016); (iii) evaluation of the lattice thermal conductivity $\kappa_{L}$ by using a unified theory that considers both diagonal terms from the standard Peierls contribution and off-diagonal terms from the coherent Wigner distribution Simoncelli et al. (2019); Isaeva et al. (2019).

For the Peierls-Boltzmann transport equation, the lattice thermal conductivity $\kappa_{p}$ can be calculated as

where $\hbar$, $k_{B}$, $T$, $V$, $N_{0}$ are the reduced Planck constant, Boltzmann constant, absolute temperature, primitive unit cell volume, and the total number of sampled phonon wave vectors in the first Brillouin zone, respectively. $n_{\lambda}$, $\omega_{\lambda}$, $v_{\lambda}$, and $\tau_{\lambda}$ are the equilibrium component of the phonon population, frequency, group velocity, and lifetime for the $\lambda$ mode (wave vector $q$ and branch index $s$), respectively. Except for $\tau_{\lambda}$, all the above parameters can be obtained from harmonic approximation (HA). Usually, $\tau_{\lambda}$ can be obtained from the perturbation theory by consideration of three-phonon scattering Debernardi et al. (1995); Gao et al. (2018).

The temperature-dependent phonon dispersion could be considered by the anharmonic phonon renormalization (APRN) at finite temperatures Li et al. (2010); Souvatzis et al. (2008); Errea et al. (2011, 2014); Xia (2018); Ravichandran and Broido (2020). Among various existing approaches, SCPH Tadano and Tsuneyuki (2015); Xia et al. (2020) approximation is one effective method that can rigorously account for the first-order correction of phonon frequencies from the quartic anharmonicity. It can better describe the soft phonon modes and strong anharmonicity. In brief, under the SCPH approximation, the temperature-dependent renormalized phonon frequency $\Omega_{\lambda}$ can be obtained from the following equation

where $\omega_{\lambda}$ is the original phonon frequency from the harmonic approximation. The scalar $I_{\lambda\lambda_{1}}$ can be obtained as,

in which $V^{(4)}$ is the fourth-order IFCs in the reciprocal representation. The phonon population $n_{\lambda}$ satisfies Bose-Einstein distribution as a function of temperature. Both Eq. (2) and Eq. (3) have parameters $I_{\lambda\lambda_{1}}$ and $\Omega_{\lambda}$ in common, and thus the SCPH equation can be solved iteratively. Note that $I_{\lambda\lambda_{1}}$ can be interpreted as the interaction between a pair of phonon modes, $\lambda$ and $\lambda_{1}$ including the temperature effects Tadano and Tsuneyuki (2015); Xia et al. (2020).

Moreover, if one considers the off-diagonal terms of the heat-flux operator, which depicts the tunneling of coherent phonons, an additional contribution of lattice thermal conductivity, $\kappa_{c}$, needs to be considered Semwal and Sharma (1972); Knauss and Wilson (1974); Srivastava and Prasad (1981). Usually, $\kappa_{c}$ is neglected in simple crystals because of well-separated phonon dispersions and slight broadening as a function of temperature. However, it could dominate in disordered and glass-like amorphous compounds where phonon and related group velocities cannot be clearly defined, and heat transfer is mediated by diffusons and locons Allen and Feldman (1989, 1993); Cepellotti and Marzari (2016).

Recent studies show that $\kappa_{c}$ is substantial for materials with ultralow thermal conductivity, such as Mn₄Si₇ with twisting phonons Chen et al. (2015), Ba_7.81Ge_40.67Au_5.33 clathrate Lory et al. (2017), and Tl₃VSe₄ Mukhopadhyay et al. (2018). Therefore, in all-inorganic halide perovskite CsPbBr₃, we incorporate $\kappa_{c}$ as follows,

where the phonon lifetime in Eq. (1) is substituted as $\Gamma^{s}_{\bm{q}}=1/\tau_{\lambda}$, including three-phonon (3ph) and four-phonon (4ph) scattering. The group velocity is replaced with a generalized form containing off-diagonal elements Simoncelli et al. (2019); Allen and Feldman (1993),

in which $\bm{e_{q}}^{s}$ and $D(\bm{{q}})$ are the polarization vector and the phonon dynamical matrix, respectively. When $s=s^{\prime}$, it stands for the phonon band diagonal terms, while $s\neq s^{\prime}$, corresponds to the off-diagonal terms.

Therefore, the total lattice thermal conductivity $\kappa_{L}=\kappa_{p}+\kappa_{c}$. Note that in order to compute the generalized group velocity correctly, we used the phase convention that accounts for atomic positions within its lattice point to construct the dynamical matrix, as adopted in earlier studies Auerbach and Allen (1984). The details of the calculation are shown within the Supplemental Material sup ; Kresse and Furthmüller (1996); Blöchl (1994); Kresse and Joubert (1999); Esfarjani and Stokes (2008); Togo and Tanaka (2015); Candes and Wakin (2008).

CsPbBr₃ is a typical phase change material of ABX₃ perovskite. One can discriminate the transition temperature by the dynamical instability of the appearance of soft acoustic phonons from anharmonic potential energy surfaces He et al. (2020); Krapivin et al. (2022). Specifically for CsPbBr₃, a second-order phase transition occurs at about 318 K, and the transition temperature may vary, up to 361 K, depending on different samples. At about 373 K, there is another first-order phase transition. The critical temperature was also found at a higher temperature of about 401 K due to different experimental conditions Natarajan and Prakash (1971); Malyshkin et al. (2020). Here, we have chosen an intermediate temperature by considering the above different experimental values Natarajan and Prakash (1971); Malyshkin et al. (2020). Accordingly, the estimated temperatures for the first-order and second-order phase transitions are 320 K and 400 K, respectively. The high-temperature phase remains a cubic symmetry. The temperature reduction induces symmetry-breaking, leading to anisotropic structures from the cubic to the tetragonal at mid-temperature and finally to the orthorhombic crystals at low temperature Lanigan-Atkins et al. (2021); Bechtel et al. (2019); Malyshkin et al. (2020); Natarajan and Prakash (1971). All three crystal structures are depicted in Fig. 1 (a-c).

Fig. 1(e) shows the effects of SCPH, 4ph, and $\kappa_{c}$ on the calculated lattice thermal conductivity of CsPbBr₃. In the following, we neglect the SCPH notation for simplicity. Different primitive cells are used to calculate the corresponding temperature range. Due to the different crystal symmetry, we find that $\kappa_{L}$ increases from the orthorhombic to the tetragonal and the cubic phases. Moreover, $\kappa_{p}^{3,4ph}$+$\kappa_{c}$ in each phase decreases as temperature increases because of the enhanced phonon scattering. Since P1 and P2 are anisotropic, we use the arithmetic mean value in the figure.

Compared with $\kappa_{p}^{3ph}$, $\kappa_{p}^{3,4ph}$ is smaller due to the additional 4ph scattering. Moreover, the gap between them ($\Delta=\kappa_{p}^{3ph}-\kappa_{p}^{3,4ph}$) is growing significantly from P1 to P2 and finally to P3 based on Eq. (1). For instance, $\Delta$ is 0.114, 0.168, and 0.361 W m^-1 K^-1 for temperature 300 K, 400 K, and 500 K, respectively. It is usually attributed to the different scaling laws of 4ph ($\tau_{4}^{-1}\sim T^{2}\omega^{4}$) and 3ph ($\tau_{3}^{-1}\sim T\omega^{2}$) scatterings in which $\tau$ is the relaxation time Feng et al. (2017). Therefore, 4ph scattering is more critical than 3ph scattering at high-temperature and $\Delta$ is proportional to the temperature.

Since the ABX₃ perovskite has ultralow $\kappa_{L}$ and off-diagonal terms contribute significantly Lee et al. (2017); Simoncelli et al. (2019), we include $\kappa_{c}$ calculation of CsPbBr₃ based on Eq. (4) and Eq. (5). At 300 K, the value of $\kappa_{c}$ of CsPbBr₃ for P1 phase is 0.158 W m^-1 K^-1. At 400 K and 500 K, $\kappa_{c}$ of CsPbBr₃ for both P2 and P3 phases are 0.138 and 0.084 W m^-1 K^-1, respectively. More details of the values can be found in Supplemental Material S3 sup .

Our results of $\kappa_{c}$ agree reasonably well with the one reported by Simoncelli, Marzari, and Mauri Simoncelli et al. (2019), while the minor deviations might come from the size of the supercell used in the calculation of harmonic phonon as well as additional effects arising from quartic anharmonicity. For the P3 phase, the value of ($\kappa_{p}^{3,4ph}+\kappa_{c}$) is 0.501 W m^-1 K^-1 at 500 K, which is quite close to the result of 0.50 W m^-1 K^-1 from the recent quasiparticle nonlinear theory (QP-NL) Tadano and Saidi (2022).

Next, we investigate the influence of anharmonic renormalization on phonon dispersion among orthorhombic, tetragonal, and cubic phases of CsPbBr₃. Phonon-phonon interaction and lattice anharmonicity are ascribable to the cubic, quartic, and even higher-order IFCs.

The phonon spectra at different temperatures in three phases are shown in Fig. 2. Unexpectedly, we notice that the acoustic and optical phonon branches become hardened as temperature increases for the low-temperature orthorhombic phase of CsPbBr₃, shown in Fig. 2(a). Nevertheless, it is indisputable for the high-temperature cubic phase that the acoustic branches stiffen, whereas the top three optical branches soften with increasing temperature, shown in Fig. 2(c). For the tetragonal phase of CsPbBr₃ in Fig. 2(b), the high-lying optical phonons is almost temperature independent. We also plot the off-diagonal term contribution for CsPbBr₃ heat transport of three phases, which can be found in the Supplemental Material sup .

We also analyze the frequency-resolved (dashed lines) and cumulative (solid lines) lattice thermal conductivity $\kappa_{L}$ at different temperatures for cubic CsPbBr₃, shown in Fig. 2(d-f). Since 4ph has proved to be additional scattering, $\kappa_{L}$ of 3ph (the upper one) is more significant than that of 4ph (the lower one).

Traditionally, acoustic phonons are the main heat carriers. However, in cubic CsPbBr₃, it is found that optical phonons, ranging from 13.0 to 16.0 meV, dominate the heat transport, no matter whether or not the 4ph scatterings processes are accounted for. Furthermore, 4ph scattering reduces $\kappa_{L}$ of the cubic CsPbBr₃ by almost 40$\%$ on top of 3ph.

It is noticed that phonons of various frequencies dominate $\kappa_{L}$ among different phases. Optical phonons above the frequency of 15.0 meV control the heat transport of 3ph and 4ph for the orthorhombic phase. In comparison, acoustic and optical phonons among the frequency of 3.0 meV and 6.0 meV are also important for 3ph transport. However, for the tetragonal phase, we can find that optical phonons among the frequency from 13.0 meV to 19.0 meV dominate the phonon transport of 3ph and 4ph as well.

Different phonon-temperature tendencies is an interesting phenomenon that has yet to be thoroughly investigated. Previous work only found consistently softened or hardened phonons Zheng et al. (2022); Zhao et al. (2021) in different materials. For the first time, we unveiled the multi-tendency of high-frequency phonon modes variation occurring in different phases of the same material. We reveal in the following that such an opposite tendencies of optical phonon modes of cubic CsPbBr₃ as a function of temperature is sourced from the interaction between the top three opticals and other phonon modes.

In order to understand the underlying physical mechanism of the optical branches with high frequencies of different CsPbBr₃ phases have opposite temperature dependence in their phonon spectra, we have systematically studied the strength of 4ph interaction matrix elements $I_{\lambda\lambda_{1}}$ that have been introduced in Eq. (2) and Eq. (3). We set $\lambda$ for the highest optical phonon mode as any mode of the three highest optical phonon branches (index of mode=13-15 for cubic and 58-60 for the orthorhombic phase) and change $\lambda_{1}$ from the lowest acoustic phonon (index of mode =1) to all other optical phonons, gradually scrutinizing the interaction between phonon population $n_{\lambda}$, quartic-anharmonicity $V^{(4)}(\lambda,-\lambda,\lambda_{1},-\lambda_{1})$, phonon frequencies $\Omega_{\lambda_{1}}$, and $\Omega_{\lambda}$, respectively. Here we use the imode parameter to label the index of the phonon branches.

Interestingly, $I_{\lambda\lambda_{1}}$ is mainly positive for the orthorhombic phase, as is shown in Fig. 3. It can be both positive and negative for the tetragonal phase, and they almost cancel with each other leading to a small net frequency change. In contrast, it is found that the highest three optical branches have strong coupling with phonons in the low-frequency region, and the interaction $I_{\lambda\lambda_{1}}$ is even negative in the cubic CsPbBr₃, shown in Fig. 3. We further examine the wave vector position behind the negative $I_{\lambda\lambda_{1}}$ and discover that most negative $I_{\lambda\lambda_{1}}$ stem from the low-energy acoustic phonon modes, especially around $M$ and $R$ high-symmetry points.

In Supplemental Material S9-S11 sup , we also show the $I_{\lambda\lambda_{1}}$ between imode=3, 57 at 100 K for the P1 phase, imode=3, 27 at 350 K for the P2 phase, and imode=3, 12 at 500 K for the P3 phase, respectively. Since $I_{\lambda\lambda_{1}}$ can be either positive or negative, the renormalized phonon frequency $\Omega_{\lambda}$ as a function of temperature can either increase or decrease according to Eq. (2). Based on Eq. (3), only when $V^{(4)}(\lambda,-\lambda,\lambda_{1},-\lambda_{1})$ is negative for the three highest optical phonons of cubic phase, leading to a negative $I_{\lambda\lambda_{1}}$ and a reduced renormalization phonon frequency based on Eq. (2). On the contrary, $I_{\lambda\lambda_{1}}$ is positive as a function of temperature for the orthorhombic phase, and finally results in an increased $\Omega_{\lambda}$. Owing to the strong interaction between low-energy phonon modes around $M$ and $R$ and high-frequency optical phonons, the three highest optical phonon frequencies, with the temperature increasing, are softening for the cubic phase while hardening for the orthorhombic phase. More details can be found in the Supplemental Material sup .

To further understand the effects of anharmonic phonon renormalization and 4ph scattering on the thermal transport properties of CsPbBr₃, we continue to examine several parameters related to the lattice thermal conductivity, i.e., phonon phase space and scattering rates, respectively.

All available 3ph and 4ph scattering phase spaces need to satisfy the energy and quasi-momentum conservation simultaneously Lindsay and Broido (2008); Gao et al. (2018), shown in Fig. 4(a). The phase space of 3ph and 4ph scatterings increases as the temperature rises from 500 K to 800 K. Since the unit of phase space of 3ph and 4ph is different, one cannot compare them directly. Nevertheless, the larger phase space means more available scattering channels. The scattering strength in each accessible channel determines the final phonon relaxation time. Therefore, by including 4ph scattering, the lattice thermal conductivity is generally smaller than that of with only 3ph scattering.

The phonon scattering rates of the cubic phase are shown in Fig. 4(b). It displays that 4ph scattering has the same order of phonon scattering strength as that of 3ph. The scattering rates and phase space results for the orthorhombic and the tetragonal structures also show the same trends in Supplemental Material S5-S8 sup .

Fig. 4(c) displays the absorption and emission processes of the 3ph and 4ph as a function of frequency at 500 K, respectively. For the 3ph scattering, we consider the phonon splitting ($\lambda{\to}\lambda_{1}+\lambda_{2}$) and combination ($\lambda+\lambda_{1}{\to}\lambda_{2}$). For the 4ph situation, we count both phonon splitting ($\lambda{\to}\lambda_{1}+\lambda_{2}+\lambda_{3}$), and combination ($\lambda+\lambda_{1}+\lambda_{2}{\to}\lambda_{3}$), as well as redistribution ($\lambda+\lambda_{1}{\to}\lambda_{2}+\lambda_{3}$) processes. In the low-frequency region wherein acoustic modes dominate, 3ph combination processes are stronger than the splitting situation, while the redistribution processes of 4ph are the dominant ones. However, in the high-frequency region that is dominated by optical modes, the splitting process of 3ph becomes more important. For 4ph scattering, the splitting process also increases to a dominating portion and has the same order as the redistribution process.

Fig. 4(d) shows the temperature effect on $v^{2}$, where $v$ is the group velocity at different temperatures for the cubic phase. Interestingly, $v^{2}$ of 800 K is higher than that of 500 K for most frequencies but almost the same for the acoustic phonons. It stems from the renormalized phonon dispersions at finite temperatures. Besides, we distinguish phonon group velocity for the cases with and without considering SCPH in Supplemental Material S4 sup .

Besides, for the ABX₃ perovskites, we find that dynamical stability is in line with the thermodynamic stability Yang and Li (2022); Sun and Yin (2017). Previous work also used finite-temperature phonon dispersion of different phases of perovskites to predict the phase transition temperature Tadano and Saidi (2022). But for other materials, dynamical stability and thermodynamical stability have no direct connection.

In general, the higher the lattice constant, the weaker the interatomic interaction in materials, which will usually lead to a lower lattice thermal conductivity. On the one hand, the calculation of three-phonon scattering including thermal expansion in lower temperature phases (orthorhombic phase and tetragonal phase) is computationally prohibitive since they have lower symmetry and more atoms in the primitive cell compared with the cubic phase. On the other hand, based on the experimental investigations, the coefficient of thermal expansion for CsPbBr₃ is 3.8 $\times$ 10^-5 K^-1, 6.5 $\times$ 10^-5 K^-1, 2.6 $\times$ 10^-5 K^-1 for orthorhombic, tetragonal, and cubic phase, respectively Haeger et al. (2020). The effect of the lattice constant on the thermal conductivity might be neglected in a moderate temperature range. Therefore, in our calculation, we neglect the thermal expansion like previous works Simoncelli et al. (2019); Lee et al. (2017); Lanigan-Atkins et al. (2021). Most recently, an effective one-body Hamiltonian that well represents the quasiparticle-peak have been developed. In this method, the thermal expansion in the calculation is included Tadano and Saidi (2022).

In summary, our study reveals the significant contributions of four-phonon scattering and the off-diagonal terms of the heat flux operators in calculating the thermal conductivity in systems with harmonic phonons exhibiting imaginary frequencies and temperature renormalization.

Specifically, our investigation of CsPbBr₃ halide perovskites in orthorhombic, tetragonal, and cubic phases yields the following key findings:

(i) For materials with ultralow lattice thermal conductivity $\kappa_{L}$, the inclusion of high-order anharmonicity and off-diagonal terms bridge the gap between experimental observations and theoretical predictions;

(ii) The strong coupling between high-frequency optical phonons and overdamped acoustic phonons ($I_{\lambda\lambda_{1}}$), provides insights into the intriguing phonon renormalization phenomena observed in strongly anharmonic systems as a function of temperature;

(iii) Beyond the conventional phonon-phonon scattering perspective, phenomena such as electron-phonon coupling, polaron formation, and entropy in halide perovskites warrant further theoretical advancements Zhou et al. (2018); Lanigan-Atkins et al. (2021).

Our study that presents an effective approach to understand the ultralow $\kappa_{L}$ observed in halide perovskites, might inspire further experimental investigations exploring materials with glass-like thermal conductivity.

We acknowledge the support from the National Natural Science Foundation of China (No.12104356 and No.52250191). Z.G. acknowledges the support of China Postdoctoral Science Foundation (No.2022M712552), the Opening Project of Shanghai Key Laboratory of Special Artificial Microstructure Materials and Technology (No.Ammt2022B-1), and the Fundamental Research Funds for the Central Universities. We also acknowledge the support by HPC Platform, Xi’an Jiaotong University. Y.X. acknowledges Portland State University Lab Setup Fund.