Non-monotonic pressure dependence of the thermal conductivity of boron arsenide

Recent experiments demonstrate that boron arsenide (BAs) is a showcase material to study the role of higher-order four-phonon interactions in affecting heat conduction in semiconductors. Here we use first-principles calculations to identify a phenomenon in BAs and a related material - boron antimonide, that has never been predicted or experimentally observed for any other material: competing responses of three-phonon and four-phonon interactions to pressure rise cause a non-monotonic pressure dependence of thermal conductivity, κ, which first increases similar to most materials and then decreases. The resulting peak in κ shows a strong temperature dependence from rapid strengthening of four-phonon interactions relative to three-phonon processes with temperature. Our results reveal pressure as a knob to tune the interplay between the competing phonon scattering mechanisms in BAs and similar compounds, and provide clear experimental guidelines for observation in a readily accessible measurement regime.


Supplementary Note 1. Pressure-dependent properties of cBN and MgO
In this section, we present additional pressure-dependent calculations on cubic Boron Nitride (cBN) and Magnesium Oxide (MgO). Supplementary Figure 1 (a) shows the evolution of the phonon dispersions in cBN with pressure. Both acoustic and optic modes increase in frequency due to the application of hydrostatic pressure, similar to MgO (shown in fig. 1 of the main text). We also present the change in frequency of the transverse optic (TO) mode at the Γ-point as the lattice is compressed due to hydrostatic pressure, and the pressure-volume curve at 300 K for cBN in the Supplementary Figures 1 (b) and (c) respectively. Our calculations are in good agreement with experiments in literature [1][2][3]. The curve in Supplementary Figure 1 (c) was rigidly shifted to match the experimental volume at zero pressure only, and the corresponding volume was used as V 0 in Supplementary Figure 1 (b) (∼ 0.8% change in the lattice constant, only for these two plots).  [1]. The x-axis is scaled by the volume at zero pressure, V0. (c) Calculated pressure-volume curve at 300 K compared with experiments from Datchi et al. [1], Knittle et al. [2] and Solozhenkho et al. [3].  [4]. The x-axis is scaled by the volume at zero pressure, V0.
For MgO, we present our calculations for the pressure-volume curve at 300 K in Supplementary Figure 2, and we get good agreement with the experiments in Ref. [4]. The curve in Supplementary Figure 2 was rigidly shifted to match the experimental volume at zero pressure only (∼ 1.2% change in the lattice constant, only for this plot).
Supplementary Figure 3 (a) shows our calculated total three-phonon and process-wise four-phonon scattering rates of cBN at 300 K and at different pressures. There are two stark differences in the pressure dependence of the scattering rates between cBN and BAs (presented in the main text). First, the three-phonon scattering rates of heat carrying phonons (in the frequency range of 5-25 THz, see Supplementary Figure 3 (b)) are an order of magnitude stronger than the four-phonon scattering rates. Hence the effect of four-phonon scattering in cBN is weak at 300 K. Second, the three-phonon scattering rates decrease with increasing pressure. Similar trends for the pressure dependence of the three-phonon and four-phonon scattering rates are also observed for MgO (shown in Supplementary Figure 4). Since the three-phonon scattering rates, which dominate the total scattering rates, decrease with pressure, both the three-phonon limited thermal conductivity, κ (3) and the 3+4-phonon limited thermal conductivity, κ (3+4) of cBN and MgO increase with pressure unlike BAs, as shown in fig. 1 (b) and (c) of the main manuscript respectively.
Supplementary Figure 3. Pressure-dependent scattering rates and spectral κ for cBN. (a) Calculated total three-phonon and process-wise four-phonon scattering rates of the acoustic modes of cBN at 300 K and at different pressures. (b) Calculated spectral contributions to κ (3) and κ (3+4) of the acoustic modes of cBN at 300 K and at different pressures. The spectral contributions to both κ (3) and κ (3+4) increase with pressure, and the percentage reduction in κ due to four-phonon scattering decreases with increasing pressure.
Supplementary Figure 4. Pressure-dependent scattering rates and spectral κ for MgO. (a) Calculated total three-phonon and process-wise four-phonon scattering rates of the acoustic modes of MgO at 300 K and at different pressures. (b) Calculated spectral contributions to κ (3) and κ (3+4) of the acoustic modes of MgO at 300 K and at different pressures. Both scattering rates and spectral κ have similar responses to increasing pressure as in cBN.
Supplementary Note 2. Pressure dependence of the thermal conductivity of natural BAs In this section, we present κ (3) and κ (3+4) for BAs with naturally occurring isotopic composition of the constituent atoms (19.9% 10 B, 80.1% 11 B; As is isotopically pure). Supplementary Figures 5 (a) and (b) show that the pressure dependencies of κ (3) and κ (3+4) for natural BAs are qualitatively similar to the isotopically pure BAs results presented in the main text. Furthermore, Supplementary Figure 5 (c) shows that the pressure at which κ (3+4) of natural BAs peaks also shows a similar temperature dependence as the isotopically pure BAs. BAs, 3-phonon lines, showing the non-monotonic behavior of κ (3+4) of natural BAs with pressure at each temperature. Also shown is the shifting position of the κ (3+4) peak on the pressure-temperature surface.

