\section{Introduction}

Vibrational thermal properties of condensed matter systems play critical roles in determining material stability and synthesizability, properties and functionalities, and resulting pathways towards building useful multicomponent, multiscale devices for a variety of technologies \cite{champierThermoelectricGeneratorsReview2017,bellCoolingHeatingGenerating2008,mooreEmergingChallengesMaterials2014,shinAdvancedMaterialsHighTemperature2019}. Often material synthesis pathways are designed to reduce crystal imperfections for improved properties such as low electrical and thermal resistances. However, sometimes defects and disorder are intentionally created to improve structural and elastic behaviors or increase thermal resistance \cite{horiTuningPhononTransport2018,beekmanInorganicCrystalsGlass2017,hanusThermalTransportDefective2021}. Alloys represent such a class of materials which have played critically important roles in technological advancements over the course of human history, e.g., bronze tools and weapons, steel construction materials, and high temperature aluminum alloys in combustion engines. 

In terms of thermal management applications, a variety of novel alloys have been developed and employed, particularly towards reducing vibrational contributions to thermal transport for more efficient thermoelectrics – example alloys include Si$_{1-x}$Ge$_x$ \cite{basuHighTemperatureSi2021}, PbTe$_{1-x}$Se$_x$ \cite{lalondeLeadTellurideAlloy2011}, half-Heusler alloys \cite{huangRecentProgressHalfHeusler2016}, and Mg$_2$Si$_{1-x}$Sn$_x$ \cite{pandelReviewMg2Si2021,farahiHighEfficiencyMg22019}. To build insights into the vibrational behaviors of these materials, density functional theory (DFT) based Peierls-Boltzmann transport models have been developed and deployed to examine intrinsic and extrinsic phonon scatterings and resulting thermal conductivities ($\kappa$) as a function of alloy concentration ($x$) \cite{gargRoleDisorderAnharmonicity2011,liThermalConductivityBulk2012,tianPhononConductionPbSe2012}. These models have been surprisingly successful despite the array of \textit{ad hoc} approximations employed: virtual crystal approximation (VCA), perturbation theory for mass disorder (Fermi’s golden rule), and lack of force disorder. In particular, the VCA ignores short range variations in the alloy systems and perturbation theory is not well-justified for systems with strong disorder, as is the case for the alloys described above. A case may be made for ignoring force disorder as the alloying elements are often isoelectronic; however, more recent DFT Green’s function calculations of phonon-disorder scattering and $\kappa$ in In$_{1-x}$Ga$_x$As have demonstrated that even this assumption is not fully valid \cite{arrigoniFirstprinciplesQuantitativePrediction2018b}. 
     
Building on previous work featuring mass disorder in spring-mass models \cite{thebaudSuccessBreakdownTmatrix2020}, here we critically examine the widely used VCA and perturbation approximations for describing phonon-mass-disorder scattering in alloyed semiconductors from DFT methods with specific application to maximally mass disordered Mg$_2$Si$_{0.5}$Sn$_{0.5}$. We confront this standard methodology with a more rigorous non-perturbative approach, the Chebyshev polynomials Green’s function method (CPGF) \cite{thebaudSuccessBreakdownTmatrix2020,bouzerarDrasticEffectsVacancies2020,muUnfoldingComplexityPhonon2020,allenRecoveringHiddenBloch2013}, thus isolating the effect of the perturbative approximation on the phonon lifetimes and on the thermal conductivity. We find that the general VCA phonon quasiparticle picture breaks down above a certain frequency but succeeds in describing low-frequency long-wavelength acoustic phonons that carry much of the heat in this system.  Moreover, the presence of an ordered Mg sublattice reduces the sensitivity to disorder of dispersive optic phonons, enabling them to play a significant role in heat transport. In light of these results, we discuss the conditions under which we expect the thermal conductivity to be well-predicted by perturbation theory.

\section{Mass disorder models}

In this section we describe calculations of the phonon spectrum and phonon transport in mass-disordered alloys from (i) the virtual crystal approximation and Fermi’s golden rule (quantum perturbation theory) and (ii) the Chebyshev polynomials Green’s function method \cite{ferreiraCriticalDelocalizationChiral2015}. 

\subsection{The virtual crystal approximation and Fermi’s golden rule}

