Dispersion considerations affecting phonon-mass impurity scattering rates

The ability to accurately predict the thermal conductivity of quantized lattice vibrations, or phonons, is of utmost importance for a wide array of science and engineering applications, especially in nanosystems where scattering rates differ from that of bulk thereby drastically affecting the thermal properties.¹ Since the seminal work by Peierls,² there has been substantial progress in developing models and expressions to accurately account for the effect of various phonon scattering mechanisms on thermal transport.^3–8 However, the accepted forms of the various phononic relaxation mechanisms resulting from these extensive studies are based on low-temperature, harmonic assumptions. Recent work has included anharmonicites in phonon-phonon scattering terms used to solve the phonon Boltzmann equation,^9,10 resulting in excellent agreement between their predictions and measurements for the thermal conductivity of isotopically pure Si, Ge, and diamond. This agreement can be ascribed to, in part, not relying on a fundamentally harmonic, elastic assumption in three phonon scattering events, which are intrinsically anharmonic. This attests to the importance of properly accounting for the underlying phonon physics and crystal properties when deriving phonon scattering rates in solids.

Although three phonon scattering is the dominant thermal resistive mechanism at elevated (above a few Kelvin) temperatures in isotopically pure elements, there are several other phononic scattering mechanisms that affect the thermal conductivity in naturally occurring crystals, alloys, and composite materials.¹¹ For example, the naturally occurring isotopes in Si create phonon scattering events that can lower the thermal conductivity from that of isotopically pure Si by 14% at room temperature.^12,13 This type of phonon-mass impurity scattering process can be used to systematically reduce the thermal conductivity in alloys,^13,14 which has immense importance for thermoelectric applications in which controlled reduction of phonon thermal conductivity is crucial for increasing thermoelectric material efficiency.¹⁵ Clearly, the ability to accurately predict the effects of phonon-mass impurity scattering is of utmost importance for understanding and engineering thermal transport processes in a wide array of material systems.

In this work, the phonon-mass impurity scattering rate is derived from time dependent perturbation theory in a manner applicable to isotropic dispersive crystals. The previously accepted form of the scattering rate for phonon-mass impurity scattering is based on a low temperature assumption; that is, a dispersionless phonon system. In dispersive phonon systems, this assumption breaks down and can lead to errors in the phonon scattering rates. The model for phonon-mass impurity scattering derived in this work predicts a departure from the previously-assumed Rayleigh-like ω⁴ scattering in highly dispersive modes. The mass impurity scattering rate determined from the model derived in this work is in excellent agreement with theoretical calculations. The scattering rate determined from the previously accepted model using low temperature assumptions under predicts the mass impurity scattering rate by over an order of magnitude.

To begin, consider the scattering rate of a phonon mode in an isotropic crystal from a state of wavevector q to that of q^′ due to scattering of an atom with a different mass, given by⁵ 

where Ω is the volume of the primitive cell, ω is the phonon frequency, e is the polarization vector, j is the polarization, |$\bar n$|$n¯$ is the equilibrium phonon distribution, and Γ = ∑_iz_i(ΔM_i/M)² where M is the atomic mass of the host atoms in the crystal and ΔM is the difference in mass between the impurity atoms of species i of concentration z and that of the host atoms. Equation (1) can be recast into the frequency domain by recognizing that |$\omega _j \left( {\bm q} \right) = v_{p,j} \left( {\bm q} \right){\bm q}$|$ωjq=vp,jqq$ and |$d\omega _j \left( {\bm q} \right) = v_{g,j} \left( {\bm q} \right)d{\bm q}$|$dωjq=vg,jqdq$⁠, where v_p,j and v_g,j are the phase and group velocities of a phonon mode of polarization j, yielding

In the Klemens limit, a Debye crystal is assumed, and therefore v_p = v_g. Since the focus of this work is on impurity scattering in an otherwise homogeneous lattice, we consider only energy conserving two phonon processes (i.e., elastic), and therefore, in the Klemens limit, q = q^′. This leads to

As assumed in the Klemens limit, |$\left( {{\bm e}_{{\bm q},j} \cdot {\bm e}_{{\bm q},j} } \right)^2 = 1/3$|$eq,j·eq,j2=1/3$⁠, and Eq. (3) reduces to

Note this differs from the traditionally used Klemens limit which has a factor of 4 in the denominator instead of the factor of 12. The scattering rates derived in this work are given on a per polarization basis, so assuming three identical modes as in the original derivation⁴ yields the better-known form of the Klemens limit for the total scattering rate which is divided by the factor of 4; i.e., a sum of Eq. (4) over j = 3 identical modes.