Supplementary Note 3. Pressure-dependent properties of BAs using PBEsol
In this section, we present the pressure-dependent κ (3) and κ (3+4) for isotopically pure BAs using ultrasoft pseudopotentials with PBEsol exchange-correlation functionals (generalized gradient approximation -GGA type) from the GBRV pseudopotential library [5]. For these calculations, kinetic energy cutoffs of 50 Ry for the wavefunctions and 200 Ry for the charge density, and a 8 3 -shifted electronic k-grid produced a convergence of 6 × 10 −4 Ry per unit cell for the total energy and 0.15 kbar per unit cell for the total stresses. For the force-displacement calculations, 5X5X5 supercells were used, for which a Γ-shifted electronic k-grid produced a convergence of less than 10 −5 Ry/au for the forces. For the thermal conductivity calculations, the converged parameters were: 9 3 q-grid for the initial density functional perturbation theory calculation to obtain the bare harmonic interatomic force constants (IFCs) and 200 snapshots for the anharmonic IFCs for each point on the pressure-temperature grid, 17 3 q-grid for the solution of the Peierls-Boltzmann equation to obtain κ (3) and κ (3+4) using 11 nearest neighbors for the harmonic IFCs, 7 nearest neighbors for the cubic IFCs and 3 nearest neighbors for the quartic IFCs.
Supplementary Figure 6 (a) shows that the frequencies of the optic and the longitudinal acoustic (LA) modes of BAs, calculated with the PBEsol pseudopotentials, increase with pressure while the lower transverse acoustic (TA) branch undergoes weak softening, similar to results using the norm-conserving local density approximation (LDA) results presented in the main text. Supplementary Figure 6 (b) and (c) show that the monotonically decreasing κ (3) with pressure, the non-monotonic pressure dependence of κ (3+4) and the temperature-dependent shift in the pressure at which κ (3+4) peaks, are all observed in the calculated results with PBEsol, similar to the LDA calculations in the main text. Furthermore, the PBEsol results of the temperature-dependent shift in the peak-κ (3+4) pressure are shown in the pressure-temperature contour in Supplementary Figure 7 (a), and the opposing responses of the three-phonon and four-phonon scattering rates to pressure rise at 300 K and 750 K are shown in Supplementary Figure 7     average of 1140 ± 15% Wm −1 K −1 from the recent experiments on natural BAs [6][7][8][9]. Thus, the unique pressure and temperature-dependent features of the κ (3+4) of BAs are robust across the two pseudopotentials used in this study, and the observed small quantitative differences are consistent with that reported in the literature for the thermal and thermodynamic properties of BAs and other materials using these two exchange-correlation functionals [10][11][12][13].