In the VCA, the various properties of the alloy (e.g., masses, lattice constants, interatomic forces) are determined by the weighted average of the same properties of the end member materials based on the alloy concentration \cite{abelesLatticeThermalConductivity1963}. That is, the alloy is treated as a perfect crystal with modified phonon dispersion and other intrinsic properties based on the end members. The phonon scattering rates due to alloy mass disorder are then calculated via quantum perturbation theory, i.e., Fermi’s golden rule (FGR) \cite{tamuraIsotopeScatteringDispersive1983,tamuraIsotopeScatteringLargewavevector1984}:
\begin{equation}
\frac{1}{\tau^\text{FGR}_{\bm{q}j}} = \frac{\pi}{2 N} \omega^2_{\bm{q}j} \sum_k g_k \sum_{\bm{q}'j'} |\bm{e}_{k,\bm{q}j} \cdot \bm{e}^*_{k,\bm{q}'j'} |^2 \delta (\omega_{\bm{q}j} - \omega_{\bm{q}'j'})
\end{equation}
with $N$ the number of cells in the crystal, $\omega_{\bm{q}j}$ the VCA frequency for phonon with wavevector $\bm{q}$ and polarization $j$, and $\bm{e}_{k,\bm{q}j}$ the normalized VCA eigenvector for the atomic site that the alloyed atoms occupy labeled by $k$. Here, $g_k = f_k (m_k - m_\text{vc})^2/m_\text{vc}^2$ is a mass variance parameter for the $k^\text{th}$ atom type of the alloyed site with $m_k$ being the mass of the $k^\text{th}$ atom type, $m_\text{vc}$ being the VCA averaged mass, and $f_k$ being the concentration of the $k^\text{th}$ atom in the alloy. The total phonon scattering rate can be obtained from Matthiessen’s rule: $\frac{1}{\tau_{\bm{q}j}} = \frac{1}{\tau^\text{FGR}_{\bm{q}j}} + \frac{1}{\tau^\text{3ph}_{\bm{q}j}}$, where $\tau^\text{3ph}_{\bm{q}j}$ is the lifetime limited by anharmonicity, most often computed from 3-phonon interactions. Within the relaxation time approximation, the thermal conductivity ($\kappa$) along a particular direction $x$ can be determined by:
\begin{equation}
\kappa = \frac{1}{\Omega} \sum_{\bm{q}j} C(\omega_{\bm{q}j}) v_{x,\bm{q}j}^2 \tau_{\bm{q}j}
\end{equation}
with $\Omega$ the volume, $C(\omega_{\bm{q}j})$ the heat capacity and $v_{x,\bm{q}j}$ the velocity of the mode in the $x$ direction.

This methodology has been widely used in the literature to predict the thermal conductivity of various semiconducting alloys \cite{tianPhononConductionPbSe2012,lindsayCalculatedTransportProperties2015,pandeyInitioPhononThermal2017,leeLatticeThermalConductivity2014,maIntrinsicThermalConductivities2016}.  
In some studies, the phonon-disorder scattering has been treated using the T-matrix approximation, which takes into account multiple scatterings off a single impurity and is therefore exact in the limit of dilute disorder \cite{kunduRoleLightHeavy2011a,arrigoniFirstprinciplesQuantitativePrediction2018b}. However, for mass-disordered binary alloys at the $50\%$ composition, multiple-occupancy corrections to the T-matrix cancel all contributions from three or more scatterings, leaving only the FGR term (see eq.~(20) of Ref.~\onlinecite{schwartzComparisonAveragetMatrixCoherentPotential1971} and Fig.~9 of Ref.~\onlinecite{elliottTheoryPropertiesRandomly1974}). Thus, we only consider the FGR in the present study. In particular, it has been very successful for Si$_{1-x}$Ge$_x$ \cite{gargRoleDisorderAnharmonicity2011} and Mg$_2$Si$_{1-x}$Sn$_x$ \cite{liThermalConductivityBulk2012}, two cases when a straightforward comparison with experiment has been possible. This success is unexpected since the FGR, as a lowest-order perturbation theory, is supposed to be valid only for weak, dilute disorder. In Si$_{1-x}$Ge$_x$ with maximal disorder ($x = 0.5$), for instance, every atom can be considered a defect, and the mass ratio between Ge and Si is close to 2.5. Therefore, the disorder is neither weak nor dilute and the success of the FGR is confounding.