As apparent from the derivation of Eq. (4), the Klemens limit is not valid for crystals with dispersion. Although this limit is also not valid for a scattering event that does not conserve wavevector – i.e., a scattering event that causes a change in polarization – this work will focus on the regime in which polarization is conserved to examine the effect of dispersion. This gives insight into a subtle inconsistency in the derivation above. Since the derivation of the Klemens limit assumes no change in polarization after the mass impurity scattering event, |$( {{\bm e}_{{\bm q},j} \cdot {\bm e}_{{\bm q},j} } )^2 = 1$|$(eq,j·eq,j)2=1$ and not the typically assumed 1/3 (because of this, the true Klemens limit should actually be a factor of 3 larger than that originally derived, that is, the prefactor should be 3/4 and not 1/4). Therefore, rederiving the scattering rate from Eq. (1), making no assumption about crystal dispersion yields

Note that, in this work, the role of strain fields and changes in binding forces from impurities on phonon scattering is not considered;⁴ only the phonon scattering rate due to an atom of a different mass (hence the term, phonon-mass impurity scattering). The example calculations focus on isotopic variations of single crystalline Si to avoid any large changes in interatomic force constants that could exist in alloys leading to these additional scattering processes. In addition, changes in the temperature trends of thermal conductivity have been observed even in dilute binary alloys indicating additional anharmonicities affecting phonon scattering,^13,14 which is not focused on in this work.

Figure 1 compares phonon-mass impurity relaxation rates predicted from the model derived in this work to the Klemens limit for the longitudinal and transverse branches of Si assuming ΓΩ = 4.3x10^-33 m³ since ∑_i(z_i(ΔM_i/M)²) = 2.17x10^-4 and Ω = 2.0x10^-29 m³ for natural Si.¹⁶ The factor Γ = ∑_iz_i(ΔM_i/M)² was calculated assuming isotopic concentrations of 92.23% of ²⁸Si, 4.67% for ²⁹Si and 3.10% of ³⁰Si (Ref. 11). We use the Si dispersion in the (100) direction calculated from lattice dynamics¹⁷ and fit the data to a 4^th degree polynomial to obtain an analytical expression for |$\omega _j \left( {\bm q} \right)$|$ωjq$⁠. For calculations of the Klemens limit, longitudinal and transverse phonon velocities are taken as 8,970 and 5,332 m s^-1, respectively.

At low frequencies in the regime where the phonon modes are dispersionless, the predictions from both models exhibit the traditionally assumed Rayleigh-like ω⁴ trend. This is intuitive since, in this regime, the phase and group velocities of the phonon modes in Si are nearly identical. However, the models show different trends and values at high frequencies when the phonon modes become very dispersive; this deviation between the models is more apparent in the transverse branch where the dispersionless assumption breaks down about half way into the Brillouin zone in Si. A clear departure from the typical Rayleigh-like frequency trends is observed at high frequencies in the dispersive regime – it is apparent in Fig. 1 that the Klemens limit predicts a lower phonon-mass impurity scattering rate than that predicted using Eq. (5). As quantified in the next section, in naturally occurring silicon, not properly accounting for dispersion can lead to under prediction of the phonon-mass impurity scattering rates by over an order of magnitude.

To quantify the effect of dispersion during phonon-mass impurity scattering, the thermal conductivity of Si is modeled using the scattering rate expressions shown in Fig. 1. Assuming an isotropic medium, the thermal conductivity is given by^7,8

where f is the Bose-Einstein distribution function and T is the temperature. In any material with isotopic impurities, the scattering time is given by |$\tau _j ( {\bm q} ) = ( {\tau _{N,j}^{ - 1} + \tau _{U,j}^{ - 1} + \tau _{I,j}^{ - 1} } )^{ - 1}$|$τj(q)=(τN,j−1+τU,j−1+τI,j−1)−1$⁠, where the Normal and Umklapp scattering rates are given by |$\tau _{N,j}^{ - 1} = A_{N,j} \omega ^2 T[ {1 - \exp [ {3T/\theta } ]} ]$|$τN,j−1=AN,jω2T[1−exp[3T/θ]]$ and |$\tau _{U,j}^{ - 1} = A_{U,j} \omega ^4 T[ {1 - \exp [ {3T/\theta } ]} ]$|$τU,j−1=AU,jω4T[1−exp[3T/θ]]$⁠, respectively,¹⁸ where θ is the Debye temperature and A_N,j and A_U,j are coefficients to be determined. Note that these expressions for Normal and Umklapp three-phonon scattering rates are of a different form than that which is typically assumed;¹¹ these forms for three-phonon scattering were recently determined by Ward and Broido¹⁸ via first-principles calculations of harmonic and anharmonic interatomic force constants using density-functional perturbation theory with an exact numerical solution of the phonon Boltzmann transport equation. From this approach by Ward and Broido, the Normal and Umklapp scattering rates are not restricted to low temperature, dispersionless modes. In fact, these scattering rates were derived specifically for diamond-type lattices, ideal for the study in this work which focuses on scattering rates in Si.