Supplementary Note 4. Thermal conductivity of BAs without three-phonon scattering
Although four-phonon scattering strongly affects κ (3+4) of BAs even at room temperature, it is still significantly weaker than three-phonon scattering in most regions of the Brillouin zone. Only within the frequency range of ∼ 4-8 THz, the four-phonon scattering rates are comparable to the three-phonon scattering rates at ambient conditions (see figs. 3(a) and 4(b) in the main text), primarily because of the unusually weak three-phonon scattering in that frequency range, and not due to particularly strong four-phonon scattering in BAs. Note that away from the 4-8 THz region, the four-phonon scattering rates are much smaller than the corresponding three-phonon scattering rates. To further elucidate this important point, Supplementary Figure 10 shows the pressure dependence of three-phonon, four-phonon and 3+4-phonon limited thermal conductivity of BAs. Thermal conductivity including only four-phonon scattering (κ (4) ) in Supplementary Figure 10  in Supplementary Figure 10 (c) at room temperature, indicating that four-phonon scattering is still weak in BAs compared to other strongly anharmonic materials like sodium chloride [12] and lead telluride [14].
At temperatures higher than 500 K and at ambient pressure, four-phonon scattering rapidly strengthens and κ (4) becomes comparable to κ (3) . For example, at 750 K and ambient pressure, However, in this high temperature region, pressure rise strengthens the three-phonon scattering processes, causing κ (3) to decrease monotonically, but weakens the four-phonon processes, thus resulting in κ (4) [P > 0] > κ (4) [P = 0] even beyond 70 GPa, as shown in Supplementary Figure 10 (b). Also, as described in the main text, pressure rise strengthens AAAA four-phonon scattering channels but weakens AAOO and AAAO four-phonon scattering channels in BAs. These opposing responses to pressure rise, even among different four-phonon scattering channels, cause a non-monotonic pressure dependence of κ (4) too, as shown in Supplementary Figure 10 (b). Note that the origin of the non-monotonic pressure dependence of κ (4) (Supplementary Figure 10 (b)) is different from that of κ (3+4) presented in the main text and in Supplementary Figure 10 (c), which is caused by the opposing responses of three-phonon and four-phonon scattering channels to pressure rise.
However, quantitative conclusions about the pressure dependence of κ (3+4) cannot be drawn just from the plots of κ (3) and κ (4) (Supplementary Figures 10 (a) and (b) respectively), since the competition between three-phonon and four-phonon processes occurs within a small frequency range (∼ 4-8 THz), while κ (3) and κ (4) are averaged over the entire phonon spectrum. For example, at 200 K and ambient pressure, κ (4) > 10 κ (3) ; yet κ (3+4) ∼ 0.5 κ (3) is still affected by four-phonon scattering. Thus, the full spectral information about the rates of different phonon scattering mechanisms and their pressure dependencies must be considered (as in fig. 3 of the main text) to predict the pressure dependence of κ (3+4) of BAs.

BAs, 4-phonon only
Supplementary Figure 10. Relative strengths of three-phonon and four-phonon scattering in BAs. Three-phonon only (a), four-phonon only (b) and 3+4-phonon (c) limited thermal conductivity of isotopically pure BAs. Figures (a) and (c) are the same as in the main text ( fig. 2 (b) and (c) respectively). Thermal conductivity including only four-phonon scattering (κ (4) ) in (b) is much larger than κ (3) (a) and κ (3+4) (c) at room temperature, indicating that four-phonon scattering is still weak in BAs compared to other strongly anharmonic materials like sodium chloride [12] and lead telluride [14]. However, at temperatures higher than 500 K and at ambient pressure, four-phonon scattering rapidly strengthens and κ (4) becomes comparable to κ (3) .