\subsection{Chebyshev polynomials Green’s function method}

In order to shed light on this puzzle, we confront the FGR with another approach based on the formalism of Green’s functions, the CPGF method, which does not require assumptions of weak or dilute disorder. Indeed, it is a nonperturbative technique that allows for a full treatment of disorder in very large systems. In this way, we can directly evaluate the validity or breakdown of lowest-order perturbation theory by comparing phonon lifetimes and resulting conductivities predicted by the two methods. Moreover, comparisons can be made with frequency resolution and without having to disentangle the effects of anharmonicity, boundary scattering, modeling parameters, and measurement uncertainties as found in experiments. 

The CPGF method \cite{ferreiraCriticalDelocalizationChiral2015} relies on a real-space representation of the phonon Green’s function of a large, disordered supercell (see Appendix~\ref{appendix_green_function}). The retarded Green’s function $\matr{G}(\omega)$ is a frequency-dependent operator that can be defined \cite{mingoClusterScatteringEffects2010,elliottTheoryPropertiesRandomly1974} as
\begin{equation}
\matr{G}(\omega) = \frac{1}{(\omega + i \eta)^2 - \matr{D}}
\end{equation}
in which $\eta$ is a real positive infinitesimal and $\matr{D} = \frac{1}{\sqrt{\matr{M}}} \matr{\Phi} \frac{1}{\sqrt{\matr{M}}}$ is the dynamical matrix of the supercell with $\matr{M}$ the diagonal matrix of the masses and $\matr{\Phi}$ the matrix of the force constants. Since we neglect force constant disorder here, $\matr{\Phi}$ is ordered while $\matr{M}$ and $\matr{D}$ are disordered. Evaluation of the Green’s function allows determination of the phonon density of states (DOS) and spectral function (see Appendix~\ref{appendix_green_function}). In the CPGF method, the phonon Green’s function is expanded on the basis of Chebyshev polynomials, with an efficient iterative evaluation of the successive terms in the expansion (see Appendix~\ref{appendix_CPGF}). Like exact diagonalization, this is a full treatment of disorder in the sense that all diagrams in the self-energy expansion are incorporated, including vertex corrections for two-particle quantities such as the thermal conductivity (see below). Unlike exact diagonalization, the favorable $O(N)$ scaling of CPGF allows for very large supercell sizes of tens of millions of atoms. For such systems, one disorder configuration is sufficient to obtain self-averaged properties. However, the method also scales linearly with the number of harmonic interatomic force constants (IFCs) per atom, rendering it impractical when long-range force constants have to be included.

Phonon-disorder scattering rates can be extracted from the full width at half-maximum of peaks in the spectral function for given VCA phonon modes, in the same manner as in Ref.~\onlinecite{thebaudSuccessBreakdownTmatrix2020} and Ref.~\onlinecite{bouzerarDrasticEffectsVacancies2020}. Large system sizes provide sufficient resolution to probe the lifetimes of low-frequency acoustic modes. Frequency-resolved thermal transport can also be calculated directly using the Kubo formalism (see Ref.~\onlinecite{flickerLatticeThermalConductivity1973,allenThermalConductivityDisordered1993} and Appendix~\ref{appendix_green_function}), which does not assume the presence of well-defined phonon quasiparticles and does not require phonon-disorder lifetimes as inputs. The thermal conductivity in the $x$ direction can be written as 
\begin{equation}
\kappa = \int_0^\infty d\omega \, W_\text{ph} (\omega) \Sigma_\text{ph} (\omega) 
\end{equation}
with $ W_\text{ph} (\omega)$ a normalized $k_B T$ window (see Appendix~\ref{appendix_green_function}) and the phonon transport distribution function (TDF) $ \Sigma_\text{ph} (\omega)$ is:
\begin{align}
\label{kubo}
\Sigma_\text{ph} (\omega) = \frac{\pi k_B ^2 T}{3 \hbar \Omega} \text{Tr} \bigg[ & \text{Im} \matr{G}(\omega + i \frac{\Gamma_\text{3ph}(\omega)}{2}) \matr{S}_x \; \\ 
& \times \text{Im} \matr{G}(\omega + i \frac{\Gamma_\text{3ph}(\omega)}{2}) \matr{S}_x \bigg] \nonumber
\end{align}
where $T$ is the temperature, $\Omega$ is the supercell volume and $\matr{S}_x$ is the heat current operator in the $x$ direction. The trace over the supercell degrees of freedom is computed via an efficient stochastic method (see Appendix~\ref{appendix_CPGF}). A frequency-dependent phonon-phonon scattering rate $\Gamma_\text{3ph}$ is included as an imaginary part in the Green’s function to account for phonon-phonon interactions in a simple way. 

