Phonon-Phonon Interactions in Strongly Bonded Solids: Selection Rules and Higher-Order Processes

We show that the commonly used lowest-order theory of phonon-phonon interactions frequently fails to accurately describe the anharmonic phonon decay rates and thermal conductivity ($\kappa$), even among strongly bonded crystals. Applying a first principles theory that includes both the lowest-order three-phonon and the higher-order four-phonon processes to seventeen zinc blende semiconductors, we find that the lowest-order theory drastically overestimates the measured $\kappa$ for many of these materials, while inclusion of four-phonon scattering gives significantly improved agreement with measurements. We have identified new selection rules on three-phonon processes that help explain many of these failures in terms of anomalously weak anharmonic phonon decay rates predicted by the lowest-order theory competing with four-phonon processes. We also show that zinc blende compounds containing boron (B), carbon (C) or nitrogen (N) atoms have exceptionally weak four-phonon scattering, much weaker than in compounds that do not contain B, C or N atoms. This new understanding helps explain the ultrahigh $\kappa$ in several technologically important materials like cubic boron arsenide, boron phosphide and silicon carbide. At the same time, it not only makes the possibility of achieving high $\kappa$ in materials without B, C or N atoms unlikely, but it also suggests that it may be necessary to include four-phonon processes in many future studies. Our work gives new insights into the nature of anharmonic processes in solids and demonstrates the broad importance of higher-order phonon-phonon interactions in assessing the thermal properties of materials.


I. INTRODUCTION
Phonon-phonon interactions arise from the inherent anharmonicity of the interatomic bonds in solids. They determine fundamental properties of crystals such as the lattice thermal conductivity (κ) as well as the infra-red, Raman and neutron scattering cross sections [1][2][3]. They impact diverse phenomena including phonon drag [4], phonon bottleneck [5,6] and hydrodynamic thermal transport [7,8].
In his pioneering work, Peierls described intrinsic thermal resistance in semiconductors and insulators using the lowest-order theory of phonon-phonon interactions, which involves three phonons [1]. The complexity of this theory forced simple model approximations to be used until the past decade, when these simple models have been supplanted by first principles computational approaches to describe anharmonic phonon decay and heat conduction [9][10][11], still based on the lowest-order theory.
But what about weakly anharmonic materials i.e. those with relatively strong interatomic bonding? How well does the lowest-order theory work for such cases? Significant failures of the lowest-order theory have been found recently in cubic boron arsenide (BAs) [13,17], boron antimonide (BSb) [19] and aluminum antimonide (AlSb) [20]. But, the underlying reasons for the failures in BAs and BSb apply to very few materials, and a full explanation for the failure in AlSb has not been given in Ref. [20], as described below. At the same time, the lowest-order theory has been found to work well for diamond [13,14], gallium nitride (GaN) [20] and cubic boron nitride (BN) [18]. The computational cost of four-phonon calculations has thus far restricted consideration to only this small number of weakly anharmonic compounds. As a result, To date, a number of fundamental questions about the nature of anharmonic phonon-phonon processes in solids remain unanswered. Specifically: Are the failures of the lowest-order theory of phonon-phonon interactions rare or common among weakly anharmonic materials? Is there a general framework to understand the interplay between three-phonon and four-phonon processes that can drive the failures? How do four-phonon scattering strengths vary across such materials? Can new trends be identified to aid in the search for materials with high κ?
To address these questions, in this work, we examine the anharmonic phonon decay rates and thermal conductivities of seventeen different zinc blende (ZB) semiconductors with widely varying phonon properties using a first principles computational approach that includes both the lowest-order three-phonon and higherorder four-phonon interactions as well as interactions between phonons and the isotopic mass disorder found in most crystals. The findings for κ are summarized in Fig. 1 and confirmed by the available measured data for these materials (shown later). It is evident from the figure that failures of the lowest-order theory are remarkably common among ZB compounds, with nine of them showing a reduction in κ of at least 20% at room temperature, and at least 40% at 750 K. For many of the seventeen compounds, anharmonic decay rates of the optic phonons from four-phonon scattering are as large or even larger than those predicted by the lowest-order three-phonon theory.
Along with the previously identified compounds, BAs, BSb and AlSb, we find a particularly catastrophic failure of the lowest-order theory to occur for indium phosphide (InP). Yet, the properties responsible for the failures in BAs and BSb cannot explain the behavior in InP, as we show below. The same is true of the significant failures seen in gallium antimonide (GaSb), indium arsenide (InAs) and indium antimonide (InSb). At the same time, the lowest-order theory works quite well for many compounds containing boron, carbon or nitrogen atoms.
To fundamentally understand these findings, we first identify a set of selection rules that follow from energy and quasi-momentum conservation on the lowest-order three-phonon scattering and are dictated entirely by features in the phonon dispersions. To the extent that a material's phonon dispersions activate these selection rules, the phase space for the corresponding lowest-order three-phonon scattering channels are weakened or frozen out, which can result in unusually long lifetimes for both acoustic and optic phonons with correspondingly large contributions to κ, predicted in the lowest-order theory. While some of these selection rules have been identified previously [31][32][33], they cannot explain the behavior found in many of the studied compounds. We have identified new selection rules that are critical for developing a full understanding of the failures of the lowest-order theory in the zinc blende compounds.
Next, we show that in striking contrast, selection rules do not affect the four-phonon scattering rates, which generally increase with frequency for all seventeen materials. Thus, whenever the selection rules on the lowestorder processes drive weak three-phonon scattering, fourphonon scattering can become more important. Most remarkably, the materials naturally divide themselves 0.0 0.5 Pure ) for seventeen zinc blende compounds at 300 K (top panel) and 750 K (bottom panel). The large thermal conductivity reductions even at 300 K show that the lowest-order theory of phononphonon interactions is failing for many of these materials.
into two groups according to the strength of the fourphonon scattering processes. In the materials not containing boron, carbon or nitrogen (non-BCN compounds) the four-phonon scattering rates are of similar strength. For materials containing boron, carbon or nitrogen (BCN compounds) four-phonon scattering rates are up to an order of magnitude smaller than those in the non-BCN compounds, a difference we trace to the unusually strong bonding in the BCN compounds. As expanded on below, the large difference in the strength of four-phonon scattering in the non-BCN compounds compared with the BCN compounds has profound implications on the phonon lifetimes and κ.