Supplementary Note 5. Pressure dependence of the scattering phase space in BAs
In this section, we present the pressure dependence of the normalized scattering phase space of various three-phonon and four-phonon scattering channels in BAs at different temperatures. The normalized three-phonon scattering phase space is defined as [13,15]: where N p is the number of phonon polarizations, V BZ is the volume of the Brillouin zone, ω j (q) is the frequency of the phonon mode (q, j) and G is a reciprocal lattice vector that brings q ± q into the first Brillouin zone. Similarly the normalized four-phonon phase space can be defined as: The three-phonon (P 3 ) and four-phonon (P 4 ) phase spaces are normalized such that the energy-unrestricted phase space, where all δ-functions are set to unity, is equal to one second for both P 3 and P 4 . Supplementary Figures 11  (a) and (b) show the three-phonon and four-phonon scattering phase spaces respectively as functions of pressure at different temperatures, broken down by individual scattering channels for the acoustic modes of BAs. The phase spaces of all three-phonon scattering channels increase with pressure, while those of four-phonon scattering channels decrease with pressure. Furthermore, the AAA scattering channel has the largest three-phonon phase space, and the AAOO scattering channel has the largest four-phonon phase space. The AOOO four-phonon scattering channel has identically zero phase space for BAs at all temperatures and pressures, since the δ-functions cannot be satisfied for this combination of phonons; thus it is not shown in Supplementary Figure 11  Quantitative comparisons between P 3 and P 4 alone cannot be used to draw conclusions about the pressuredependence of κ (3+4) of BAs, since the corresponding scattering rates, which eventually determine κ (3+4) , also contain additional factors (Bose factors and the anharmonic scattering matrix elements), that (i) can have complex pressure dependencies different from those of the δ-functions in the phase space, and (ii) are, in fact, orders of magnitude different for three-phonon and four-phonon scattering channels. Furthermore, P 3 and P 4 are averaged over the entire phonon spectrum ( Supplementary Equations 1 and 2), while the competition between three-phonon and four-phonon processes occurs only within a small frequency range for the acoustic modes of BAs, as described in the main text. Thus, as mentioned in the main text, the phase space -being a harmonic quantity, can only present a part of the overall picture describing the pressure-driven competition between the three-phonon and four-phonon scattering strengths in BAs.

Supplementary Note 6. Temperature dependence of phonon scattering rates in BAs
In this section, we describe the temperature dependence of three-phonon and four-phonon scattering rates in BAs. Supplementary Figure 12 (a) and (b) show that four-phonon scattering rates have stronger temperature dependence than three-phonon processes at 0 GPa and 50 GPa respectively, consistent with the previous findings at ambient pressure [16]. This stronger temperature dependence of four-phonon scattering is the root cause of the temperature dependence of the non-monotonic behavior of κ (3+4) with pressure rise, as observed in the main text. For example, at both 0 GPa and 50 GPa, four-phonon scattering increases ∼ ten-fold when the temperature increases from room temperature to 1000 K, while three-phonon scattering only undergoes a two to four-fold enhancement (Supplementary Figure 12). Thus, the pressure-driven competition between three-phonon and four-phonon scattering channels becomes strongly temperature-dependent, as described in the main text.
Supplementary Figure 12. Temperature dependence of three-phonon and four-phonon scattering rates of BAs. Three-phonon and four-phonon scattering rates of different acoustic phonon modes in the Brillouin zone as functions of temperature at (a) 0 GPa and (b) 50 GPa. Each color represents a separate phonon mode. Three-phonon scattering rates show weaker temperature dependence at both low and high pressures compared to four-phonon scattering rates, consistent with previous findings at ambient pressure [16].