It is informative to compare our present approach to thermal conductivity calculations with theories beyond Peierls-Boltzmann as recently proposed by Simoncelli, Caldarelli, and Isaeva. Simoncelli’s paper establishes a Wigner formalism for complex crystals that separates two terms in the thermal conductivity: diagonal propagative transport described by Boltzmann-Peierls theory and off-diagonal diffusive transport involving phonon modes mixing via mode-dependent broadening due to anharmonicity and disorder \cite{simoncelliUnifiedTheoryThermal2019}. In addition, a very recent paper by Caldarelli extends this formalism to the overdamped regime by using the Green-Kubo formula on the basis of the phonon modes and, neglecting vertex corrections (dressed bubble approximation), writing it in terms of the phonon spectral function \cite{caldarelliManybodyGreenFunction2022}. Isaeva’s approach to thermal transport in disordered systems uses the Green-Kubo formula expressed on the basis of disordered normal modes of the supercell and introduces a constant anharmonic lifetime in the Green’s function for each normal mode \cite{isaevaModelingHeatTransport2019}. Our approach is essentially equivalent to Isaeva’s, but the Green-Kubo formula is implemented in real-space (not using the normal modes) and the anharmonicity is included through a frequency-dependent lifetime, allowing us to study very large systems (tens of millions of atoms) through the CPGF algorithm. In Caldarelli’s language, disorder is treated at the ‘FSF’ level (actually beyond FSF, because vertex corrections are included) and anharmonicity at the ‘LSFA’ level  \cite{caldarelliManybodyGreenFunction2022}.  Thus, our approach can describe systems simultaneously featuring arbitrarily strong disorder (phonons overdamped by disorder, including localization effects) and wave-like tunneling through anharmonic overlap of normal modes.

\subsection{DFT calculations: Mg$_2$Si$_{0.5}$Sn$_{0.5}$}

We focus our study to the case of maximally disordered Mg$_2$Si$_{0.5}$Sn$_{0.5}$ as a representative of the class of semiconducting alloys that are well-described by the FGR. It crystallizes in the antifluorite structure (see Fig.~\ref{fig1}(a)) with Si (\SI{28.09}{amu}) and Sn (\SI{118.71}{amu}) atoms (mass ratio $\sim 4$) randomly placed in a disordered fcc sublattice. The Mg (\SI{24.31}{amu}) atoms occupy the tetrahedral sites, forming an ordered sublattice. 

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth]{fig1}
\caption{(a) Atomic structure of Mg$_2$Si$_{0.5}$Sn$_{0.5}$. (b) Calculated phonon density of states (DOS) from the virtual crystal approximation (VCA, dashed blue curve) and from the Chebyshev polynomials Green’s function method (CPGF, solid black curve). The DOS projected on Mg, Si, and Sn atoms are also shown (green, cyan, and red curves, respectively). (c) Phonon spectral function calculated by the CPGF method and the VCA dispersion (dashed blue curves). The intensity scale is in \si{(rad.THz)^{-1}}.}
\label{fig1}
\end{figure}

The second-order IFCs of the end members Mg$_2$Si and Mg$_2$Sn were computed from density functional perturbation theory \cite{baroniPhononsRelatedCrystal2001}. To make the CPGF calculations tractable, we neglected the long-range Coulomb corrections and applied a \SI{6}{\angstrom} cutoff for atomic interactions before enforcing the acoustic sum rule. The frequency-dependent three-phonon scattering rates $\Gamma_\text{3ph}(\omega) = 1/\tau^\text{3ph}(\omega)$ were determined by averaging the mode-dependent lifetimes $\tau^\text{3ph}_{\bm{q}j}$ over phonon wavevectors and branches.  The $\tau^\text{3ph}_{\bm{q}j}$ were built from harmonic and third-order IFCs previously obtained in Ref.~\onlinecite{liThermalConductivityBulk2012}. Further details can be found in Appendix~\ref{appendix_mg2sisn_model}.