To determine the Normal and Umklapp scattering rates, Eq. (6) is used with only the Normal and Umklapp scattering times and fit with the coefficients A_N,j and A_U,j to data from measurements on isotopically pure Si (²⁸Si) from 100 – 350 K (this high temperature range available in the experimental data is chosen for the fit to avoid any complication from boundary scattering);¹² from this, A_N,L = 6.5x10^-20 s K^-1, A_U,L = 8.5x10^-49 s³ K^-1, A_N,T = 2.8x10^-19 s K^-1, A_U,T = 2.8x10^-48 s³ K^-1. In this fit, θ is slightly adjusted to achieve a better fit over the temperature (θ = 635 K). This is in excellent agreement with the textbook values of the Debye temperature of Si (θ = 625 – 645 K).^16,19 Note these coefficients are different than those determined by Ward and Broido¹⁸ since these coefficients are highly dependent on the dispersion used for the calculations. The fit used to determine these coefficients along with the ²⁸Si data are shown in Fig. 2.

Given the Normal and Umklapp scattering rate, Eq. (6) is then fit to the measured thermal conductivity on natural Si (n-Si).¹² For this fit, the prefactor ΓΩ in Eqs. (4) and (5) is iterated to achieve a best fit using the two different scattering rates. The best fit ΓΩ for the two different cases – Klemens limit (Eq. (4)) and the model derived in this work (Eq. (5)) – are 2.26x10^-32 m³, and 3.76x10^-33 m³, respectively. The best fit ΓΩ using the model derived in this work is in excellent agreement with the theoretically determined value (ΓΩ = 4.3x10^-33 m³, as previously mentioned). The fit with the Klemens limit is nearly an order of magnitude too large indicating that phonon-mass impurity scattering rates are under predicted by the Klemens limit by over an order of magnitude. In other words, if the Klemens limit is used to fit to the thermal conductivity of naturally occurring silicon, the resulting best fit scattering rate for phonon-mass impurity scattering will over predict the actual scattering rate in silicon by over an order of magnitude. This elucidates the importance of using the model derived in this work to properly account for phonon-mass impurity scattering.

To further illustrate the accuracy of this new expression for phonon-mass impurity scattering, calculations of Eq (6) using the theoretically determined ΓΩ for n-Si and the scattering rates determined from ²⁸Si are shown in Fig. 2; that is, Fig. 2 shows calculations for k using the theoretically determined ΓΩ = 4.3x10^-33 m³ (i.e., there is no fitting of Eq. (6) to n-Si). The model derived in this work shows excellent agreement between the predicted κ and measured data for the thermal conductivity of n-Si. Note that the sound agreement between κ predicted using Eq. (5) and the experimental data is obtained with no fitting parameters. The slight difference between the prediction and the experimental data are most likely due to the approximate dispersion assumed in this work. Calculations of Eq. (6) using the Klemens limit are not included since the magnitude of the scattering rate for mass-impurity scattering is so low the thermal conductivity prediction essentially overlaps with the predictions for ²⁸Si (i.e., the effect of mass-impurity scattering predicted from the Klemens limit is negligible compared to the three phonon scattering rates).

In conclusion, a phonon-mass impurity scattering rate is derived in a manner applicable to dispersive crystals. The model for phonon-mass impurity scattering derived in this work predicts a departure from the previously-assumed Rayleigh-like ω⁴ scattering in highly dispersive modes. The mass impurity scattering rate determined from the model derived in this work is in excellent agreement with theoretical calculations where the scattering rate determined from the previously accepted model using assumptions in the Klemens limit under predicts the mass impurity scattering rate by over an order of magnitude.

I am grateful for funding from the National Science Foundation (CBET-1134311).