Supplementary Note 7. Pressure dependence of the thermal conductivity of BSb
In this section, we present additional results for the pressure-dependent properties of Boron Antimonide (BSb). Supplementary Figure 13 (a) shows that the frequencies of the optic and the LA modes of BSb increase with pressure while the lower TA branch undergoes weak softening, similar to BAs. For the isotopically pure BSb, κ (3) decreases with pressure and κ (3+4) shows a non-monotonic pressure dependence, as shown in Supplementary Figures 13 (b) and (c) respectively, also similar to BAs.
On the other hand, BSb with naturally occurring isotopic composition of the constituent atoms (19.9% 10 B, 80.1% 11 B, 57.2% 121 Sb and 42.79% 123 Sb) shows a qualitatively different pressure dependence of κ (3) compared to BAs and isotopically pure BSb, as shown in Supplementary Figure 14 (a). Below 500 K, even κ (3) of natural BSb shows a non-monotonic pressure dependence, which disappears above 500 K. This non-monotonic behavior at low temperature is also carried over to κ (3+4) , as shown in Supplementary Figure 14 (b). Above 500 K, κ (3+4) shows a similar non-monotonic behavior as BAs and isotopically pure BSb, with the percentage enhancement of peak κ (3+4) from 1000 K, κ 0 =311 Wm -1 K -1 750 K, κ 0 =424 Wm -1 K -1 500 K, κ 0 =663 Wm -1 K -1 300 K, κ 0 =1129 Wm -1 K -1 200 K, κ 0 =1687 Wm -1 K -1 the zero pressure value increasing with increasing temperature. The presence of a peak κ (3) at low temperature also results in the peak pressure position of κ (3+4) decreasing at first, then increasing with increasing temperature, as shown in Supplementary Figure 14   The principal reason for the non-monotonic pressure dependence of κ (3) at low temperature is the strong phononisotope scattering in natural BSb. Similar to BAs, the largest contribution to κ (3+4) of isotopically pure BSb comes from the frequency range where three-phonon scattering is the weakest (4-6 THz). Supplementary Figure 15 shows that, in this frequency range, phonon-isotope scattering rates are stronger than the three-phonon and four-phonon scattering rates in BSb at 200 K. Furthermore, the phonon-isotope scattering rates show weak reduction with increasing pressure at 200 K. Thus, the total phonon scattering rates, including three-phonon, four-phonon and phonon-isotope scattering, initially decrease with increase in pressure in the 4-6 THz frequency range at 200 K and between 0-20 GPa. This reduction in total scattering rates, along with the increasing group velocities of the LA mode and one of the TA modes with pressure, increases κ (3+4) of natural BSb at 200 K from 0-20 GPa.
As the pressure increases beyond 20 GPa, Supplementary Figure 15 shows that the three-phonon scattering rates increase and dominate over the phonon-isotope scattering rates, due to the weakening of the bunching effect described in the main text for BAs. Thus, κ (3+4) of natural BSb at high pressure is dominated by three-phonon scattering, resulting in a decreasing κ (3+4) with increasing pressure beyond 20 GPa at 200 K, similar to the isotopically pure BSb.
Supplementary Figure 15. Effect of isotope scattering on the thermal conductivity of natural BSb. Three-phonon (blue triangles), four-phonon (red stars) and phonon-isotope scattering rates (green squares) at 200 K (top three panels) and 750 K (bottom three panels) from zero to 40 GPa. From zero to 20 GPa, phonon-isotope scattering dominates total scattering rates at 200 K in the 4-6 THz frequency range, where the largest contribution to isotopically pure κ (3+4) comes from. Thus, the opposing responses of three-phonon and four-phonon scattering rates to pressure is not important at 200 K between 0 to ∼20 GPa, and so both natural κ (3) and natural κ (3+4) increase with pressure, similar to MgO. However, the pressure dependence of natural κ (3) and natural κ (3+4) is similar to that of isotopically pure BSb and BAs at 200 K beyond ∼20 GPa and at 750 K for all pressures, since the isotopic scattering rates are weaker than 3+4-phonon scattering rates under these conditions. At temperatures higher than 500 K, both three-phonon and four-phonon scattering rates are stronger than the temperature-independent phonon-isotope scattering (see Supplementary Figure 15 for scattering rates at 750 K). Therefore, the pressure-dependent behaviors of both κ (3) and κ (3+4) for natural BSb are similar to the isotopically pure BSb, as shown in Supplementary Figure 14. It is worth noting that, while the isotopic compositions of Boron in both natural BAs and natural BSb are the same, Antimony has a large isotopic composition while Arsenic is isotopically pure. Hence, the phonon-isotope scattering is significantly stronger in natural BSb and reduces κ (3+4) by 58% from the isotopically pure value at 200 K, while the corresponding reduction is only 20% in natural BAs.

