Exposing the hidden influence of selection rules on phonon-phonon scattering by pressure and temperature tuning

Using ab initio calculations, we show that the hidden influence of selection rules on three-phonon scattering can be exposed through anomalous signatures in the pressure ($P$) and temperature ($T$) dependence of the thermal conductivities, $\kappa$, of certain compounds. Boron phosphide reveals such underlying behavior through an exceptionally sharp initial rise in $\kappa$ with increasing $P$, which may be the steepest of any material, and also a peak and decrease in $\kappa$ at high $P$. These features are in stark contrast to the measured behavior for many solids, and occur at experimentally accessible conditions. Similar anomalous behavior is predicted for silicon carbide and other related materials.

The thermal conductivity, κ, is a fundamental transport parameter that governs the efficiency of heat conduction through solids. In insulating crystals, the heat is transported by phonons, and intrinsic thermal resistance comes from phonon-phonon scattering processes, which arise from the anharmonicity of the interatomic bonding potential [1][2][3].
Lowest-order processes involving the mutual interaction among three phonons typically dominate this intrinsic resistance.
The influence of this selection rule on κ was long thought to be negligible since inter-branch phonon scattering processes dominate thermal transport, particularly at higher temperatures.
However, it was recently pointed out that, even around and beyond room temperature (RT; 300 K), this selection rule can be activated for high frequency phonons in materials where the three acoustic phonon branches come close to each other in some region of the Brillouin zone (BZ), thereby mimicking the behavior of phonon decay within a single branch [9,10,12]. The phase space available to a heat carrying, high frequency acoustic (A) phonon in such a region of the BZ, decaying into two others (AAA process) can then become quite small, thus resulting in its long lifetime and large contributions to κ.
Perhaps the most interesting example of the impact of this AAA selection rule is cubic boron arsenide (BAs). First principles calculations [9] predicted that unusually weak AAA scattering, resulting from activation of the AAA selection rule, contributes to a κ for BAs that should be comparable to that of diamond, the crystal having the highest measured κ of all materials at RT. While inclusion of higher-order four-phonon scattering was found to reduce the BAs κ [17,18], it still achieved a predicted RT κ value of around 1300 Wm −1 K −1 , by far the highest value of any naturally occurring material, behind only the carbon crystals (e.g. diamond and graphite). These theoretical findings have been confirmed by measurements [18][19][20], adding strong support to the impact of selection rules and corresponding reductions in the three-phonon scattering phase space on κ.
Recently, we have identified several other compounds, such as cubic boron phosphide (BP) and silicon carbide (SiC), where the AAA selection rule is activated, giving anomalously weak AAA scattering channels [12]. However, unlike in BAs, for these compounds, the weak AAA scattering is masked by another strong three-phonon scattering channel involving two acoustic phonons and an optic (O) phonon (AAO scattering channel), thus drastically reducing its impact on κ. In this letter, we show that the otherwise-hidden weak AAA scattering in BP and SiC can be exposed by hydrostatic pressure (P ) and temperature (T ) tuning of κ. We demonstrate that the evolving interplay between AAA and AAO scattering channels with varying P and T gives rise to two remarkable features in the κ of BP around and above RT -an unusually sharp rise in κ as P increases from ambient pressure and, a peak and subsequent drop in κ at high P . These features are contrary to the measured κ(P , T ) of many other insulating crystals that instead show a roughly linear increase with P far away from phase transitions [21][22][23][24][25][26][27][28], consistent with simple theories [26,29].
FIG. 1. RT κ of BP with naturally occurring B isotope concentration, isotopically pure 11 B, and isotopically pure 10 B as a function of P , with (solid) and without (dashed) fourphonon scattering. All curves show a peak and drop in κ with increasing P , unlike many other materials.
The κ(P , T ) of BP is calculated using a recently developed unified ab initio theoretical framework, which combines the use of density functional theory to obtain phonon dispersions and phonon scattering rates, with a numerical solution of the phonon Boltzmann transport equation [30]. The theory has no adjustable parameters, and it has demonstrated good agreement with the measured κ of several materials over a broad temperature and pressure range [12,18,[30][31][32], including BP and SiC. Both three-phonon and higher-order four-phonon scattering processes are included in the calculations along with phonon scattering by mass disorder arising from the natural isotope mixture on the constituent atoms [33]. The details of this first principles approach have been published in Ref. [30] and are summarized in the Supplemental Material (SM) section S1 [34]. We present the results for BP here, while those for SiC can be found in the SM section S4 [34]. Figure 1 shows the calculated RT κ(P ) of BP, with and without four-phonon scattering, for three cases: BP with (1) naturally occurring B isotope mix ( Nat BP), (2) B atoms isotopically enriched to 100 % 11 B ( 11 BP) and (3) B atoms isotopically enriched to 100 % 10 B ( 10 BP). All curves peak at around 55-60 GPa and subsequently decrease with increasing P , contrary to the typically measured linearly increasing behavior for many other materials. In addition, the κ(P ) curves rise rapidly with increasing P , and there is a significant increase in the slopes, dκ/dP , going from Nat BP to 11 BP to 10 BP. The peak values of κ are extremely large, around 1150 Wm −1 K −1 ( Nat BP), 1400 Wm −1 K −1 ( 11 BP) and 1600 Wm −1 K −1 ( 10 BP), about 2.5 times higher than their values at ambient pressure. It is important to note that, the effect of four-phonon scattering on κ is relatively weak at all pressures in BP. Thus, the non-monotonic κ(P ) behavior for BP is unrelated to that predicted for BAs [31], as discussed below.
To understand these pressure-dependencies of κ, we first examine the phonon-phonon scattering processes in BP at ambient pressure. The light atoms and stiff bonding of BP result in a high κ of around 500-600 Wm −1 K −1 at RT and ambient pressure [32,35,36]. BP crystallizes in the zinc blende structure and has two atoms in each unit cell, so the phonon dispersions consist of three A and three O phonon branches. Only three types of three-phonon scattering processes can occur, which involve combinations of A and O phonons: AAA, AAO, and AOO. Energy conservation forbids all other processes (e.g. OOO) [12]. Figure 2 (b) shows the RT AAA, AAO and AOO scattering rates for BP at various pressures, along with those for four-phonon scattering. The sharp dip in the AAA scattering rates at P =0 in Fig. 2 (b) results from the AAA selection rule driven by the close proximity of the three acoustic phonon branches in a small region of the BZ, as seen in Fig. 2 (a), along Γ → K direction. Similar dips driven by the same selection rule occur in other cubic compounds such as BAs, BSb, SiC, diamond, c-BN, and GaN [10,12] and certain transition metal carbides [37]. Some examples are shown in Fig. S1 in the SM. Apart from the close proximity of the acoustic phonons, other features in the phonon dispersions of BP and SiC activate additional selection rules. The mass difference between constituent atoms of these compound semiconductors causes a frequency gap, ∆ω A−O , between acoustic and optic phonons (A-O gap). Energy conservation mandates that only those acoustic phonons with frequencies greater than ∆ω A−O can participate in AAO scattering processes [9]. Similarly, energy conservation restricts the participation of acoustic phonons in AOO scattering processes to those with frequencies less than the optic phonon bandwidth, ∆ω O [11,12]. These restrictions are illustrated in Fig. S2 in the SM.
As shown in Fig. 2 (b), the AAO scattering rates in BP dominate over almost the entire frequency region of the dip in AAA scattering rates. The relatively small mass difference between boron and phosphorus atoms (M P /M B = 2.6) gives a ∆ω A−O that is not large enough to freeze out the AAO scattering channel for acoustic phonons above ∼6 THz. This behavior is in stark contrast to large mass ratio compounds like BAs (M As /M B = 7), where ∆ω A−O is large enough to almost completely remove AAO processes for acoustic phonons, thereby fully exposing the sharp dip in the AAA scattering rates, and resulting in significantly enhanced three-phonon limited κ [9,38]. To illustrate the importance of the AAO processes in BP, we perform calculations by artificially turning them off, and find that the RT κ of 11 BP at ambient pressure jumps from 530 , including three-phonon and four-phonon (three-phonon only) scattering. These values are about twice the corresponding values in BAs.
In BP, application of hydrostatic pressure introduces three intertwined behaviors that affect κ: (i) the longitudinal acoustic (LA) and optic phonons shift to higher frequencies, which increases LA phonon group velocities and weakens the AAO scattering rates due to decreased optic phonon occupations.
Both of these changes contribute to the increasing κ with P found in many materials. In BP, two additional features are critical to explain the anomalous behavior seen in Fig. 1 -(ii) optic phonons stiffen faster than do acoustic phonons [see Fig. 2 (a)], which increases ∆ω A−O and shifts the onset of AAO processes to higher frequencies, thereby exposing increasingly large portions of the AAA dip. This behavior significantly increases the rate of rise in the κ(P ) of BP in the low P range, as discussed below; and (iii) with increasing P , the LA phonons stiffen, while the transverse acoustic (TA) phonons weakly soften. The resulting increased separation between LA and TA phonon branches gradually removes the impact of the AAA selection rule, and so, increases the AAA scattering rates, relative to those at ambient pressure [see Fig. 2 (b)]. This behavior acts to drive κ lower, as pointed out previously [31,39].
Peak and decrease in κ -At RT and low P , the influence of trends (i) and (ii) is stronger than trend (iii) resulting in an increasing κ(P ). Beyond 50 GPa, the continued shift of AAO processes towards higher phonon frequencies almost fully exposes the AAA dip, and the rising AAA scattering rates from trend (iii) eventually dominate the behavior causing κ to decrease with increasing P . As a result of this evolving interplay between AAA and AAO processes, the RT κ(P ) achieves a peak value at around 60 GPa, which is around 2.5 times that of its value at P = 0. Figure 3 (a) shows the κ(P ) for 11 BP scaled by its zero pressure value, κ 0 , for different T . Around and above RT, the curves for each T roughly overlap, indicating that the P and T dependences are separable. Such separability is also predicted from empirical theory [25,26] and is found in conventionally behaving compounds such as MgO (see Fig. S4 in the SM). That it also occurs in BP with its anomalous non-monotonic κ(P ), is a consequence of the independent changes in AAA and AAO scattering rates induced by changing T or P . Increasing T strengthens AAO processes relative to AAA processes, while increasing P shifts the AAO processes to higher frequencies relative to AAA processes. Above RT, AAO scattering rates in BP dominate in magnitude over AAA scattering rates until very high P [ Fig. 3 (b)], just as they do at RT [ Fig. 2 (b)]. Thus, the evolution of κ(P ) above RT is qualitatively the same as that at RT. However, for T well below RT, as AAO scattering rates weaken relative to AAA scattering rates [see Fig. 3 (b)], the rise in κ from trends (i) and (ii) is overcome by trend (iii) at lower P , resulting in a weaker P dependence to κ(P ), which eventually decreases below κ 0 , as shown in Fig. 3 (a) for T = 100 K. Similar behavior of κ(P ) also occurs in SiC (see SM section S4).
Previously, we predicted a non-monotonic pressure dependence of κ for BAs [31], which was connected to an evolving competition between three-phonon and four-phonon scattering processes with increasing P . This competition is responsible for a significant suppression of the measured BAs κ-values compared with the lowest-order three-phonon prediction [17,18]. In Process-wise three-phonon and total four-phonon scattering rates for 11 BP at 100 K and 1000 K, and at different pressures. stark contrast, however, the non-monotonic behavior of κ(P ) in BP arises from fundamentally different physics involving competition only among three-phonon scattering channels. This competition contributes to a significantly larger peak-κ value in BP than in BAs (2.5x in BP vs. 1.1x in BAs at RT), and a much sharper increase in the BP κ at low P (discussed next), thus making it more accessible for experimental validation. We note that the effect of four-phonon scattering on the κ(P ) of BP and SiC is dominated by the increasing strength of the AAAA scattering channel with increasing P (see SM sections S3 and S4), but it is generally weak compared with the dominant three-phonon scattering channels, and so does not cause qualitative differences in the observed κ(P ) trends [see Figs. 1 and 3 (a)].
Unusually sharp rise in κ(P ) above P=0 -The calculated slopes, dκ/dP , of the κ(P ) curves in Fig. 1 around ambient pressure are exceptionally large. Significant contributions to these large values come from the rapid increase in the intrinsic lifetimes of acoustic phonons whose frequencies lie just below the onset of the AAO scattering. With increasing P , this onset shifts to higher frequencies, exposing more of the AAA dip. Phonons in the newly exposed region of the AAA dip have reduced scattering rates and so give correspondingly larger contributions to κ. This behavior is illustrated in Fig. 4, which shows, κ(ν), the spectral contributions to κ as a function of phonon frequency, ν, for 10 BP and 11 BP at different pressures and RT. The effect is enhanced in 10 BP compared with 11 BP because the lighter 10 B atom shifts the onset of AAO scattering processes to higher frequencies even at P = 0 (Fig. S3). This difference helps explain the much larger calculated RT and ambient pressure κ of 10 BP (630 Wm −1 K −1 ) compared with that for 11 BP (530 Wm −1 K −1 ), consistent with measured results [32,35]. Wm −1 K −1 GPa −1 ( Nat BP), 19.1 Wm −1 K −1 GPa −1 ( 11 BP) and 24.1 Wm −1 K −1 GPa −1 ( 10 BP) at low P . Measurements and calculations of dκ/dP given in the literature are only for materials with much lower κ, such as those of interest for geophysical studies [23][24][25][26][27], e.g. MgO, as well as alkali halides [21,22,28]. For these materials, dκ/dP is less than 3 Wm −1 K −1 GPa −1 , far smaller than our calculated values for BP. Therefore, to put the BP values in better context, we have calculated the RT dκ/dP for diamond, whose κ is higher than that of BP. We find RT dκ/dP values of 21 Wm −1 K −1 GPa −1 and 17 Wm −1 K −1 GPa −1 for isotopically pure and naturally occurring diamond respectively. Remarkably, the calculated dκ/dP value for 10 BP is higher than that for diamond, even though the absolute RT κ is five times lower. The present findings suggest that the dκ/dP for 10 BP may be the highest value achievable for any material, a striking signature of the influence of the AAA selection rule.
In summary, through ab initio calculations for BP, we have shown that a hidden weak three-phonon scattering channel arising from a selection rule on anharmonic decay of acoustic phonons can be revealed through examination of the pressure dependence of the thermal conductivity, which displays anomalous features including an exceptionally fast rise and a peak with increasing pressure. We predict that such anomalous behavior will also occur in other materials similarly affected by the same selection rule, which we have demonstrated for cubic silicon carbide (SiC) in the SM (see Figs. S6-S9). Other recently identified phase space selection rules should also give rise to unusual pressure dependence of phonon lifetimes and κ in the affected materials [12]. The present work gives new insights into the impact of selection rules on phonon-phonon scattering and demonstrates novel behavior of phonon thermal transport in solids under pressure.
This work was supported by the Office of Naval Research under a Multidisciplinary University Research Initiative, Grant No. N00014-16-1-2436.

First principles theoretical approach
We solve the linearized Boltzmann transport equation for phonons for the non-equilibrium distribution function, n λ = n 0 λ + n 1 λ , resulting from a small applied temperature gradient, ∇T , where n 0 λ is the Bose distribution function, and n 1 λ is the small deviation from n 0 λ generated by ∇T . Here, λ ∼ (j, q) designates the phonon mode, where j is the phonon branch and q is the phonon wave vector. The linearized Boltzmann equation is: where v λ is the phonon velocity of the mode λ, and the right-hand side is the collision term, which includes threephonon scattering, four-phonon scattering and phonon scattering by the isotopic mass disorder. Harmonic interatomic force constants (IFCs), and anharmonic third-and fourth-order IFCs are determined using the Density Functional Theory (DFT) as implemented in Quantum Espresso. Expressing n 1 λ = n 0 λ n 0 λ + 1 F λ · (−∇T ) allows Eq. 1 to be recast into an integral equation for the function, F λ , which is solved numerically. The thermal conductivity is then obtained as: where k B is the Boltzmann constant, V is the crystal volume, and κ αβ is the thermal conductivity tensor for heat current flow along the Cartesian direction, α, resulting from a temperature gradient along the direction, β. For the cubic structures considered in the present work, the thermal conductivity tensor is diagonal: κ αβ = kδ αβ .
Further details of the computational approach including computation scheme for IFCs, expressions for the phonon scattering rates and implementation of the solution of the phonon Boltzmann equation are given in Ref. [30]. All the calculations in this study are performed using norm-conserving pseudopotentials with the local density approximation (LDA) for the exchange correlation. The converged parameters used in the first principles calculations in this work for BP and SiC, such as the energy cutoffs for the DFT calculations, q-grids used in the solution of the Boltzmann equation and the cutoffs for the harmonic, cubic and quartic IFCs, can be found in Ref. [12]. The pressure in our calculations is calculated by obtaining the derivative of the Helmholtz free energy, correct to fourth-order in anharmonicity (F 4 th −order ), with respect to crystal volume (V ) at each temperature (T ) using the expression: , where a is the crystal lattice constant and ∆a ∼ 0.05% of a. The complete expression for the Helmholtz free energy, correct to fourth-order in anharmonicity of the crystal potential, can be found in section 10 of the supplemental materials in Ref. [31].  The slightly stiffer optic phonons in 10 BP compared to 11 BP shifts AAO processes to higher frequencies, thereby exposing more of the AAA phase space dip, causing weaker AAO scattering rates and, as a result, higher κ contributions from a narrow frequency region (shown by the blue shaded region in (c)) in 10 BP.
Increasing strength of four-phonon scattering in BP with increasing pressure Figs. 1 and 3(a) in the main text show that, as P increases, inclusion of four-phonon scattering increasingly suppresses the κ(P ) of BP, and that, this suppression becomes even larger at higher T . This behavior is driven by strengthening of AAAA four-phonon scattering with increasing P , as shown in Fig. 9. Hydrostatic pressure stiffens the chemical bonding, resulting in increased interatomic force constants (IFCs). These increases in the fourth-order IFCs produce corresponding increases in the four-phonon scattering matrix elements acting to increase the scattering rates of all four-phonon channels. At the same time, the increases in second-order IFCs shift optic phonons to higher frequencies, thereby reducing their populations and reducing the scattering rates of all four-phonon processes involving optic phonons, such as AAAO, AAOO and AOOO processes. As a result, AAAA four-phonon scattering rates are increased preferentially over all other four-phonon scattering rates, which contain at least one optic phonon.  9. Process-wise four-phonon scattering rates for BP at T = 300 K at different pressures. AAAA scattering rates increase with increasing P and dominate at high P .
Hidden influence of the selection rules on the κ of SiC The results for the pressure and temperature-dependent phonon dispersions, three-phonon and four-phonon scattering rates, and κ for SiC are presented in this section. Figure 10 (a) shows the stiffening of the longitudinal acoustic and all optic phonon branches, and the weak softening of the transverse acoustic phonon branch of SiC with pressure. Figure 10 (b) shows the evolution of the process-wise separated three-phonon and total four-phonon scattering rates with pressure at room temperature. The two plots show the striking similarity of the trends found for SiC and for BP (in the main text). Similarly, application of hydrostatic pressure causes a peak and drop in κ in SiC (Fig. 11) as in BP (in the main text). Furthermore, the observed temperature-and pressure-dependence of the competition between AAA and AAO scattering rates in BP is also seen in SiC (Fig. 12). Finally, the dominant all-acoustic AAAA four-phonon scattering channel (see Fig. 13) drives the reduction in κ of SiC at high P , also similar to BP.
FIG. 12. Process-wise three-phonon and total four-phonon scattering rates for SiC at T = 100 K and T = 1000 K, and different pressures.
FIG. 13. Process-wise four-phonon scattering rates for SiC at T = 300 K and different pressures. AAAA scattering rates increase with increasing P and dominate at high P , qualitatively similar to the behavior in BP.