The aforementioned averaging procedure and truncated range of the harmonic IFCs are approximations made so that numerical implementation of the CPGF method is tractable. This introduces some deviation from measured vibrational frequencies and lattice thermal conductivities, though still allows for direct comparison of the FGR and CPGF methods, as the approximations are made in both cases.

\section{Results and discussion}

Fig.~\ref{fig1}(b) gives the total and projected phonon density of states (DOS) of Mg$_2$Si$_{0.5}$Sn$_{0.5}$ calculated with the CPGF method and compared with that from the VCA, while Fig.~\ref{fig1}(c) gives the CPGF phonon spectral function along high-symmetry lines compared with the VCA calculated dispersion. There are significant differences between the VCA and CPGF spectra. Most prominently, there is a frequency range (from \SI{15}{rad.Thz} to \SI{30}{rad.THz}) corresponding to the top of the acoustic branches where the phonon modes are so broadened by disorder that they can no longer be described as well-defined quasiparticles and have a diffuson character in Allen and Feldman's terminology \cite{allenDiffusonsLoconsPropagons1999}. As a result, the lowest Van Hove (VH) singularity at \SI{13}{rad.THz} is partially smoothed out, while the higher VH singularity at \SI{26}{rad.THz} is entirely destroyed. The CPGF spectrum also features two flat branches at the high end of the spectrum (\SI{55}{rad.THz}) that are absent in the  VCA. The vibrational character of these branches is almost entirely on the light atoms Mg and Si, which is consistent with their high frequency. The other optic modes, from \SI{30}{ rad.THz} to \SI{50}{rad.THz}, are dominated by Mg vibrations. Since the Mg sublattice is not disordered, many of these modes are not sensitive to disorder. Indeed, one CPGF branch is identical to its VCA counterpart on the $\Gamma-X-U-\Gamma$ high symmetry path, though not so along the $\Gamma-L$ segment. Not surprisingly, the mode broadenings become vanishingly small close to the $\Gamma$ point as they are protected by translation invariance.  

Figure 2 gives the phonon-disorder contribution to the inverse lifetimes $\Gamma_\text{dis}$ for the acoustic modes from separate FGR and CPGF calculations. These are also compared with the calculated anharmonic scattering rates $\Gamma_\text{3ph}(\omega)$ at \SI{300}{K} in the same frequency range. Both the FGR and the CPGF method predict that $\Gamma_\text{dis}$ follows an $\omega^4$ Rayleigh power law at low frequencies, as expected from point defects \cite{klemensScatteringLowFrequencyLattice1955}. As a result, disorder scattering is smaller than anharmonic scattering for phonons below \SI{6}{rad.THz}, but dominates for acoustic phonons above this frequency.  Surprisingly, FGR and CPGF acoustic linewidths below \SI{13}{rad.THz} are in agreement; however, there is a significant discrepancy between calculated  acoustic linewidths above \SI{13}{rad.THz}: FGR linewidths are roughly 50\% larger than those from the CPGF method. As mentioned in Ref.~\onlinecite{thebaudSuccessBreakdownTmatrix2020}, the unexpected success of the FGR to describe strong disorder at low frequencies is due to the particular form of the mass perturbation term in the phonon spectral function. The Green's function enters its definition~\eqref{def_spectral_function} as 
\begin{equation}
\sqrt{\frac{\matr{M}_\text{vc}}{\matr{M}}} \matr{G}(\omega) \sqrt{\frac{\matr{M}_\text{vc}}{\matr{M}}} = \frac{1}{(\omega + i \eta)^2 - \matr{D}_\text{vc} + \frac{\Delta \matr{M}}{\matr{M}_\text{vc}} \omega^2}
\end{equation}
where $\matr{D}_\text{vc} = \frac{1}{\sqrt{\matr{M}_\text{vc}}} \matr{\Phi} \frac{1}{\sqrt{\matr{M}_\text{vc}}}$ is the VCA dynamical matrix and $\Delta \matr{M} = \matr{M} - \matr{M}_\text{vc}$ is the mass perturbation. The perturbation term is itself proportional to $\omega^2$, therefore this type of disorder is weak at low frequencies. By contrast, an IFC perturbation term would be of the form $\frac{1}{\sqrt{\matr{M}_\text{vc}}} \Delta \matr{\Phi} \frac{1}{\sqrt{\matr{M}_\text{vc}}}$, in which the frequency does not appear explicitly. As a consequence, it is possible to estimate at which frequency the FGR for mass disorder can be expected to fail by a straightforward evaluation of the dominant higher order term in the standard perturbative expansion of the phonon self-energy (see Appendix~\ref{appendix_diagrams} for more details). In the specific case of Mg$_2$Si$_{0.5}$Sn$_{0.5}$, the breakdown of the FGR occurs around the frequency of the lowest VH singularity at \SI{13}{rad.THz}. The sharp upturn in the DOS is fully reflected in the FGR scattering rates, but has much less of an impact on the CPGF lifetimes, leading to an overestimation of the scattering by the FGR.

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig2}
\caption{Acoustic phonon linewidths as a function of angular frequency calculated via CPGF (orange circles) and FGR (blue curve). Room temperature anharmonic phonon-phonon scattering rates are also shown (green curve).}
\label{fig2}
\end{figure}