Supplementary Note 8. Anharmonic phonon renormalization in BAs and BSb
In this section, we describe the effect of anharmonic phonon renormalization on the phonon and thermal transport properties of BAs and BSb. The renormalization procedure used in this work has been described in our prior work [12], and includes the effects of anharmonicity at finite temperature and the zero-point motion of atoms. Supplementary  Figure 16 (a) and (b) show that the effect of renormalization is small on the phonon dispersions of BAs at ambient pressure and 75 GPa at 300 K. From Supplementary Figure 16 (c) we find that, for BAs, the renormalization procedure causes 12% difference in κ (3) at ambient pressure throughout the temperature range considered, while the difference in κ (3+4) is 8%. For BSb, the renormalization procedure causes 15% difference in κ (3) , while the difference in κ (3+4) is 9%. Thus, the effect of phonon renormalization is weak in BAs and BSb, as expected for high κ, weaklyanharmonic materials. In this section, we describe the effect of anharmonic phonon renormalization and point defect scattering on κ (3+4) of MgO. We find from Supplementary Figure 17 (a) that the effect of phonon renormalization is weak in MgO and cannot explain the discrepancy with experiments from Ref. [17]. The observed discrepancy could be attributed to the presence of impurities in the sample used in Ref. [17]. For example, addition of Fe impurities (500 parts per million [ppm] as mass defects on the Mg site) to our calculation improves the qualitative and quantitative agreement with the experiments, as shown in Supplementary Figure 17 (b). Addition of Fe impurities lowers the κ (3+4) of natural MgO by ∼ 9.5% at 0 GPa and ∼ 15% at 60 GPa. Pressure-dependent weakening of three-phonon and four-phonon scattering rates (Supplementary Figure 4 (a)), and the pressure-driven increase in phonon group velocities cause a larger absolute effect of Fe impurity scattering on κ (3+4) of natural MgO at 60 GPa compared to that at 0 GPa in Supplementary Figure 17  In this section, we present the expressions for various phonon scattering probabilities that appear in the Peierls-Boltzmann equation (PBE) for phonon transport (equation 3 of the main manuscript) and the expression for the fourth-order anharmonic free energy used to determine the pressure. The three-phonon scattering probabilities (W (±) λλ1λ2 ) are given by, where −λ represents a phonon mode ((−q) s) when λ represents a phonon mode (qs) with wavevector q and polarization s, and Ψ λλ1λ2 ≡ Ψ qsq s q s and Ψ λλ1λ2λ3 ≡ Ψ qsq s q s q s are the third-order and fourth-order matrix elements given by: Ψ λλ1λ2 = Ψ qs,q1s1,q2s2 = ( /2) 3/2 1/N 1/2 0 [ω qs ω q1s1 ω q2s2 ] −1/2 × N P µνπ αβγ Ψ αβγ (0µ, N ν, P π) (M µ M ν M π ) −1/2 × e iq1·R(N ) e iq2·R(P ) × w α (qs, µ) w β (q 1 s 1 , ν) w γ (q 2 s 2 , π) and, Ψ λλ1λ2λ3 = Ψ qs,q1s1,q2s2,q3s3 = ( /2) 2 (1/N 0 ) [ω qs ω q1s1 ω q2s2 ω q3s3 ] −1/2 × N P Q µνπρ αβγη Ψ αβγη (0µ, N ν, P π, Qρ) (M µ M ν M π M ρ ) −1/2 × e iq1·R(N ) e iq2·R(P ) e iq3·R(Q) × w α (qs, µ) w β (q 1 s 1 , ν) w γ (q 2 s 2 , π) w η (q 3 s 3 , ρ) respectively. Here, Ψ αβγ and Ψ αβγη are the real-space third and fourth-order interatomic force constants (IFCs), α, β, γ and η are the Cartesian indices, N, P and Q are the indices for the lattice positions with lattice vectors R (N ) , R (P ) and R (Q) respectively, µ, ν, π and ρ are the indices for the basis atoms, [w (qs, µ)] α = w α (qs, µ) is the α th component of the eigenvector for the phonon mode (qs), M µ is the mass of the µ th basis atom, N 0 is the total number of q-points in the Brillouin zone and = h/2π with h being the Planck's constant. The phonon-isotope scattering probabilities (W iso λλ1 ) are given by [18], where g (σ) = 1/M 2 σ ζ f ζσ M ζσ −M σ 2 is a mass variance parameter with f ζσ and M ζσ being the concentration and mass of the ζ th isotope of the σ th atom respectively andM σ is the average mass of the σ th atom. The energy-conserving δ-functions in Supplementary Equations 5, 6 and 7 are computed using the analytical tetrahedron method [19].
The expression for the fourth-order anharmonic Helmholtz free energy (F 4 th −order ), which is used to determine the pressure in our calculations, is given by: [12,20]  (ω qs + ω q s + ω q s ) p + (2n qs n q s − n qs n q s + n q s ) (ω qs + ω q s − ω q s ) p + 2Ψ qs,−qs,q s Ψ q s ,−q s ,−q s n qs n q s + n qs + 1 4 (ω q s ) p where Ψ 0 , F H , F 3 and F 4 are the electronic, harmonic, third-order and fourth-order parts of the total anharmonic free energy.