{'fraction_non_alphanumeric': 0.03531026173103956, 'fraction_numerical': 0.014291451904411272, 'mean_word_length': 4.727237876174047, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Success and breakdown of the T-matrix approximation for phonon-disorder scattering

We examine the validity of the widely used T-matrix approximation for treating phonon-disorder scattering by implementing an unfolding algorithm that allows simulation of disorder up to tens of millions of atoms. The T-matrix approximation breaks down for low-energy flexure phonons that play an important role in thermal transport in two-dimensional materials. Furthermore, insights are developed into the success of the T-matrix approximation in describing maximally mass disordered systems. To achieve this, the phonon unfolding formalism is generalized to describe mass disorder and strongly nonperturbative features of the spectrum are connected to the Boltzmann quasiparticle picture.


I. INTRODUCTION
Heat management has become an important technological and environmental issue, now critically considered for further electronics miniaturization and for realistic solutions of global energy challenges. As two thirds of the energy consumed worldwide is dissipated as waste heat, 1 the development of more efficient heat management strategies with regards to thermal insulation, transport, and conversion will drive transformative industrial, economic, and environmental impacts. [2][3][4][5][6] However, scientific questions remain regarding thermal transport in materials with defects, which can limit thermal functionalities, but also provide an important tuning mechanism for engineering thermal transport, particularly for thermal insulation and thermoelectric applications. [7][8][9][10][11][12][13][14] Widely adopted first principles transport methods 15,16 have been coupled with combinations of secondorder perturbation theory (i.e. Fermi's golden rule [17][18][19] ) and Green's function methods (i.e. the Tmatrix approximation 20,21 ) to evaluate the influence of isotopes, [22][23][24] vacancies, 18,25,26 atomic substitutions, 17,27 and mass and force constant disorder 28-36 on thermal transport in nanoscale, bulk, and alloy systems. The T-matrix approximation, in particular, is computationally accessible, even from first principles density functional theory methods, includes resonant scattering effects, and is not limited by linewidth resolution, which can be crucial for calculating thermal transport contributions from low energy, long mean free path phonons. However, success of perturbative methods and the Tmatrix approximation in determining thermal conductivities of maximally-disordered alloys 30,33,36 is surprising due to their strong mass variations and defect concentrations well beyond the dilute limit. Furthermore, phonon unfolding methods in alloys with large mass differences have demonstrated strongly altered phonon spectra, including the appearance of heavy and light branches where only acoustic branches exist in the virtual crystal approximation (VCA). 37 Such strong distortions of the vibra- tional structure are difficult to reconcile with approaches based in the phonon quasiparticle picture. An important class of crystals for which the validity of the independent defect scattering approximation may be questionable is two-dimensional (2D) and van der Waals (vdW) layered materials. These systems have received significant attention over the past decade due to their interesting electronic and vibrational properties and potential for device applications. [38][39][40] Monolayer materials feature peculiar out-of-plane modes characterized by quadratic dispersion near the Brillouin zone center Γ, reflecting the ease with which the atomic sheet can vibrate out of the plane. This can be illustrated by a simple two-dimensional spring-mass model sketched in Fig. 1, featuring effective out-of-plane springs between the nearest and third-nearest neighbors. The latter negative spring constants are required to satisfy the rotational sum rules (which are equivalent to the absence of residual forces and stress) that give rise to quadratic dispersion. 41 Without these, the sum rules are broken and the flexure branch has linear dispersion as in a strained system. These quadratic out-of-plane modes can carry a large proportion of the heat in 2D and vdW systems, [42][43][44] and both T-matrix and phonon unfolding studies of defect scattering in graphene suggest that they are sensitive to disorder, particularly vacancies. 26,45 However, evaluation of the validity of the T-matrix approximation for phononvacancy scattering of the quadratic modes has not been carried out. Such an evaluation is important as prominent monolayer compounds such as MoS 2 can have large vacancy concentrations (beyond the dilute limit); a few percents are not atypical depending on synthesis methods, and up to 20% is possible if engineered. [46][47][48] More generally speaking, developing insights into the range of validity (strength of perturbations and defect concentrations) for prominent T-matrix methods is critically important for building confidence in defect-limited first principles transport predictions.
In this article, the validity of the T-matrix approximation is evaluated in generic two and three-dimensional models with vacancies and strong mass disorder. To this end, the phonon spectrum of large disordered supercells of tens of millions of atoms is unfolded using the non-perturbative Chebyshev polynomials Green's function method. 49 We demonstrate that low-energy out-ofplane vibrations in 2D monolayers are especially sensitive to multiple-impurity scattering effects. We find that the usual unfolding procedure to calculate the spectral function must be generalized to treat mass disorder, and that strong mass variance gives rise to disorder-induced dispersive branches. We connect these results to the Boltzmann transport picture to understand the success of perturbative and T-matrix methods in predicting the thermal conductivity of bulk alloys despite strong distortions of their phonon spectra.
This manuscript is organized as follows: the T-matrix approximation and phonon unfolding method are briefly explained in section II. Section III focuses on phonondefect scattering in atomic monolayers featuring vacancies. In section IV is discussed the case of mass disordered alloys. We summarize our study in section V.

II. METHODS
We focus here on three approaches for calculating phonon-defect scattering rates, Fermi's golden rule (FGR), the T-matrix approximation (TMA), and the supercell phonon unfolding method (SPU). For other technical approaches we refer to Ref. 35,50-52. The T-matrix method approximates the self-energy of the Green's function by a subset of diagrams corresponding to successive phonon scatterings with a single impurity. 53 As such, it is most appropriate in the limit of low defect concentrations and is expected to break down when multiple-impurity scatterings become important. For TMA calculations presented here, the scattering rate of the Bloch mode with wavevector q is given by 20,21,54 where r runs through the atoms, x is the impurity concentration, ω 0 q is the frequency of mode q in the clean system, and T is the T-matrix: T (ω) = P (ω)[1 − G 0 (ω)P (ω)] −1 . P is the perturbation of a single impurity: P rr (ω) = ∆φ rr /m 0 − ω 2 δ rr ∆m r /m 0 , with ∆φ rr the change in the interatomic force constant (IFC) between r and r , and ∆m r = m r − m 0 is the mass variation of r. G 0 (ω) = [(ω + iη) 2 − D] −1 is the Green's function of the clean system, with D rr = φ rr /m 0 the dynamical matrix and η an infinitesimally small positive imaginary part. The direction indices have been dropped as different directions are decoupled in our models. Due to the local nature of the perturbation, the dimensions of P and T are relatively small and matrix inversions and multiplications are straightforward. To obtain the scattering rate from FGR, the T-matrix in Eq. (1) is replaced with the second-order perturbation correction: P (ω) + P (ω)G 0 (ω)P (ω).
In recent years, headway has been made in the study of vibrations in strongly disordered materials via the supercell phonon unfolding method. 25,37,45,50,[55][56][57][58][59][60][61] In this non-perturbative numerical technique, the phonon spectral function is calculated for a disordered supercell according to: where r runs through the N atoms in the supercell, Q is the wavevector in the supercell first Brillouin zone that is equivalent to q in the supercell reciprocal lattice, ω QJ is the supercell eigenmode frequency of band J at wavevector Q, and V r,QJ is the corresponding eigenvector. For the case of force constant disorder the phonon lifetimes can be extracted from the full-width at half-maximum of the spectral peaks given by Eq. (2) (see appendix B). For the case of mass disorder, this correspondence breaks down as will be discussed below. With increasing supercell size, the unfolded spectral function converges to the exact spectral function of the fully disordered system. In practice, the spectral function is often computed through an exact diagonalization of the supercell dynamical matrix. However, the unfavorable N 3 scaling of this approach constrains the supercell size to around 10 4 atoms or less. The resulting finite-size effects in turn limit the energy resolution attainable. This is especially problematic for phonons because the linewidths of the acoustic modes become vanishingly small as q → Γ and these often contribute significantly to the thermal conductivity. To circumvent this limitation, the spectral function is evaluated here via the Chebyshev polynomial Green's function method. 49,62,63 In this approach, the Green's function is expanded on the Chebyshev polynomial basis (see appendix D), yielding a linear scaling with respect to N . In this way, we are able to reach very large supercell sizes of tens of millions of atoms thus enabling investigation of the low-lying acoustic modes. For such system sizes, averaging of multiple disorder configurations and going beyond Q = 0 is not necessary.

III. VACANCIES IN 2D SYSTEMS
We consider vacancies in the 2D spring-mass model ( Fig. 1) that exhibits a generic quadratic dispersion. It features atoms of mass m on a square lattice, nearestneighbor springs of stiffness k and third-nearest neighbor springs of stiffness −k/4. All angular frequency and inverse lifetime values are presented in units of k/m. A concentration x of vacancies is introduced by removing the springs between the missing atom and its neighbors. Springs around the vacancies are modified in order to satisfy the rotational sum rules for the IFCs. (see appendix A) Fig. 2 gives the ratio τ SPU q /τ TMA q , where τ TMA q and τ SPU q are the disorder-induced phonon lifetimes calculated by the TMA and SPU, respectively, for vacancy concentrations x = 1%, 2%, 5%, and 10%. The inset shows the inverse lifetimes divided by x. Both methods give a quadratic power law at low frequency, as expected from FGR. 64 The T-matrix predicts the lifetime at angular frequencies above 0.5 reasonably well, apart from the discrepancies at 2 due to a Van Hove singularity and at ≈ 2.8 due to the vanishing of the density of states. However, the T-matrix fails at lower frequencies, especially in the 0.15 − 0.4 range. This region, marked by a green shading in Fig. 2, corresponds to the bottom 15% of the out-of-plane spectrum, which typically represents a significant part of the total thermal conductivity in 2D materials (for instance, see appendix E for the case of MoS 2 ). More precisely, in this highly-relevant region, SPU lifetimes can be as low as 55% of the TMA values for x = 1%, and 40% for x = 2%. The ratio can reach 25% for x = 5% and 15% for x = 10%. This suggests that the T-matrix approximation breaks down in 2D materials with high vacancy concentrations. Vacancy concentrations larger than 1% are not uncommon in transition metal dichalcogenides, for instance, which may result in a significant overprediction of defect-limited thermal transport from the TMA. In such cases, accurately predicting the thermal conductivity requires consideration of multiple-impurity scatterings, relying for instance on the SPU.

IV. MASS DISORDER
Mass disorder from isotopes or atomic substitutions presents an important phonon thermal resistance in materials, particularly in alloys such as (Si,Ge) or Mg 2 (Si,Sn). Mass perturbations are diagonal in real space and proportional to ω 2 , which implies that perturbative methods and the TMA should become exact at low energy. 21,65 The questions remain: what is the range of validity for these approaches, and does the low frequency behavior explain the success of these techniques in computing thermal conductivities of strongly disordered alloys?
Regarding the SPU method, the phonon unfolding formula (2) used throughout the literature 37,55-58,60,61 does not probe the spectrum (and thus the lifetime) of a Bloch vibrational mode in the presence of mass disorder. Because the eigenmodes are defined using mass-rescaled atomic displacements, Eq. (2) corresponds to a mode in which the true atomic displacements are proportional to e iq·r √ N 1 √ mr , (see appendix B) with m r being disordered. The atomic displacements corresponding to actual Bloch modes with ordered masses should be used instead. In the present context, it is helpful to look to the Green-Kubo formula, which is exact beyond the quasiparticle picture in the linear regime without phonon-phonon interactions. Neglecting vertex corrections, the thermal conductivity κ can be expressed (see appendix C) as q and v 0 q the Bloch mode frequencies and velocities along direction x and is the generalized phonon unfolding formula for the spectral function of mass disordered systems that captures the lifetime of the Bloch mode q (see appendix B). m 0 is the mass in the ordered reference system. Here we assess the effects of mass disorder in bulk materials by calculating the spectral function of a threedimensional (3D) spring-mass model on a cubic lattice featuring nearest-neighbor springs k and clean atomic masses m, with pure mass defects (no IFC perturbations). The mass differences in thermally-relevant alloys can reach a factor four, as in Mg 2 (Si,Sn) or (In,Al)N alloys. Therefore, we introduce impurities with mass m = 4m as a reasonable upper limit to the strength of mass fluctuations and focus on maximally-disordered alloys (half the atoms are substituted). For the reference system we use the VCA: m 0 = m+m 2 . Fig. 3 (a) gives the spectral functionÃ(q, ω) obtained from the generalized phonon unfolding formula, Eq. (3). An artificial Lorentizan broadening of 0.003 has been kept for visualization purposes. Strikingly, the single branch of the clean system separates into an acoustic-like branch at low-energy and an optic-like branch at high energy. Despite these strong disruptions of the vibrational structure, the VCA predicts the phonon dispersion correctly for the acoustic-like branch below ω = 0.5 due to the decrease of the mass perturbation as ω approaches zero. Fig. 3 (b) gives the calculated projected density of states (pDOS) for the light and heavy atoms via exact diagonalization of the dynamical matrix of 100 supercells containing 400 atoms on average. The acoustic-like branch is governed primarily by the heavy atom vibrations, while the optic-like branch arises from vibrations of the light atoms against the matrix of heavy atoms (this is qualitatively consistent with calculated vibrational structures of certain binary alloys 37 ). Instead of a simple broadening of the VCA branch, two distinct branches are present with equally strong quasiparticle peaks. As these starkly non-perturbative features of the vibrational spectrum demonstrate, the phonon quasiparticle picture that the Boltzmann transport formalism is built upon breaks down, calling into question the validity of previous transport calculations in strongly disordered alloys.  We assume for simplicity that the spectral function near the resonances can be reasonably approximated by separate Lorentzians:Ã(q, ω) = l F ql 1 π 1/2τ ql (ω−ω ql ) 2 +(1/2τ ql ) 2 , where l = 1, 2 designates the heavy and light branch, respectively, and F ql are quasiparticle weights as given by the spectral function. Using the approximations necessary for derivation of the kinetic equation from the Green-Kubo formula, 66 the thermal conductivity can be further simplified:

DOS
This is similar to the kinetic particle description of phonon conductivity for a system with one acoustic branch and one optic branch. Indeed, this equation connects the SPU formalism to the Boltzmann transport theory in the case where strong disorder generates multiple branches. It allows the different peaks to be interpreted as quasiparticle excitations and their broadenings 1/τ ql as inverse lifetimes. Remarkably, the heat flux is given by the products of the VCA frequencies ω 0 q and velocities v 0 q renormalized with the quasiparticle weights F ql , not by the quasiparticle frequencies ω ql and velocities ∂ω ql /∂q. The broadening of the acoustic-like branch at low energy is well below 0.01, much smaller than in the rest of the spectrum, including the whole optic-like branch whose broadening is of order 1. We extracted the inverse phonon lifetimes from the width at half-maximum of the acoustic-like spectral peaks in the high-symmetry segments Γ − X, Γ − M and Γ − R of the Brillouin zone. The ratios τ SPU q /τ TMA q and τ SPU q /τ FGR q are shown in Fig. 4, where τ FGR q is the phonon lifetime calculated by FGR 19 (the perturbation is defined with respect to the VCA). The inverse lifetimes themselves are shown in the inset. At high frequencies, the phonon lifetimes are significantly underestimated by the TMA and FGR compared with the SPU values. However, all techniques give an ω 4 power law at low energy, as expected for point defect scattering. There is good agreement (within 15%) between all methods below ω = 0.4, which roughly corresponds to the bottom 20% of the VCA frequency range, as indicated by the shaded region in Fig. 4. These low frequency modes give a significant contribution to heat transport in disordered alloys. In Si 0.5 Ge 0.5 for example, we find from first-principles calculations that 90% of the thermal conductivity is derived from the bottom 20% of the frequency range (see appendix E). Although definitive pronouncements on any specific material require more realistic calculations, this strongly suggests that the ω 2 dependence of the mass perturbation is the main explanation for the unexpected success of T-matrix and perturbative approaches for calculating alloy thermal conductivity. As we examined the extreme case of large mass ratio and maximal disorder, it is expected that agreement of the different techniques will only improve with weaker disorder.

V. SUMMARY
The validity of the T-matrix approximation has been evaluated by comparing disorder-induced lifetimes in generic mass-spring models using the phonon unfolding Chebyshev polynomials Green's function method. This non-perturbative technique allows the study of disordered supercells containing tens of millions of atoms and fully includes multiple-impurity scatterings. In order to include mass variance, we have generalized the unfolding formalism for the spectral function and connected the strongly distorted vibrational structure to the phonon quasiparticle picture using the Green-Kubo formula as a guide. Second-order perturbation theory and the Tmatrix approximation are valid at low frequencies in systems dominated by mass disorder, explaining the success of thermal conductivity predictions in various alloys. However, these methods break down in a similarly relevant frequency range for out-of-plane vibrations in two-dimensional systems containing 1% of vacancies or more. In these cases, the phonon unfolding methodology presented in this work captures the relevant multipleimpurity scattering effects.
Acknowledgements. This research was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division. We used resources of the Compute and Data Environment for Science (CADES) at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. In two-dimensional materials, in-plane (xy) and outof-plane (z) vibrational modes have different characteristic behaviors. The former have a linear dispersion around Γ, while the latter exhibit a quadratic dispersion. This quadraticity reflects the greater ease with which a two-dimensional atomic layer can vibrate out of the plane than in the plane, and disappears if strain is applied to the material 67 . In practice, the interatomic force constants (IFC) must satisfy a set of rotational sum rules called the Born-Huang conditions and the Huang invariance 68 , that are equivalent to the absence of residual forces and strains in the material 41 . To preserve the quadraticity of the phonon dispersion, the out-of-plane IFCs of our two-dimensional model must therefore obey the following constraints: where α, α = x, y, φ rr is the out-of-plane IFC between atoms r and r , and r α is component α of the position of atom r. We consider a spring-mass model on a square lattice, with atoms of mass m and springs of stiffness k between the first nearest-neighbors. To satisfy condition (A2), it is necessary to introduce a third nearest neighbor interaction with a spring stiffness −k/4. Note also that the IFCs satisfy by construction the translational sum rule r φ rr = 0 ensuring the existence of zero-energy modes at Γ. When a vacancy is introduced in the 2D model, equations (A1) and (A2) must still be satisfied. Therefore, we add two modifications of the IFCs: the third-neighbor springs across the vacancies are cut and the spring constants between the four vacancy nearest-neighbors and their first neighbor in the direction opposite the vacancy are reduced by k/2. The rotational sum-rules being satisfied, the dispersion relation calculated through the Chebyshev polynomials Green's function method (CPGF) remains quadratic even at high vacancy concentration.

Appendix B: Phonon unfolding and the spectral function
In a clean, ordered system containing N sc N unit cells of a crystal with Born-von Karman periodic boundary conditions, each unit cell contains several atoms i of mass m 0 i , and r denotes the sites of the Bravais lattice. The displacement of atom i in cell r along direction α from its equilibrium position is given by the quantum operator u α ir , and its momentum is given byp α ir . (φ 0 ) αα ir,i r denotes the IFC between atom i in cell r along α and atom i in cell r along α . The Hamiltonian is with the reduced coordinatesÛ 0 α ir = m 0 iû α ir andP 0 α ir = p α ir / m 0 i and the clean dynamical matrix (D 0 ) αα ir,i r = (φ 0 ) αα ir,i r / m 0 i m 0 i . Going to the Fourier domain and considering a wavevector q in the first Brillouin zone (BZ) of the lattice yields a q-dependent dynamical matrix, which can be diagonalized: with j the branch index and e α i,qj the polarisation vectors, i.e. the normalized eigenvectors of D 0 q with (ω 0 qj ) 2 the corresponding eigenvalues. The Hamiltonian is diagonal in terms of the Bloch displacement and momentum eigenoperators: with the momentum eigenoperators defined similarly. E α ir,qj = e iqr e α i,qj / √ N sc N is a normalized eigenvector of the real-space dynamical matrix (D 0 ) αα ir,i r . Then, let us divide the system in N sc supercells (SC) with positions R containing N units each. We introduce mass and IFC disorder (m 0 → m, φ 0 → φ, D 0 → D) identically in every SC, so that the system is unchanged by translations from one supercell to another. The Hamiltonian is now with the new reduced coordinatesÛ α IR = √ m Iû α IR and P α IR =p α IR / √ m I , and I running over the atoms in a supercell. The Fourier transformed dynamical matrix can again be defined as a function of Q, a wavevector in the reduced BZ of the new Bravais lattice of the supercells: with J the new branch index and V α I,QJ the normalized eigenvectors of D Q with ω 2 QJ the corresponding eigenvalues. The displacement eigenoperators arê and the Bloch operators can be expressed in terms of the disordered eigenoperators: (B7) where r runs over the unit cells contained in the central supercell only, and Q(q) is the wavevector in the reduced BZ of the supercell lattice that is equivalent to q.
The phonon retarded Green's function of the disordered system can be defined on the eigenbasis from the correlation function of the displacement operators 69 : with Θ the Heaviside function and |0 the phonon ground state of the disordered system. Its Fourier transform is G QJ,Q J (ω) = dtG QJ,Q J (t)e i(ω+iη)t : η being a positive infinitesimal value. The decay of the Bloch mode qj in the disordered system is described by the auto-correlation function of the Bloch displacements: From here, the generalized unfolding formula for the spectral function is which yields for ω > 0: Note that if the factor m 0 i mir is omitted in equation (B12) in the presence of mass disorder, then what is computed is the spectrum of the disordered vibration modê The spectral function is then defined as yielding for ω > 0: (B15) For the simple models considered in our study, equations (B15) and (B12) reduce to equations (2) and (3) of the manuscript, respectively.
Appendix C: The Green-Kubo formula and the spectral function In the absence of many-body interactions between phonons, the lattice thermal conductivity is given by the Green-Kubo formula. It can be expressed as (see eq. with S 0 αα ir,i r the heat current matrix in the clean system. The Green's function matrix is diagonal in the disordered eigenbasis: and the clean heat flux is diagonal in the Bloch basis (we neglect interband terms 70 ): where ω 0 qj and v 0 qj are the frequency and the group velocity of the eingenmode qj in the clean system. Thus, taking the trace on the Bloch vector basis E qj : where M IR,I R = δ RR δ II m I is the matrix of the disordered masses and M 0 ir,i r = δ rr δ ii m 0 i is the matrix of the clean masses. If we neglect the cross-terms q = q (i.e. we neglect the interferences between the two Green's functions, that is to say the vertex corrections), then the thermal conductivity along the normalized vector u becomes: where v 0 q = v 0 q · u. Noting that Im 2ω (ω+iη) 2 −ω 2 QJ = −πδ(ω − ω QJ ), we find: Appendix D: The Chebyshev polynomial Green's function method The Chebyshev polynomials Green's function (CPGF) approach has been reviewed in Refs. 49,63 for the case of electrons, and has been adapted in Ref. 45 to investigate phonon lifetimes in defected graphene. Here a brief overview of the method will be given. The phonon Green's function of a large disordered supercell is expanded on the Chebyshev polynomial basis: where the bar indicates that the spectrum has been rescaled to [−1, 1], the g n (z) are known complex functions:  Figure 5: Calculated room temperature normalized thermal conductivity accumulation as a function of frequency for Si (black curve), Ge (purple curve), and a Si0.5Ge0.5 alloy (red curve). The thermal conductivities are scaled to the corresponding calculated room temperature bulk values (including phonon-isotope scattering): 144.9 W/m K for Si, 60.7 W/m K for Ge, and 9.7 W/m K for Si0.5Ge0.5. The frequencies are scaled by the highest calculated frequency in each system: 15.6 THz for Si, 9.0 THz for Ge, and 11.2 THz for Si0.5Ge0.5 within the virtual crystal approximation. and the T n (D) are the Chebyshev polynomials evaluated at the dynamical matrix, that follow the recursion relation T n+1 (D) = 2DT n (D) − T n−1 (D).
The spectral function (B12) for the Bloch mode qj can be expressed as so in practice, the quantities to be calculated are the socalled moments µ n,qj : which are computed using the recursion relation between the T n (D). The number of moments necessary for the sum (D1) to converge is roughly equal to 1/2ωη. Becauseη is an artifical broadening and should be smaller than the disorderinduced spectral linewidth, probing modes closer and closer to Γ requires more and more polynomials to be included.  for Si, Ge, and a Si 0.5 Ge 0.5 . Here, ω qj is the frequency of phonon mode with wavevector q and polarization j, C qj is the mode specific heat, v qj is the phonon velocity in the direction of an applied heat flux, τ qj is the phonon transport time in this direction, and Θ is the Heaviside step function. The conductivities and frequencies are normalized as described in the figure caption. For Ge and Si, the harmonic and anharmonic interatomic force constants that determine the phonon properties in Eq. (E1) are determined from density functional theory as described in Ref. 22. For these calculations, full solution of the Boltzmann transport equation 22,71 determines the phonon lifetimes as limited by three-phonon and phonon-isotope interactions 19 . For Si 0.5 Ge 0.5 all properties (harmonic and anharmonic force constants, masses, lattice constants, etc.) of separate Si and Ge systems were averaged as in the virtual crystal approximation as described in Ref. 33. Additional phonon scattering from Si/Ge mass disorder was included via perturbation theory exactly as done for phonon-isotope interactions 19,33 . Fig. 6 gives the calculated contribution from out-ofplane phonons to the thermal conductivity accumulation as a function of phonon frequency for monolayer MoS 2 and monolayer MoS 1.96 (2% S vacancies). The contribution from out-of-plane phonons is calculated by multiplying the contribution from each mode by the component of the polarization vector along the out-of-plane direction. Details of the first principles thermal conductivity calculations for MoS 2 are given in Ref. 25