As will be seen below, the optic phonons contribute significantly to the thermal conductivity of Mg$_2$Si$_{0.5}$Sn$_{0.5}$, so it is worth examining their scattering. Figure \ref{fig3} gives the disorder-induced inverse lifetimes of selected optic phonon modes between \SI{30}{rad.THz} and \SI{50}{rad.THz}, along with the room-temperature anharmonic scattering rates. The $x$-axis corresponds to VCA frequencies for the FGR data and to disorder-renormalized frequencies (shift of the phonon mode due to the presence of disorder) for the CPGF data. The modes with high Mg vibrational weight (>95\%) predominantly sit on the ordered Mg sublattice, thus their disorder scattering rate is comparable to or weaker than the anharmonic scattering rates and have very small frequency renormalization. By contrast, modes with a lower Mg vibrational character have significant weight on the Si/Sn sublattice, so phonon-disorder scattering dominates over anharmonic scattering and the frequencies are substantially renormalized. For most phonons, including many modes with high Mg vibrational character, the FGR breaks down, resulting in large errors compared to the CPGF method. We note that the errors are irregular: the FGR tends to underestimate $\Gamma_\text{dis}$ in the \SIrange{30}{40}{rad.THz} range (though not always), while it tends to overestimate $\Gamma_\text{dis}$ in the \SIrange{40}{50}{rad.THz} range.

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig3}
\caption{Phonon linewidths for selected optic modes calculated by the CPGF method as a function of the renormalized mode angular frequency (circles) and calculated by Fermi’s golden rule (FGR) as a function of the VCA angular frequency (squares). Cyan and orange data correspond to modes with a high Mg vibrational character ($> 95\%$), while blue and red data correspond to modes with lower Mg vibrational character ($ < 90\%$). The room temperature anharmonic phonon-phonon scattering rates are also shown (green curve).}
\label{fig3}
\end{figure}

Next, we investigate how these results propagate into calculations of the thermal conductivity. Figure~\ref{fig4} gives the phonon TDF at \SI{300}{K} and the corresponding accumulated thermal conductivity (inset). The phonon TDF gives the spectral thermal conductivity contributions, which when integrated give the overall thermal conductivity. The curve labeled ‘FGR’ was calculated using the FGR phonon-disorder scattering rates and the Peierls-Boltzmann transport equation within the relaxation time approximation. The curve labeled ‘CPGF’ was calculated using the Kubo formalism (eq.~\eqref{kubo}) directly, without reliance on a phonon quasiparticle picture or the results of Figs.~\ref{fig2} and \ref{fig3}. Note that the FGR and CPGF calculated room temperature thermal conductivities presented here are significantly higher than measurements (Fig.~\ref{fig4} inset). This is partly due to approximations made for numerically tractable CPGF calculations: a range cutoff on the force constants, not including long-range polar corrections, and the averaging of the mode-dependent phonon lifetimes. In the future, this difficulty might be overcome by a modified mixed-space CPGF algorithm in which the VCA part of the supercell dynamical matrix is represented in reciprocal space.