II. SELECTION RULES ON PHONON-PHONON PROCESSES
We begin by describing six selection rules for the lowest-order phonon-phonon scattering and discuss their interplay with the higher-order interactions. Some of these selection rules have been discussed previously [31][32][33]. We focus on crystals with two different atoms in the unit cell, which encompasses a large number of semiconductors and insulators including the technologically important zinc blende compounds studied in the present work, which have inherently weak anharmonicity of the crystal bonds. This fact is demonstrated by the small thermal expansion and weak effects of anharmonic renormalization on their phonon frequencies (see Supplementary Materials [SM] section S1 for phonon dispersions and SM section S2 for details on thermal expansion for all seventeen compounds). At the same time, the breadth of the considered materials provides simple yet sufficiently different phonon dispersions that enable exploration of the effects of various selection rules, co-existing to different extents in different materials, on their anharmonic phonon decay rates and thermal conductivities.

A. Lowest-Order Processes
We divide the six phonon polarizations (branches) into three acoustic (A) and three optic (O) branches. There are four possible combinations of the lowest-order anharmonic processes that involve the decay of a phonon into two others or the coalescence of two phonons into a third. These are: (i) (AAA); (ii) (AAO); (iii) (AOO); (iv) (OOO). These processes must satisfy conservation of phonon energy and quasi-momentum [34]: Here, the three phonons participating in the scattering process have wave vectors q, q and q , branches, j, j and j , and frequencies, ω j (q), ω j (q ) and ω j (q ) respectively, and G is a reciprocal lattice vector to account for normal (G = 0) and Umklapp (G = 0) processes.
For each phonon mode, (j, q), and fixed j , j , the phase space for the lowest-order scattering is defined by the set of wavevectors -q and q , that satisfy Eqs. 1 and 2 [9,10,[35][36][37]. Features in the phonon dispersions that cause the phase space to vanish or significantly weaken define selection rules on these processes.
First, we note a selection rule on OOO processes: energy conservation freezes out OOO processes if the optic phonon bandwidth, ∆ω O , is smaller than the lowest optic phonon frequency. This is the case in all compounds studied in the present work, so we will not discuss the OOO selection rule further. We have identified five selection rules on AAA, AAO, and AOO processes that follow directly from Eqs. 1 and 2, combining those previously found [31][32][33] with two new ones. These selection rules divide into two types that complement each other. The first two selection rules (classified as Type 1 selection rules) follow directly from energy conservation (Eq. 1): 1. Type 1 selection rules 1. AOO#1 selection rule: AOO processes cannot occur involving acoustic phonons whose frequencies exceed the optic phonon bandwidth. This selection rule establishes a cut-off frequency (∆ω O ) in the acoustic phonon spectrum beyond which an acoustic phonon cannot participate in AOO processes. This cut-off has been pointed out previously [33].
2. AAO selection rule: An AAO process cannot occur if one of the participating acoustic phonons has a frequency smaller than the frequency gap between acoustic and optic phonons. This selection rule establishes a cut-off frequency for the acoustic phonons, ∆ω A−O -the A-O frequency gap, below which an acoustic phonon cannot participate in AAO processes. It also establishes a cut-off frequency for the optic phonons, 2∆ω A (∆ω A is the acoustic phonon bandwidth), above which optic phonons cannot participate in AAO processes. From this selection rule, it follows that for materials with A-O gaps larger than the highest acoustic phonon frequency, ∆ω A , AAO processes are completely forbidden [32].  Fig. 2 (e). If frequency windows are created by the AAO and AOO#1 selection rules, then the following three selection rules (classified as Type 2 selection rules) can drive anomalously weak three-phonon scattering: 2. Type 2 selection rules 3. AAA#1 selection rule: An anharmonic process of any order in which one phonon decays into a set of other phonons with higher phase velocity cannot occur. This selection rule was proven by Lax, Hu and Narayanamurti [31].
4. AAA#2 selection rule: As the group velocity of longitudinal acoustic (LA) phonons increases compared with transverse acoustic (TA) phonon velocities, the phase space of low frequency AAA processes decreases. To our knowledge, this selection rule has not been pointed out before. Unlike the other selection rules, this one does not completely forbid any particular three-phonon scattering channel except in the limit of infinite LA phonon velocities. However, the trend it expresses is critical in explaining the anomalously small acoustic threephonon scattering rates in several zinc blende compounds discussed in the next section. This selection rule is examined in more detail in the SM section S8(B).
5. AOO#2 selection rule: An AOO process is forbidden when both optic phonons derive from the same phonon branch, provided the group velocities of the optic phonons are smaller than those of the acoustic phonons. To our knowledge, this selection rule has not been pointed out before. Typically, the group velocities of optic phonons are smaller than those of the acoustic phonons that can participate in an AOO process while satisfying energy conservation. We have confirmed numerically that this selection rule is satisfied for all seventeen materials in this study. Exceptions can occur in materials with larger optic phonon group velocities, such as the lead chalcogenides. Further discussion is given in the SM section S8(A).

B. Minimal Influence of Selection Rules on Four-Phonon Processes
Ziman has suggested that selection rules on fourphonon processes are less restrictive than those on three-phonon processes [34]. Selection rules do exist for some four-phonon processes. For example, the general selection rule noted in Ref. [31] forbids A ↔ A + A + A processes when all phonons derive from the same branch. However, it does not forbid the corresponding processes. Our quantitative first principles calculations show that the AOAO processes dominate in all seventeen compounds studied in this work, with AAAA and OOOO processes also being important for acoustic and optic phonons respectively, and the influence of selection rules on four-phonon processes is minimal.
C. Interplay of lowest-order selection rules and four-phonon scattering The above selection rules impose strict constraints on the phase space for specific three-phonon processes, which are governed completely by a material's phonon dispersions a harmonic quantity. As described earlier, the selection rules complement each other: For materials whose constituent atoms have different masses (so as to create finite ∆ω A−O ) and relatively low ionicity, (so that ∆ω O is relatively small), AAO and AOO#1 selection rules create frequency windows where only AAA scattering can occur for acoustic phonons, and only AOO scattering can occur for optic phonons. Then, if features in the material's phonon dispersions activate AAA#1, AAA#2 or AOO#2 selection rules, regions of small total three-phonon phase space result. In principle, if hypothetical phonon dispersions existed with a collection of features described in the above selection rules, then the lowest-order theory would predict infinite intrinsic lifetimes and dissipationless transport for the affected phonon modes. For example, if the three acoustic phonon branches of hypothetical phonon dispersions coincided throughout the Brillouin zone (BZ), then the AAA#1 selection rule requires that the phase space for AAA scattering vanish. Then, if ∆ω A−O > ∆ω O , a frequency window would be created with zero phase space for lowestorder phonon-phonon scattering for acoustic phonons, according to AOO#1 and AAA#1 selection rules. This scenario is illustrated in Figs. 2 (a) and (b). Similarly, if all of the optic phonon branches of hypothetical phonon dispersions coincided throughout the BZ, the optic phonon velocities were much smaller than those of the acoustic phonons available that satisfy energy conservation, and ∆ω A−O > ∆ω A , then the optic phonons could not decay at all due to AAO and AOO#2 selection rules (no separate plot is shown for the optic phonon phase space, since it would be uniformly zero in this scenario). Although these scenarios do not occur in real materials, insofar as features in the real dispersions approach those in the hypothetical ones, selection rules will be activated that reduce the phase space for the lowest-order processes. For example, if acoustic branches have similar frequencies in some region of the BZ, then the phase space for some AAA processes involving phonons with frequencies in this region becomes small. If AAO and AOO#1 selection rules create a frequency window where the AAA processes are weak, the total three-phonon scattering phase space will be small, as shown in Fig. 2 (d). This is what happens in BAs and BSb, as has been explained previously [32]. Similarly, if phonons from different optic branches become nearly degenerate in some region of the BZ and their frequencies are larger than ∆ω A − ∆ω A−O , then the phase space for the corresponding AOO processes becomes weak. This is illustrated in Fig. 2 (e). We will show in the next section, that this scenario explains the large contributions to the thermal conductivity from optic phonons in BAs and AlSb. We note that the degree of activation of the type 2 selection rules is qualitative and cannot be fully assessed by visual inspection of phonon dispersions. Quantitative ab initio calculations are required to determine the impact of the selection rules on the anharmonic decay rates and thermal conductivity of a given material.
Higher-order four-phonon scattering rates are generally weaker than their lowest-order counterparts in weakly anharmonic materials at low-to-moderate temperatures. However, here we emphasize that when features in the phonon dispersions cause the abovedescribed selection rules to impose severe constraints on the phase space for the lowest-order phonon-phonon processes, the effects of four-phonon scattering can become quite important. In the next section and in the SM section S3, the calculated thermal conductivities of the studied seventeen different materials are presented and compared to measured data, and the failures and successes of the lowest-order theory are explained in terms of the selection rules and their interplay with the varying strengths of four-phonon scattering.
Three-phonon phase space Acoustic phonon frequency frequency window with zero phase space Three-phonon phase space Three-phonon phase space

III. RESULTS
The thermal conductivities of the seventeen compounds with zinc blende structure -GaAs, AlP, AlAs, GaP, GaSb, InAs, InSb, AlSb, InP, BN, BP, BAs, BSb, SiC, GaN, InN and AlN -are obtained as functions of temperature by solving the Boltzmann equation for phonon transport including three-phonon, four-phonon and phonon-isotope scattering using a recently-developed first principles approach [14]. A brief description of this approach is included in the Methods section and Appendix. Of the seventeen materials considered here, all but five (AlP, GaN, InN, AlN and BSb) have measured thermal conductivity data available for comparison. Six such comparisons are presented in Fig. 3 for BAs, AlSb and InP and Fig. 4 for GaSb, InAs and InSb, while the remaining (except BN) are given in the SM section S3. The case of BN has been extensively described in the main text and the supplementary information of our recent publication [18], so we do not present those plots in the SM. Calculated thermal conductivities are given for the following four conditions: lowest-order three-phonon scattering only (κ Pure ), three-phonon scattering and phonon scattering from the natural isotopic disorder (κ  SM Tables II and III. In many of the cases, κ (T ) calculated from the lowest-order theory significantly overestimates the measured κ (T ), while including four-phonon interactions brings the calculated κ (T ) curves into very good agreement with the measured data [17,18,[38][39][40][41][42][43][44][45][46][47][48][49]. As already seen in Fig. 1, four-phonon scattering suppresses the κ Pure for nine of the seventeen compounds by at least 20% at room temperature (RT, 300K) and at least 40% at 750K. The larger suppressions at 750K are explained by the more rapid increase of four-phonon scattering rates with temperature compared with three-phonon scattering rates [13,18]. Thus, it is clear that the lowest-order theory is failing for many of these compounds.
As explained in the previous section, activation of the Type 1 selection rules can create frequency windows where only AAA can occur for acoustic phonons and only AOO processes can occur for optic phonons. Pockets of small three-phonon scattering phase space result if Type 2 selection rules are simultaneously activated. The most significant effects from this combined activation of Type 1 and Type 2 selection rules occurs in the four compounds -BAs, BSb, AlSb and InP, those falling in the extreme classification in Fig. 1. Since different Type 2 selection rules are required to explain the extreme behavior in these materials, in the results discussed below materials are grouped according to which of the Type 2 selection rules (AAA#1, AAA#2 and AOO#2) are activated. We note that the degree of activation of Type 2 selection rules is somewhat qualitative but will be evident from the three-phonon phase space plots. The seven compounds -BAs, BSb, BP, BN, SiC, GaN, and AlN have acoustic branches that become bunched together, thereby activating the AAA#1 selection rule. The signature of this activation is the sharp dip in the AAA scattering phase space (see Figs. 5 (a), S6, S7 (a), S8 (a)). BAs and BSb have been discussed in depth in the context of ultra-high κ where the AAA#1, AAO and AOO#1 selection rules conspire to make three-phonon scattering unusually weak [32]. The large pnictide to B mass ratios of BAs and BSb (∼7 in BAs, ∼11 in BSb) give large A-O gaps that completely freeze out AAO scattering, and the narrow optic phonon bandwidths constrain AOO scattering of the acoustic phonons to low frequencies (see Figs. S1 (c) and (d) for phonon dispersions, Fig. 5 (a) and Fig.   S6 (c) for the three-phonon phase space for BAs and BSb respectively). As seen in these figures, these two selection rules expose a large frequency window of only AAA scattering processes. This combined with the activation of the AAA#1 selection rule gives a small phase space for the lowest-order processes and results in significantly weaker three-phonon scattering rates and large contributions to κ Pure (see Figs. 5(a), 5(d), 5(g) for the results for BAs, and Figs. S6 (c), S13 (c) and S16 (c) for those of BSb). Upon including the four-phonon scattering the RT κ Pure of BAs is sharply reduced from 3100 Wm −1 K −1 to around 1400 Wm −1 K −1 . In spite of this large reduction, we will show below that the four-phonon scattering in BAs is, in fact, much weaker for the acoustic phonons in the frequency range where the AAA three-phonon scattering rates are small, than in many other materials. We note that our calculations of κ including three-phonon, four-phonon and phonon-isotope scattering are in good agreement with the measured data [17,38,39] (Fig. 3 (a)).
The much smaller heavy-to-light atom mass ratios in BN, SiC, BP and AlN result in smaller A-O gaps in these materials. As a result, AAO scattering dominates in the same frequency region as the dips in the AAA scattering rates. The implications of this masking of weak AAA scattering are discussed below. The large κ Pure in these materials occurs for the conventional reasons: light atoms and stiff chemical bonds. Measured values [18,[50][51][52][53] are close to those predicted by the lowest-order theory. GaN has a relatively large Ga to N mass ratio of 5, which gives a large enough A-O gap to partially expose the dip in the AAA scattering rates (Fig. S8(a)). This contributes to a relatively large RT κ Pure ∼ 400 Wm −1 K −1 , consistent with prior calculations [54], but far lower than the RT κ Pure ∼ 3000 Wm −1 K −1 obtained for BAs.
In BP, BN, SiC, GaN, InN and AlN, inclusion of four-phonon scattering has only a weak effect on κ (3%, 5%, 9%, 4%, 13%, and 15% reductions in κ Pure at RT for SiC, BP, BN, AlN, InN, and GaN, respectively) as seen in Fig. S5 and in Ref. [18]. This is explained in part by the remarkably weak four-phonon scattering rates, which are found to be characteristic of the compounds with B, C and N atoms, as discussed in detail below.

B. Materials influenced by AAA#2 selection rule:
InP, InAs, InSb, GaSb Fig. 6 ranks the ratio of LA to TA phonon group velocities, averaged over the center quarter of the BZ around the Γ-point for all seventeen compounds considered in this work. Note that the choice of 1/4 th of the BZ is somewhat arbitrary and the same ordering would be obtained choosing e.g. 1/5 th of the BZ. The four compounds -InP, InAs, InSb, and GaSb show the largest differences between the averaged LA and TA phonon group velocities, so more strongly activate the AAA#2 selection rule. The most pronounced effect occurs in InP. The large In to P mass ratio of 3.7 in InP creates a large A-O gap, and the optic phonon bandwidth is small (see Fig. S4(b)). Thus, AAO and AOO#1 selection rules create a frequency window (1.5-3.5 THz) where only AAA scattering processes occur (red stars in Fig. 5  (c)). There, a significantly reduced phase space for AAA scattering results in the low frequency region (black oval in Fig. 5(c)). The striking effects of this restriction are seen in the spectral contributions to κ Pure at RT shown in Fig. 7 (i) where a large peak is seen around 1.5-2 THz.
In the 1.5-2 THz frequency region of weak AAA scattering, the RT four-phonon scattering rates are found to be comparable to or even exceed those of the three-phonon processes (Fig. 7 (f)). This results in a large suppression in by around 45% at 300 K, which increases to almost 70% at 750 K. These remarkably large reductions are supported by the measured data for InP [42,43] as shown in Fig. 4(c), which is in excellent agreement with the results including four-phonon scattering over a wide temperature range.
Large differences in the LA and TA velocities in InAs, InSb and GaSb (Fig. 6) also activate the AAA#2 selection rule, which gives the low frequency dips in the AAA phase space identified by the black ovals in Figs. 7(a)-(c). However, the smaller heavy-to-light atom mass ratios compared to InP result in smaller frequency windows, only in a small portion of the frequency range where the Pure for AlSb (h) also shows an anomalously large contribution from the optic phonons, which is completely suppressed by four-phonon scattering. dips in the AAA phase space occur are exposed. Since the dominant contributions to κ from the lowest-order theory occur in this frequency region, where four-phonon scattering is also strong (see Fig. 7), a significant suppression in κ Pure results. Four-phonon scattering suppresses the κ Pure for each compound in this group, ranging from around 25-30% at RT up to around 40-50% at 750 K. (b)-(c)). However, the relatively small heavy-to-light atom mass ratios in AlAs, AlP, BP, GaAs, GaSb, GaP, InAs, and InSb give no frequency window to expose the weak AOO scattering rates. Thus, the effects of the AOO#2 selection rule on optic phonon lifetimes are masked in these compounds. The most extreme effects occur in BAs and AlSb because the large A-O gaps create frequency windows where AAO scattering is frozen out, thus exposing the regions of weak AOO scattering. The unusually weak ionicity in BAs [55] results in a near degeneracy of the three optic branches near the Γ-point of the BZ (see SM Fig. S1(c)). Then, the actions of AAO and AOO#2 selection rules result in exceptionally small phase space for AOO scattering for the longitudinal optic (LO) and transverse optic (TO) phonons at the Γ-point ( Fig. 5(a), red box). As a result, large room temperature contributions to κ of around 130 Wm −1 K −1 are predicted from the lowest-order theory ( Fig. 5(g)), significantly exceeding the total κ of many common semiconductors. In contrast, BSb has a larger separation between the LO and TO phonons throughout the BZ (e.g., see SM Fig. S1(d)); thus the AOO#2 selection rule has little influence on the optic phonons, which contribute negligibly to κ.
As in BAs, the large A-O gap in AlSb (Fig. S2 (d)) resulting from the large Sb to Al mass ratio creates a frequency window with no AAO scattering and exposes the weak AOO scattering phase space from the AOO#2 selection rule (Fig. 5(b), red box). However, in stark contrast to BAs, this occurs near the degeneracy of LO and TO 2 optic phonon modes in regions of the BZ far from the zone center and away from high symmetry directions (see SM Fig. S10 (d)). In these regions, TO 2 + A ↔ LO processes become vanishingly small, while LO ↔ TO 1 + A and TO 2 ↔ TO 1 + A scattering is weak. The resulting small lowest-order optic phonon scattering rates (Fig. 5(e)) give anomalously large contributions to κ Pure of around 90 Wm −1 K −1 (Fig. 5(h)).
Four-phonon scattering in BAs and AlSb completely suppress the optic phonon lifetimes (Figs. 5 (d), (e), (g) and (h)), in striking contrast to predictions from the lowest-order theory. Fig. 5(e) shows that four-phonon scattering rates for the optic phonons are much larger than their three-phonon counterparts at 300 K for most of the optic phonon bandwidth. Surprisingly, this is true even at a low temperature of 100 K, where the influence of four-phonon scattering is expected to be weak even for strongly anharmonic materials [14]. The AOAO and OOOO four-phonon processes dominate the scattering as shown in the (SM section S9, Figs. S23(c) and S24 (d)). For AlSb, optic phonon contributions dominate the calculated RT κ Pure indicating that the lowest-order theory is failing drastically. In fact, our results indicate that there is very little hope for identifying materials with large contributions to κ coming from the optic phonons, since (a) at moderate-to-high temperatures, four-phonon scattering dominates any frequency region of weak three-phonon scattering for the optic phonons driven by the AAO and AOO#2 selection rules, thereby significantly lowering while (b) at low temperatures, where four-phonon scattering rates can get weaker than their three-phonon counterparts, the heat capacity of the optic phonons becomes vanishingly small owing to their high frequency scale (see SM section S7).
An interesting additional failure of the lowest-order theory for AlSb is its prediction of a large isotope effect (enhancement in κ upon isotopic enrichment), which is seen in the difference between the dashed black and dashed red curves in Fig. 3(b). This is quantified by P = 100 (κ Pure /κ Nat. − 1), which remains above 50% between 100 K and 300 K ignoring four-phonon interactions. In compounds where the heavy-to-light ratio of the masses of the two constituent atoms is large, phonon-isotope scattering is weak throughout the BZ if the vibrating atom is isotopically pure [32,56]. This is the case for AlSb (m Sb /m Al = 4.5) where the light, isotopically pure Al atoms dominate the vibrational component of the optic modes while the heavy Sb atoms, which have a large isotope mix (57% 121 Sb, 43% 123 Sb), hardly move. The resulting weak phonon-isotope scattering shows a striking contrast to the opposite behavior in BAs, where the heavy atom (As) is isotopically pure but the light atom (B) is not (compare phonon-isotope scattering rates in Figs. 5(d) and (e) ). Nevertheless, this weak phonon-isotope scattering strongly suppresses κ (3) Pure in AlSb because of the already weak AOO scattering. When four-phonon interactions are included, the P of AlSb drops to between 10% and 1% in the same temperature range. Note that the calculated κ (3+4) Nat. is in much better agreement with the measured κ data in Fig. 3(b) than κ   We note that large suppression of anharmonic optic phonon lifetimes and of thermal conductivity by fourphonon scattering has been predicted in prior work for BAs [13] and AlSb [20]. These works attributed the large lowest-order optic phonon lifetimes to freezing out of AAO scattering by the large A-O gaps. While it is true that AAO scattering is limited, this does not explain why AOO scattering is weak in BAs and AlSb. As discussed above BSb also has a large A-O gap that removes AAO processes but has much smaller contributions to κ (3) Pure than in either BAs or AlSb because its larger LO-TO splitting results in larger minimum AOO scattering rates. The new AOO selection rule identified in the present work provides the missing explanation for why AOO scattering becomes weak in both BAs and AlSb but not in BSb.

D. Other compounds: GaP, AlAs, GaAs, AlP
The heavy-to-light atom mass ratios in GaP (2.3) and AlAs (2.6) are large enough for AAO and AOO#1 selection rules to create frequency gaps where only AAA scattering occurs. Even though no Type 2 selection rules are activated in these regions, the lack of AAO and AOO scattering in these frequency windows give relatively weak three-phonon scattering, and correspondingly large contributions to κ (3) result. Four-phonon scattering suppresses κ (3) by around 25% at 300K and around 45% at 750K.
Relatively small reductions are found in the calculated values of κ (3) in AlP and GaAs upon inclusion of four-phonon scattering -around 10-15% at RT. The constituent elements of the two compounds have similar masses resulting in small A-O gaps, so AAO scattering extends throughout the range of the acoustic and the optic phonons (Figs. S8(c) for GaAs, S7(b) for AlP). Neither of the AAA selection rules are activated by features in the phonon dispersions for the two materials.
The results for AlP and GaAs highlight the fact that anomalous suppression of κ (3) by four-phonon scattering is unlikely to occur for strongly-bonded compounds with similar masses since no frequency window is opened by the Type 1 selection rules. A frequency window will also not occur in strongly polar materials having large LO-TO splitting and large optic phonon bandwidths. An example of this behavior is found is MgO [19], which has no frequency window opened by the selection rules. We have calculated the four-phonon scattering rates in MgO and find them to be comparable to those in the non-BCN zinc blende compounds. Since no selection rules are activated, three-phonon scattering is also strong. A reduction in κ (3) of about 15% at is found at RT, consistent with that found in GaAs and AlP.

IV. FOUR-PHONON SCATTERING AND ITS WEAKNESS IN COMPOUNDS WITH BORON, CARBON OR NITROGEN
Comparison of the four-phonon scattering rates for the seventeen materials studied shows several striking features. First, the four-phonon scattering rates generally increase with increasing frequency, being the strongest for the optic phonons. In particular, they are far larger than the small lowest-order optic phonon scattering rates in BAs and AlSb responsible for the erroneously large contributions to κ. As noted above, this finding suggests that it is not likely to be possible to get large contributions from optic phonons in any compound, and any findings to the contrary in the lowest-order phonon-phonon theory are likely incomplete and point simply to the need to consider four-phonon interactions to obtain accurate results. Nevertheless, the dominance of four-phonon anharmonic decay rates for the optic phonons compared with their three-phonon counterparts that we find in BAs, AlSb, InP, AlAs and BSb around Γ (see Figs. 5(d)-(f), S14(c) and S13(c) respectively) suggests that these signatures should be observable in e.g. temperature dependent Raman linewidth measurements.
Second, the four-phonon scattering rates for the nine non-BCN compounds i.e. those that do not contain B, C, or N atoms (GaAs, GaP, GaSb, AlSb, AlP, AlAs, InP, InAs, InSb) are remarkably similar in magnitude. An example of this is seen in Figs. 5 (e) and (f) for AlSb and InP, respectively. Comparisons for the other compounds can be made in Figs. 7 (d)-(f), and Figs. S14 (b)-(c) and S15 (b)-(c) in the SM. These similarities across materials are in striking contrast to the analogous lowest-order scattering rates. This highlights a fundamental difference between the lowest-order phononphonon interactions and those of higher order. The lowest-order processes are constrained by selection rules; features in the phonon dispersions that activate these selection rules can significantly reduce the lowest-order phonon-phonon scattering rates. In contrast, the phase space for four-phonon scattering shows remarkably little structure across materials, being immune to shifts in the A-O gap, or to the bunching together of the acoustic or optic branches. The four-phonon scattering rates of the acoustic phonons for all compounds are dominated by AOAO and AAAA processes, and for those, fourth-order bond anharmonicity and phonon frequency scales act in opposition, maintaining relative uniformity in the scattering strengths. This is discussed in more detail in the SM section S10. Based on the relatively uniform strength of the four-phonon scattering rates and lack of selection rules influencing them, it seems reasonable to expect that higher-order phonon-phonon interactions beyond four-phonon scattering (five-phonon and beyond) should play a negligible role in the thermal conductivity of this group of materials.
The most remarkable feature is that the four-phonon scattering rates for the BCN compounds, i.e. those containing B, C or N atoms (BAs, BSb, BP, BN, SiC, GaN, AlN and InN) are much smaller than those for the non-BCN compounds. This is clearly seen by comparing Figs. 5(d)-(f), which show that the four-phonon scattering rates in BAs are an order of magnitude smaller than those for InP and AlSb. Figs. S13, S14 and S15, and Fig. 3B in Ref. [18] show that the four-phonon scattering rates for BSb, BP, BN, SiC, GaN, AlN and InN are similarly small compared to their non-BCN counterparts. The four-phonon scattering rates at 750 K show the same trends as those identified at RT. We have traced this behavior to fundamental bonding differences between BCN and non-BCN compounds. BCN compounds have stronger bonding, reflected in their smaller lattice constants and higher phonon frequency scales compared to the non-BCN compounds. However, the fourth-order anharmonic component of the bonds in the BCN compounds tends to be weaker relative to the harmonic component, compared with those in the non-BCN compounds. The implications of this finding are discussed below. Interestingly, recent calculations of thermal transport in ZrC [21], a metal with rock salt structure, also show unusually weak four-phonon scattering rates and correspondingly small reductions in the phonon-phonon limited κ, consistent with our findings for the zinc blende BCN compounds.
Exploiting our finding that the AOAO processes are dominant in all the considered materials, we have developed a simple model that captures the weaker fourphonon scattering strength in the BCN compounds compared with the non-BCN materials. This model is discussed in the SM section 10, and findings are plotted in Fig. S27.

V. DISCUSSION
From the above findings we make the following additional observations: Significance of weak four-phonon scattering in BCN compounds: While the κ of BAs is significantly suppressed by four-phonon scattering, it still achieves an ultra-high RT value of κ (3+4) Nat. of around 1300 Wm −1 K −1 , confirmed by measurements [17,38,39], about two to three times larger than other high κ compounds such as SiC, BP, copper and silver. To test the importance of its unusually weak four-phonon processes, we artificially increased the RT four-phonon scattering rates in BAs by a factor of ten, making them comparable to those in the non-BCN compounds. The resulting RT κ of BAs is reduced by over 70% to 330 Wm −1 K −1 . Similarly increased four-phonon scattering rates in GaN and BP yield 50% and 30% decreases in the RT κ Pure values for these compounds. This highlights one of the fortuitous benefits that nature gives to the BCN compounds, without which their κ's could be much lower.
High κ in non-BCN compounds unlikely: The larger four-phonon scattering found in the non-BCN materials suggests that finding a material with high κ among them will be challenging. Indeed, note that of the eight compounds above, showing small deviations from the predictions of the lowest-order theory upon inclusion of four-phonon scattering, only two (GaAs and AlP) are non-BCN compounds. Conventional guidelines for achieving high κ, as described by Slack [57], require the combination of light atoms and exceptionally strong bonding, with diamond being the prototype. These criteria already rule out the non-BCN compounds. The alternative paradigm exemplified by BAs [32] in principle allows for the possibility of a high κ material composed of non-BCN constituent atoms.
Such a material would require phonon dispersions for which the selection rules give small three-phonon scattering phase space in a frequency region where the four-phonon scattering is weak. This does not occur in any of the studied non-BCN compounds, where instead, frequency windows with small total three-phonon phase space coincide with strong four-phonon scattering regions.
Masking of weak AAA and AOO scattering: The four compounds, BAs, BSb, AlSb and InP, showing the catastrophic failures of the lowest-order theory have the largest A-O gaps, which remove the AAO processes. This feature is clearly a critical requirement to achieve frequency windows in which anomalously large phonon lifetimes are possible. To highlight the importance of this requirement, consider BP and SiC compared with BAs. The AAA, AOO and four-phonon scattering rates in SiC, BP and BAs are similar. Thus, the reason SiC and BP have lower κ than BAs stems from the AAO scattering rates, which are much larger than the AAA scattering rates in the frequency range of the AAA minima. To illustrate how important these are, we have also calculated the κ of SiC and BP with AAO processes artificially removed. We find that the RT κ While it is unfortunately not possible to remove these processes in reality, they could in principle be tuned. For example, the application of hydrostatic pressure can preferentially shift the optic phonons to higher frequencies relative to the acoustic phonons. This would weaken AAO processes and possibly give unusually rapid increase of κ in SiC and BP. Similar behavior is found for the optic phonons. The compounds, AlAs, AlP, BP, GaAs, GaSb, GaP, InAs and InSb have regions of weaker AOO scattering rates due to the AOO#2 selection rule than does AlSb. Thus, even larger contributions from optic phonons would be predicted from the lowest order-theory for these compounds if not for their comparably small A-O gaps, which give strong AAO scattering where the AOO optic phonon scattering is small. As a result, optic phonons in these compounds contribute negligibly to κ even in the lowest-order theory.
Complementary features activating AAA#1 and AAA#2 selection rules: Figure 6 shows that those compounds that more strongly activate AAA#1 or AAA#2 selection rules lie on opposite ends of the velocity ratio plot. Thus, compounds that have acoustic branches that are bunched together away from the center of the BZ (BAs, BSb, BP, BN, SiC, AlN, GaN) and so activate the AAA#1 selection rule, generally have smaller differences between near-zone center LA and TA velocities. In contrast, compounds with large differences in their near-zone center LA and TA velocities (e.g. InSb, InAs, InP, GaSb), and so activate the AAA#2 selection rule, do not show the acoustic branch bunching behavior or do not have it exposed by the Type 1 selection rules.
Importance of new Type 2 selection rules: Prior studies have already identified the importance of the AAA#1 Type 2 selection rule in explaining the anomalously weak three-phonon phase space in BAs and BSb. Yet, this selection rule cannot explain the unusually large optic phonon contributions to κ (3) Pure in BAs and AlSb, nor can it explain the large acoustic phonon contributions to κ (3) Pure in InP, InAs, InSb and Gasb. The two new Type 2 selection rules, AAA#2 and AOO#2, are required for full understanding of the behavior in the studied ZB materials.
Restricted spectral contributions to thermal conductivity: In all of the studied compounds, we find that the dominant contributions to κ occur below the highest TA phonon frequency, ν max TA . Above ν max TA , LA phonons exhibit strong decay via AAA processes as well as participating in AAO processes resulting in small contributions to κ in this frequency region. The ν max TA are marked in the spectral κ-contribution figures for each material (Figs. 5(d)-(f), 7 (d)-(f) and S16-S18) to illustrate this behavior.
Distinction between phase space and matrix element selection rules: The selection rules discussed in the present work are connected to the phase space for threephonon scattering. They depend on specific features being present in phonon dispersions, as described in Section II. It is important to distinguish these from a different kind of selection rule, which can cause the matrix elements for phonon-phonon scattering to vanish. Such selection rules have been identified in systems of reduced dimensionality and are connected to the underlying symmetry in some crystals [58][59][60][61]. Their incorporation can improve efficiency in schemes for computing phonon-phonon processes, as recently discussed [61].
Effect of impurities on the measured κ: The large differences seen between the κ values calculated only using the lowest-order theory and the measured data for many of the studied materials cannot be explained by possible impurities in the measured samples. The presence of impurities gives a weaker temperature dependence to κ compared to the predictions from the lowest-order theory, since phonon-impurity scattering is temperature independent, as described previously [17]. In contrast, the measured κ (T ) shows a stronger temperature dependence than seen in many of the lowest-order calculations. The calculations including both three-phonon and four-phonon scattering processes also gives stronger T-dependence [13,17] and are both qualitatively and quantitatively in agreement with the experimentally measured temperature-dependent κ.
Impact of selection rules in other materials: The selection rules on three-phonon scattering are completely general and are not restricted to ZB compounds. Creation of frequency windows in Fig. 2 benefits from materials that have low ionicity (e.g., several ZB compounds in this study), which gives relatively small optic phonon bandwidths, and large heavy-to-light mass ratios. For materials with more than two atoms in the unit cell, the additional scattering channels resulting from the increased number of optic phonon branches will make the needed frequency windows harder to achieve for acoustic phonons, so the selection rules are less likely to affect the acoustic phonons in such compounds. But, large frequency gaps between bands of optic phonons can be found in many compounds [62] so erroneously large anharmonic phonon decay rates and contributions to κ from some optic phonons could occur from activation of selection rules such as AOO#2 and OOO. For such compounds, the selection rule picture can provide useful guidance for understanding measurements giving lower κ and stronger temperature dependence of anharmonic phonon decay than predicted from calculations including only three-phonon scattering. In crystals with one atom in the unit cell, only AAA processes occur. We are not aware of any such materials that are semiconductors or insulators, and in metals, κ is typically dominated by the electronic contribution. For the lattice contribution, inclusion of electron-phonon interactions could mask any influence of the selection rules on the lowest-order phonon-phonon interactions.

VI. SUMMARY
Ab initio theories of the anharmonic properties of materials continue to rely on descriptions that include only the lowest-order phonon-phonon interactions. The vastly different results obtained here for phonon lifetimes and thermal conductivity with and without four-phonon interactions for the many considered materials contradict prevailing understanding of anharmonic phonon decay and thermal conduction in weakly anharmonic crystals, where four-phonon interactions have frequently been assumed to be unimportant.
We have identified new selection rules on threephonon scattering dictated entirely by features in the phonon dispersions and have shown that the full set of selection rules listed in Section II and the trends found in four-phonon scattering strengths provide a useful framework to understand the anharmonic phonon decay and thermal conductivity of weakly anharmonic insulating crystals. In particular, the new selection rules are critical in explaining the failure of the lowest-order theory in describing thermal transport in many zinc blende compounds.
In compounds that do not contain boron, carbon or nitrogen, four-phonon scattering suppresses the erroneously large phonon lifetimes and thermal conductivities predicted by the lowest-order theory. Conversely, in compounds that do contain boron, carbon, and nitrogen atoms, large phonon lifetimes and thermal conductivities predicted in the lowest-order theory are less affected, because of the exceptionally weak fourphonon interactions. These findings help explain the high thermal conductivities found in the technologically important compounds such as BAs, BP, and SiC, and provide critical guidance in the search for new materials with high thermal conductivity.
Inclusion of four-phonon interactions in ab initio anharmonic decay rates and transport calculations is computationally challenging, and such calculations are currently performed by only a small number of groups [12][13][14][15][16][17][18][19][20][21]. Recent work has shown that for strongly anharmonic materials, it is essential to include four-phonon processes [14][15][16]. The present study shows that, even in weakly anharmonic materials that do not contain B, C or N atoms, much stronger four-phonon scattering occurs, and points to the possible need to include four-phonon scattering in future thermal transport calculations even in such strongly bonded compounds. While this suggestion has been made before [13], the large number of compounds studied in the present work and comparably larger four-phonon scattering strengths found in the non-BCN compounds adds strong quantitative and statistical support to this proposition. Further calculations will be required to fully establish this point.

VII. METHODS
For all materials, temperature dependent lattice parameters are obtained by minimizing the Helmholtz free energy F = Φ 0 + F H + F A , where Φ 0 is the energy of the lattice atoms in their equilibrium positions, and F H is the harmonic part and F A is the anharmonic part, which includes both third and fourth order terms [14]. To obtain phonon modes, the second-order interatomic force constants (IFCs) are calculated using density functional perturbation theory (DFPT) as implemented in the QUANTUM ESPRESSO package [63]. All calculations in this work were performed using norm-conserving pseudopotentials with the local density approximation exchange correlation functional. Converged parameters for the density functional theory (DFT) calculations are provided in the SM section S14.
The third and fourth-order anharmonic IFCs required to obtain three-phonon and four-phonon scattering rates are calculated using a thermal stochastic snapshot technique. The method has been described in detail in Ref. [14].
Temperature dependent anharmonic renormalization of phonon modes [14] was included but did not significantly affect phonon modes, scattering rates and thermal conductivities for the considered materials and temperature ranges. To calculate the thermal conductivity of a material, the Boltzmann equation for phonon transport including three-phonon, four-phonon and phonon-isotope scattering was solved using an iterative approach for the non-equilibrium phonon distribution function created from an assumed small temperature gradient across the considered materials. For all materials and all temperatures, the phonon Boltzmann equation was solved on a 51X51X51 q-grid. Four-phonon scattering rates were obtained on 17X17X17 q-grids and interpolated to the finer 51X51X51 q-grids (see SM section S13 for details). Full calculations (without interpolation) were also performed on 21X21X21 q-grids to check the convergence of the four-phonon scattering rates (see SM section S12). The phonon Boltzmann equation and the expressions for the scattering rates are given in the Appendix.

VIII. ACKNOWLEDGMENTS
where, 1/τ The three-phonon scattering probabilities are: The four-phonon scattering probabilities are: The three-phonon and four phonon matrix elements are Φ λλ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 , π) (A.11) 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 , ρ) (A.12) respectively, where w α (qs, µ) is the α th component of the eigenvector w (λ, µ) = w (qs, µ) for a phonon with wavevector q and polarization s and for the basis atom µ, and N 0 is the number of atoms in the supercell (or equivalently, the number of q-points in the commensurate first BZ). The phonon-isotope scattering probabilities are: where g 2 (b) = 1/M 2 b a f ab M ab −M b 2 is the mass variance parameter with f ab and M ab being the concentration and mass of the a th isotope of the b th atom, respectively, andM b is the average mass of the b th atom. We find that the momentum-relaxing Umklapp processes are comparable to or stronger than the momentumconserving normal processes for the four-phonon interactions for all seventeen materials and at all temperatures considered in this work (see SM section S11). As a result, treating the four-phonon scattering within a relaxation time approximation (RTA), i.e., removing the last bracketed term in Eq. A.3, gives negligible difference compared to the full solution of Eq. A.3. We have confirmed that (a) solving the full phonon Boltzmann equation, which explicitly differentiates normal and Umklapp three-phonon and four-phonon interactions and (b) solving the partial phonon Boltzmann equation with fourphonon scattering under the RTA by removing the last bracketed term in Eq. A.3 results in negligible difference in κ (3+4) at 300 K for all seventeen materials in this study. Therefore, the calculations at all other temperatures presented in this work are obtained by solving the phonon Boltzmann equation with the four-phonon scattering terms treated under the RTA, to reduce the computational cost of these temperature-dependent calculations for the seventeen materials. After solving the phonon Boltzmann equation for F λ , the phonon thermal conductivity is calculated as: where C λ = (1/V ) k B ∂n 0 λ /∂T is the mode specific heat.