\begin{figure}
\centering
\includegraphics[width=1.0\columnwidth]{fig4}
\caption{Room temperature phonon transport distribution function (TDF) calculated by the Boltzmann transport equation under the relaxation time approximation with the FGR phonon-disorder scattering rates (blue curve) and calculated by the Green-Kubo formalism with the CPGF method (orange circles). The error bars on the circles give an estimate of the uncertainty associated with the stochastic evaluation of the trace (see Appendix~\ref{appendix_CPGF}). The inset gives the corresponding accumulated thermal conductivities. The total measured value from Ref.~\onlinecite{liThermalConductivityBulk2012} is also shown (dashed black line).}
\label{fig4}
\end{figure}

As expected from the spectral function behavior and the acoustic phonon lifetimes, both methods agree below \SI{13}{rad.THz} and disagree between \SI{13}{rad.THz} and \SI{30}{rad.THz}. In particular, the FGR predicted TDF is much lower around \SI{26}{rad.THz} due to the enhanced scattering from the VH singularity in the VCA DOS. More surprisingly, the FGR predicted TDF is quite close to that predicted by the CPGF method for the optic modes (between \SI{30}{rad.THz} and \SI{60}{rad.THz}), which account for almost a third of the total thermal conductivity. This can be explained as follows. First, the FGR $\Gamma_\text{dis}$ is sometimes overestimated and sometimes underestimated in this frequency regime, which leads to some overall compensation. Second, the optic modes with mostly Mg vibrational character are much less scattered by disorder than modes with a higher Si/Sn vibrational character. As a result, modes with high Mg character contribute much more to the thermal conductivity (see Appendix~\ref{appendix_mg2sisn_model}). However, the phonon-disorder linewidths $\Gamma_\text{dis}$ are small relative to anharmonic scattering ($\Gamma_\text{tot} = \Gamma_\text{dis} + \Gamma_\text{3ph}$) for these heat-carrying optic modes with high Mg vibrational character. Thus, the failure of the FGR to predict $\Gamma_\text{dis}$ has little impact on determining their thermal conductivity contributions, yielding a generally accurate prediction of the TDF. 

Given the previous discussion, we expect the FGR to accurately describe thermal transport in mass-disordered alloys when optic modes have only a minor contribution to heat transport. This may be the case in compounds without compositionally distinct sublattices, such as in Si$_{1-x}$Ge$_x$, compounds with an ordered heavy sublattice supporting low-frequency acoustic modes with long mean free paths, or in situations where the optic modes are not thermally populated as occurs at temperatures significantly lower than the upper limit of the phonon spectrum. As in Mg$_2$Si$_{1-x}$Sn$_x$, we expect the FGR to yield good predictions of the thermal conductivity in mass-disordered alloys where there are disorder-insensitive optic modes sitting on an ordered light atom sublattice. However, perturbation theory may fail in compounds possessing both a disordered heavy sublattice (to suppress the conduction of acoustic modes) and a disordered light sublattice (to make disorder relevant at high frequencies). Perhaps more importantly, the FGR is expected to break down with significant IFC disorder in the alloy. Notably, it mostly fails to describe the measured thermal conductivity of III-V alloys (In,Ga)As, (Al,Ga)N, (In,Ga)N, and (In,Al)N \cite{arrigoniFirstprinciplesQuantitativePrediction2018b,maIntrinsicThermalConductivities2016}. 

\section{Summary}

To conclude, we explained the surprising success of lowest order quantum perturbation theory (i.e., Fermi’s golden rule) in describing phonon-disorder scattering and thermal conductivities in predominantly mass-disordered alloys. By confronting this standard computational approach with the non-perturbative Chebyshev polynomials Green’s function method, we demonstrated that the phonon quasiparticle picture breaks down above a certain frequency, but Fermi’s golden rule is valid for low frequency heat-carrying acoustic modes due to the $\frac{\Delta \matr{M}}{\matr{M}_\text{vc}} \omega^2$ form of the mass perturbation. In our reference system Mg$_2$Si$_{0.5}$Sn$_{0.5}$, there is a substantial contribution to the thermal conductivity from optic phonons with high Mg vibrational character, though these modes are scattered more significantly by anharmonicity at room temperature, thus masking the failure of Fermi’s golden rule for these modes. We discussed the conditions for validity of the perturbative predictions of thermal transport. In particular, we expect the perturbative approach to succeed in most mass-disordered alloys, but to fail for alloys with significant force constant disorder and in systems with significant optic mode thermal conductivity contributions derived from a disordered light atom sublattice. 

\section*{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.


