Disruption of superlattice phonons by interfacial mixing

Molecular dynamics simulations and lattice dynamics calculations are used to study the vibrational modes and thermal transport in Lennard-Jones superlattices with perfect and mixed interfaces. The secondary periodicity of the superlattices leads to a vibrational spectrum (i.e., dispersion relation) that is distinct from the bulk spectra of the constituent materials. The mode eigenvectors of the perfect superlattices are found to be good representations of the majority of the modes in the mixed superlattices for up to 20% interfacial mixing, allowing for extraction of phonon frequencies and lifetimes. Using the frequencies and lifetimes, the in-plane and cross-plane thermal conductivities are predicted using a solution of the Boltzmann transport equation (BTE), with agreement found with predictions from the Green-Kubo method for the perfect superlattices. For the mixed superlattices, the Green-Kubo and BTE predictions agree for the cross-plane direction, where thermal conductivity is dominated by low-frequency modes whose eigenvectors are not affected by the mixing. For the in-plane direction, mid-frequency modes that contribute to thermal transport are disrupted by the mixing, leading to an underprediction of thermal conductivity by the BTE. The results highlight the importance of using a dispersion relation that includes the secondary periodicity when predicting phonon properties in perfect superlattices and emphasize the challenges of estimating the effects of disorder on phonon properties.


I. INTRODUCTION
Superlattices are nanostructures built from periodic alternating layers of dissimilar materials.2][3][4][5] The existence of a secondary periodicity suggests that bulklike phonons will not exist in short-period superlattices.Instead, phonons related to the secondary periodicity, which we will refer to as superlattice phonons, are the vibrational modes of interest.To design a superlattice with a tailored thermal conductivity, a rigorous examination of the interplay between the superlattice period thickness and the interfacial mixing present in any experimental sample, which may disrupt the secondary periodicity, is necessary.This need provides the impetus for an analysis to elucidate the effects of the secondary periodicity and interfacial mixing on the properties of superlattice phonons.
Thermal transport in superlattices has been studied using molecular dynamics (MD) simulations and the Boltzmann transport equation (BTE).Previous MD studies used the equilibrium Green-Kubo (GK) 6 technique, the nonequilibrium direct method, [6][7][8][9][10] or imposed a spatial temperature perturbation and monitored the relaxation to equilibrium 11,12 to predict thermal conductivity.Bottom-up studies using the BTE relied upon the validity of bulk phonon properties in each layer 13,14 and approximations for the specularity and conductance of the internal interfaces. 15While these two approaches can predict trends in cross-plane and in-plane thermal conductivity versus period length, the effects of the secondary periodicity and interfacial mixing on phonon properties cannot be directly obtained.
The effect of interfacial mixing on the cross-plane thermal conductivity, but not on individual phonon modes, of Si/Ge superlattices modeled using the Stillinger-Weber potential was examined by Landry and McGaughey. 10They showed that for perfect Si/Ge superlattices, thermal conductivity decreased with increasing period length before leveling out.For superlattices with interfacial mixing, the thermal conductivity, which was always lower than the corresponding perfect case, increased with increasing period length before leveling out.Savic et al. used Monte Carlo integration and lattice dynamics calculations (harmonic and anharmonic) to predict the phonon properties and cross-plane thermal conductivities of perfect Si/Ge superlattices (i.e., inclusion of the secondary periodicity without interfacial mixing) modeled using the Tersoff potential. 16Their theoretical thermal conductivities overpredicted experimental measurements.Garg et al. used density functional perturbation theory (DFPT) and lattice dynamics calculations (harmonic and anharmonic) to examine phonon properties in perfect Si/Ge superlattices (i.e., inclusion of the secondary periodicity without interfacial mixing). 17Similar to Savic et al., their calculations overpredict experimentally measured cross-plane thermal conductivities.They attributed this discrepancy to the exclusion of interfacial mass-defect scattering in their calculations, which is expected to be present and important in experimental samples.In a follow-up study, Garg and Chen adopted Tamura elastic mass defect scattering theory 18 to modify the DFPT-predicted lifetimes through the Matthiesen rule for Si/Ge superlattices (i.e., inclusion of the secondary periodicity and interfacial mixing). 19For shortperiod superlattices, they found a tenfold decrease in the crossplane thermal conductivity, consistent with the predictions from Landry and McGaughey. 10This same approach was used by Luckyanova et al. 20 in DFPT-driven calculations on GaAs/AlAs superlattices.In a similar way, Hepplestone and Srivastava used lattice dynamics calculations and perturbative methods to study the effects of mass defects and dislocations at the interfaces in GaAs/AlAs and Si/Ge superlattices. 21,22hey found that the phonon scattering in Si/Ge superlattices is dominated by interactions with mass defects and dislocations, but that phonon scattering in GaAs/AlAs superlattices is dominated by phonon-phonon interactions.Thomas and Srivastava used the perturbative method to vary the concentration of mass defects, offering insight into experimental measurements of the thermal conductivity of short-period Si/Ge superlattices. 23hen et al. found, through a combination of experiment and theory, that interfacial mixing is responsible for the reduction in cross-plane thermal conductivity below the homogenous alloy value in Si/Ge superlattices. 24n this paper, we explore the relationship between secondary periodicity, interfacial mixing, and superlattice phonon properties.Molecular dynamics based normal mode decomposition (NMD) is used to predict the full spectrum of phonon properties in unstrained Lennard-Jones (LJ) superlattices with a mass ratio of three for perfect and mixed interfaces.Molecular dynamics has the advantage over reciprocal spacebased lattice dynamics techniques in that disorder can be explicitly included.Furthermore, the trends in superlattice thermal conductivity predicted from MD-based approaches are the same as those obtained from density functional theory based calculations. 9,10,17,19,20he rest of the paper is organized as follows.In Sec.II, the superlattice geometry is defined and the NMD algorithm is reviewed.In Sec.III, superlattice phonon dispersion, phonon lifetimes, and thermal conductivities predicted from GK, NMD, and Tamura theory are presented.In Sec.IV, the results are put in context with the concept of phonon coherence.

A. Superlattice structure and interactions
The superlattices are built by placing atoms on a facecentered cubic lattice, with the two species only differentiated by their masses.The atomic interactions are modeled using the LJ potential, for argon with an energy scale of 1.67 × 10 −21 J, a length scale σ of 3.4 × 10 −10 m, and a mass scale m of 6.63 × 10 −26 kg.This interaction is nonzero for any pair of atoms with an interatomic distance of r less than or equal to the cutoff of 2.5σ .The lighter species has a mass of m and the heavier species has a mass of 3m.The temperature of all simulations is 20 K, for which the zero-pressure lattice constant a is 5.315 Å. 25 We present results in dimensionless LJ units unless otherwise noted.
Each superlattice is identified by its unit cell, which consists of L/2 conventional four-atom unit cells of each species.The unit cell therefore contains 4L atoms.As shown in Fig. 1(a), one period of a 4 × 4 superlattice has eight monolayers (four of each species).The Brillouin zone is a rectangular prism with boundaries at 2π/(La) in the cross-plane direction and 2π/a in the in-plane directions.We consider 2 × 2, 4 × 4, 8 × 8, and 14 × 14 superlattices.Interfacial mixing is introduced to a superlattice by flipping the masses of randomly selected atoms in the monolayers adjacent to the interfaces until the desired concentrations are reached. 10A 4 × 4 superlattice with 80/20 interfacial mixing (the notation corresponds to the concentration of original/foreign species within the monolayers adjacent to the interface) is shown in Fig. 1(b).

B. Thermal conductivity prediction
As in some previous superlattice studies, 16,17,19,20 a solution to the phonon BTE under the single mode relaxation time approximation 26 is used to predict the diagonal components of the thermal conductivity tensor Here, c ph ( κ κ κ ν ) is the volumetric specific heat, v g,α ( κ κ κ ν ) is the component of the group velocity vector in the α direction, and τ ( κ κ κ ν ) is the lifetime of the phonon mode with wave vector κ κ κ and polarization branch denoted by ν.The summation is over the total number of polarization branches, 12L, and the number of unit cells in the MD simulation, N .A quantity of interest for nanostructure design purposes 27 is the phonon mean-free path (MFP), ( κ κ κ ν ), defined as the average distance traveled between scattering events, 26 where v v v g ( κ κ κ ν ) is the group velocity vector.To obtain the required inputs for Eqs. ( 2) and (3), we follow the NMD procedure outlined by McGaughey and Kaviany, 25 Turney et al., 28 and Larkin et al., 29 in which atomic velocities obtained from MD simulation are projected onto the normal mode eigenvectors obtained from harmonic lattice dynamics calculations.The mode-dependent specific heat is set to be k B /V , where k B is the Boltzmann constant and V is the volume of the MD domain, because MD simulations are classical and obey Maxwell-Boltzmann statistics.As temperature increases, the anharmonicity of the atomic interactions causes the specific heat to deviate from k B /V , but the effect is small (less than 3%) for LJ systems at the studied temperature of 20 K. 25 The unit cells for the perfect superlattices, as depicted in Fig. 1(a), are used as inputs to the harmonic lattice dynamics calculations, which are performed using GULP,30 to obtain harmonic frequencies, ω H ( κ κ κ ν ), and eigenvectors, e e e( κ κ κ ν ).The calculations are conducted at the allowed wave vectors, which are specified from Here, b b b α are the cubically orthogonal reciprocal lattice vectors.
To ensure that all wave vectors are in the first Brillouin zone, where n α are integers and N α are constant even integers corresponding to the number of unit cells in the α direction in the MD domain.Based on the convergence study described in the Appendix, we use N x = 8 and N y = N z = 6 (N x = 10 is used for the 2 × 2 superlattices).Group velocities are calculated using finite differencing about a given harmonic frequency from 26,31 The eigenvectors, from harmonic lattice dynamics, and the atomic velocities, from MD simulations performed using LAMMPS with a time step of 4.285 fs, 32 are used as inputs to obtain the trajectories of the time derivative of the normal mode coordinates, q( κ κ κ ν ,t), at time t from In Eq. ( 6), uα ( l b ,t) is the α component of velocity of atom b in the lth unit cell with equilibrium position r 0 ( l 0 ) and e * ( κ κ κ b ν α ) denotes the complex conjugate of the α-component for atom b of the eigenvector for mode ( κ κ κ ν ).While the harmonic lattice dynamics calculations were performed using the unit cells of perfect superlattices, we use the same set of eigenvectors to obtain q( κ κ κ ν ,t) for both perfect and mixed superlattices.The effects of mixing are thus captured through the differences in the atomic velocities between perfect and mixed MD domains.The validity of this assumption for the mixed cases (i.e., projecting onto an approximation of the normal mode) will be assessed in Sec.III.
By taking the Fourier transform of the autocorrelation of Eq. ( 6), the mode kinetic energy power spectrum is obtained: 33 Before evaluating Eq. ( 7), the MD system is equilibrated by velocity rescaling for 10 5 time steps followed by an NVE (constant mass, volume, and total energy) ensemble for 2.5 × 10 5 time steps.The Fourier transform sampling window τ 0 was set to depend upon the superlattice system and the mode frequency.The number of time steps for the Fourier transform sampling window and total number of time steps in the data collection period are given in Table I.The lag between velocity samples is 2 5 time steps, which is sufficient to capture the dynamics of the highest-frequency modes.The power spectrum was averaged over the Fourier transform sampling windows and over the number of independent MD simulations, with the initial atomic velocities sampled from a Gaussian distribution using a random seed (see Table I).
Further averaging was conducted by imposing the symmetry of the irreducible Brillouin zone.
In accordance with anharmonic theory, 34 the power spectrum given by Eq. ( 7) can be approximated to be a Lorentzian function centered at ω A ( κ κ κ ν ) [which is shifted from ω H ( κ κ κ ν ) on average by less than 3% in bulk LJ systems at a temperature of 20 K 28 ], with a full width at half-maximum ( κ κ κ ν ) of the form Here, C 0 is a fitting parameter.The phonon lifetime is 34 Fitting Eq. ( 8) to Eq. ( 7) was done by considering points within three orders of magnitude of the maximum value.The initial guess for ( κ κ κ ν ) was 0.01 and for ω A ( κ κ κ ν ), the frequency at the maximum value of T ( κ κ κ ν ,ω) was used.The GK method, which makes no assumptions about the nature of the vibrational modes responsible for thermal transport, was also used to predict the thermal conductivity for both perfect and mixed superlattices.Landry et al. previously applied the GK method to LJ superlattices, finding agreement with predictions from nonequilibrium MD simulations and an application of the Fourier law. 10 Ten independent MD simulations were performed for each superlattice.For the 2 × 2 and 4 × 4 superlattices, the total simulation length was 10 6 time steps with a correlation window of 5 × 10 4 time steps.For the 8 × 8 and 14 × 14 superlattices, the total simulation length was 10 6 time steps with a correlation window of 10 5 time steps.In order to minimize the uncertainty in the GK predictions, the converged value of the thermal conductivity was specified using the first-avalanche method described by Chen et al. 35

A. Dispersion and participation ratio
Phonon dispersion curves for the perfect 4 × 4 superlattice are shown in Figs.2(a)-2(c).Figure 2(a) corresponds to the [100] (cross-plane) direction, Fig. 2(c) corresponds to the [010] (in-plane) direction, and Fig. 2(b) corresponds to the [111] direction.Frequency gaps emerge at the Brillouin zone boundaries as a consequence of branch folding. 36,37The flat branches for frequencies greater than 15 in Fig. 2(a) indicate low cross-plane group velocities.The branches in Figs.2(b) and 2(c), on the other hand, vary strongly with frequency at most wave vectors.From Eq. ( 2), we note that differences between in-plane and cross-plane components of group velocity are solely responsible for the directional dependence of the thermal conductivity.Dispersion curves for other superlattices show similar features, with more branches and decreasing length of the cross-plane dimension of the Brillouin zone as the period length increases (the length of the Brillouin zone in the in-plane direction remains constant at 2π/a).The superlattice density of states, plotted as the solid blue line in Fig. 2(d), shares the ω 2 dependence at low frequencies with that of the bulk of the heavier species (dotted black line).The high-frequency portion of the superlattice density of states follows similar variations to that of the density of states of the lighter species (dashed orange line).
Spatial localization, which has previously been invoked to explain the period-length dependence of superlattice thermal conductivity, 38 can be estimated by calculating the participation ratio, p( κ κ κ ν ), defined as 39 1 which is a measure of the number of atoms that participate in a given mode.For completely delocalized modes, 1/p ∼ 1/(4L), and for spatially localized modes, 1/p ∼ 1. 39 The inverse participation ratios plotted in Fig. 2(e), where the vertical line corresponds to 1/(4L), indicate that there is no spatial localization in the 4 × 4 superlattice.The frequency dependence of the participation ratio and density of states was not found to vary with superlattice period length.To further investigate the nature of the superlattice vibrational modes, we adopt a modified version of the participation ratio, which we call the partial inverse participation ratio (PIPR).The PIPR is calculated by splitting the components of the eigenvector into portions corresponding to the light (m) and heavy (3m) species, renormalizing the resulting vectors of length 6L, and then applying Eq. (10).The PIPR, plotted for the perfect superlattices in Fig. 3, can be interpreted as a measure of how much one species participates in a given vibrational mode.For the 2 × 2 and 4 × 4 superlattices, the light and heavy PIPRs are comparable to each other across the entire frequency spectrum, indicating that both species participate nearly equally in almost all vibrational modes.For the 8 × 8 and 14 × 14 superlattices, however, the heavier species dominates modes with frequencies less than the maximum frequency of the heavier bulk system (ω 3m,max = 14.3), while the lighter species dominates modes with frequencies greater than ω 3m,max .This observation suggests that as the superlattice period length increases, some vibrational modes may start to resemble the modes of the respective bulk systems.

B. Power spectra
The power spectra, Eq. ( 7), for the nine labeled points in Figs.2(a)-2(c) are plotted in Fig. 4 for both perfect and mixed (80/20 and 60/40) superlattices.While all peaks appear to be Lorentzian centered about a single frequency, there are minor signatures at other frequencies.For perfect superlattices, the amplitude of these minor signatures are two orders of magnitudes smaller than the main peak; the largest being found for mode A. We attribute these minor signatures in the perfect superlattices to the assumption that the normal modes of the harmonic system are representative of the vibrational modes of the true anharmonic system.In mixed superlattices, the intensity of the minor signatures becomes amplified (notably for modes B, E, and H), an indication of the mutual implications of elastic scattering from point defects (random masses with linear springs) and anharmonicity (ordered masses with nonlinear springs) on phonon scattering. 40Modes around a frequency of 12 experience the largest disruptions for all dispersion directions.We attribute this result to the large density of states around this frequency [see Fig. 2(d)], such that there are many channels available for elastic scattering from point defects. 18itting a Lorentzian function [see Eq. ( 8)] to obtain the lifetimes reported in Fig. 4 was deemed suitable for the 80/20 superlattices since the coefficient of determination value 41 for the most affected modes was 0.9.For completeness, the power spectra and lifetimes of 60/40 superlattices are also included in Fig. 4, in which the peaks for modes B, E, and H are further disrupted to a point where the Lorentzian form begins to disappear.This disruption is evidence that perfect superlattice modes are not always good descriptions of modes in the mixed superlattices.Since the same number of modes are present in both perfect and mixed systems, the superlattice phonons that emerge from the secondary periodicity are effectively disrupted as the crystal symmetry is broken by the interfacial mixing.For the remainder of this work, the 80/20 superlattices are used to discuss the effects of interfacial mixing on phonon properties.

C. Lifetimes
The phonon lifetimes are plotted as a function of the harmonic frequencies in Fig. 5.As the period length increases, we maintain the same number of unit cells, such that the total number of atoms increases.Consequently, the minimum frequency decreases and the longest lifetime increases with increasing period length.Overall, the magnitudes of the lifetimes for a given frequency do not vary significantly from one superlattice to another.The lifetimes for all perfect superlattices exhibit ω −2 scaling at low frequencies.This result is consistent with theoretical predictions for phonon-phonon scattering of modes in the Debye regime, where the density of states scales as ω 2 [see Fig. 2(d)]. 42For the 8 × 8 and 14 × 14 superlattices, the perfect systems have two distinct trends, which deviate from the ω −2 scaling at intermediate frequencies and then terminate at the maximum frequencies of the corresponding bulk systems (vertical lines in Fig. 5).The lifetimes in these two regions are comparable to the bulk lifetimes at the corresponding frequencies.This observation is consistent with the emergence of bulklike modes for the longer period superlattices as demonstrated by the PIPRs in Fig. 3.
][45] As mixing is introduced to the superlattices, a ω −4 scaling is observed at intermediate frequencies for the 2 × 2 and 4 × 4 superlattices but not for the 8 × 8 and 14 × 14 superlattices.The complicated dispersions of the superlattices, particularly for the 8 × 8 and 14 × 14 structures, where there is a significant amount of branch folding, are not Debye-like at the intermediate frequencies, leading to a deviation from the ω −4 scaling.The lifetimes of low-frequency modes for all mixed superlattices are not affected and follow a similar ω −2 scaling as seen in the perfect superlattices.
The introduction of interfacial mixing broadens the power spectra (see Fig. 4) and shifts the phonon lifetimes downward (see Fig. 5), particularly at the intermediate and high frequencies.For the 2 × 2 and 4 × 4 mixed superlattices, the lifetimes of some high-frequency modes fall below the Ioffe-Regel limit, τ IR = 2π/ω, where a mode has a lifetime equal to its period of oscillation.The normal modes of a perfect superlattice have a plane-wave structure.Under the assumption that the normal modes of a perfect superlattice are representative of the mixed superlattice, reaching the Ioffe-Regel limit is therefore not an indication of spatial localization but rather of temporal localization.This statement is supported by the inverse participation ratios, which, as shown in Fig. 2(e), are all below the localization limit of unity.Modes that are below the Ioffe-Regel limit can thus be considered to be nonpropagating delocalized modes (i.e., diffusons). 46,47A similar trend in the variation of lifetimes with frequency, dropping below the (red) superlattices.The black line corresponds to the Ioffe-Regel criterion, 2πω −1 .The brown line corresponds to ω −2 scaling.The green line corresponds to ω −4 scaling.Vertical grey lines correspond to the maximum frequency observed in the bulk systems: for the lighter species, ω m,max = 24.7 and for the heavier species, ω 3m,max = 14.3.Lifetime is nondimensionalized using /σ 2 m.
Ioffe-Regel limit at intermediate frequencies and then rising above at higher frequencies, has also been predicted for LJ alloys. 48

Perfect superlattices
The GK method does not require assumptions about the nature of thermal transport or the types of modes present in the superlattices. 49As such, we consider the GK estimate of thermal conductivity to be a benchmark for comparison.With reference to Tables II and III, the trends and magnitudes in cross-plane and in-plane thermal conductivity predictions from NMD and GK for perfect superlattices are in good agreement.Justification for the reported uncertainties is provided in the Appendix.The in-plane thermal conductivity increases with increasing period length, while the cross-plane thermal conductivity first decreases then increases with period length, consistent with previous MD superlattice studies. 6,9,11The inplane thermal conductivity is nearly a factor of two larger than cross-plane thermal conductivity for all superlattices because of the larger group velocities across the entire frequency spectrum for the in-plane direction [see Figs.2(a) and 2(c)].
A thermal circuit model prediction for cross-plane thermal conductivity using the relation 6 is presented in Table II.The boundary resistance for a perfect interface, estimated from the nonequilibrium direct method (R int = 1.4 × 10 −8 Km 2 W −1 ), 50 was combined with bulk thermal conductivities, obtained from NMD, of the lighter material (k m = 1.2 W/m K) and the heavier material (k 3m = 0.7 W/m K) through their respective layer thickness, to obtain an effective thermal conductivity.The thermal circuit model underestimates the thermal conductivity for all superlattices, indicating that the bulk phonons of the constituent species are not an accurate description of superlattice phonons.The relative difference decreases with increasing period length, however, suggesting that this model may become representative of the nature of thermal transport at large enough period lengths.The diffusive limit for the in-plane direction for superlattices with layers of equal length is 6 and is presented in Table III.The in-plane thermal conductivity predictions are expected to approach this limit as the superlattice period length is increased.

Mixed superlattices
From Tables II and III, the in-plane and cross-plane thermal conductivity predictions for the 2 × 2 and 4 × 4 mixed superlattices are reduced from their corresponding perfect systems and approach the 50/50 alloy limit (0.21 W/m K from GK using N x,y,z = 6).The mixed 2 × 2 superlattice loses much of its anisotropy between the in-plane and cross-plane directions as there is mixing in all the atomic layers.From Table II, the cross-plane predictions for mixed superlattices from NMD (using perfect eigenvectors) and GK are in good agreement and follow similar trends, with increasing thermal conductivity with increasing period length.9][20] In this approach, the lifetimes predicted for perfect superlattices, τ perfect ( κ κ κ ν ), are modified using the Matthiesen rule 1 where 1 Here, g 2 (b) is the coupling term for atom b in the unit cell that defines the strength of the mass disordering, where the summation is over the possible species at that atomic position in the unit cell with concentration c μ (b), mass m μ (b), and average mass m(b).Given that there are two atom types in the superlattice unit cell, the lighter atom can be considered to be a mass defect of the heavier portion of the superlattice, and vice versa.g 2 (b) is zero if atom b is unmixed (i.e., for atoms that do not reside within one monolayer of the interface).The δ function in Eq. ( 14) is broadened into a Lorentzian function with width on the order of the frequency level spacing (≈0.1) imposed by the finite size of the systems. 46he cross-plane and in-plane thermal conductivity predictions from NMD and Tamura theory are in good agreement, thus yielding a comparable discrepancy with the GK predictions for the in-plane thermal conductivity [see Table III].While the discrepancy manifests in the in-plane direction, the effects occur at the mode level.This discrepancy is not observed for the cross-plane thermal conductivities because the most disrupted modes in the mixed superlattices (see Figs. 4 and 5) exist at intermediate frequencies that have a near-zero component of group velocity [see Fig. 2(a)].The in-plane difference between NMD and Tamura theory with GK increases with period length because the number of branches with nonzero components of group velocity increases.The thermal conductivity predictions for the mixed superlattices indicate that the use of eigenvectors from perfect superlattices to represent modes in mixed superlattices in the NMD approach is equivalent to the perturbative approximation of Tamura theory.This perturbative approximation is not always valid for mixed superlattices, particularly for modes with a large density of states, where the most disruption is observed.

IV. SUPERLATTICE PHONONS
The term "coherence" has been used in two contexts in relation to phonons.First, to describe the modes that emerge from a secondary periodicity (i.e., in a superlattice 20 or a thin film with a periodic arrangement of holes 51,52 ).Second, to describe the excitation of long-wavelength phonons, usually by femtosecond time-resolved pump-probe techniques, 53,54 that do not carry significant thermal energy and are not found in the MD simulations studied here.The first context and its relation to thermal transport are the focus of this section.In a crystalline solid, the phonons that carry thermal energy belong to the dispersion relation of that material, which depends on the geometry, the harmonic force constants, and the constituent masses.The introduction of a secondary periodicity modifies the dispersion such that modes emerge that do not exist in the composing bulk materials.These nonbulk phonons (in our case, superlattice phonons) propagate and scatter in the periodic structure in a similar manner to what a bulk phonon experiences in a bulk material.
The existence of superlattice phonons has been argued based on modeling and experimental evidence.The predicted minimum in cross-plane thermal conductivity as a function of period length 6,7,9,11,38,55 has been attributed to a transition of thermal transport dominated by superlattice phonons in shortperiod superlattices to thermal transport dominated by bulklike phonons in longer-period superlattices. 56,57A thicknessdependent thermal conductivity of finite-size GaAs/AlAs superlattices measured experimentally was explained in terms of ballistic superlattice phonons. 20y using MD simulations, we do not impose any restrictions on the phonon dynamics.The system is allowed to move through its phase space naturally and thus all classical effects, including coherence and the emergence of superlattice phonons, should be captured.Our results indicate that using the superlattice normal modes captures the physics of thermal transport in a perfect superlattice, emphasizing the importance of using the correct dispersion relation.One cannot use the bulk material phonon properties to predict thermal transport in these non-bulk-like systems.This effect is clearly present for all the systems studied here, as evidenced by the discrepancy between the bulk-based thermal circuit model and the NMD and GK predictions (see Tables II and III).We note, however, that the PIPRs plotted in Fig. 3 suggest that lengthening the superlattice period length generates modes that are becoming localized to the individual layers.
As a means to study how superlattice phonons contribute to thermal conductivity, cross-plane and in-plane thermal conductivity accumulation functions are plotted in Fig. 6.Also plotted are the accumulation functions for the two bulk species.The vertical coordinate of any point on the accumulation function represents the thermal conductivity that comes from phonons with MFPs less than the horizontal coordinate of that point. 58he in-plane accumulation functions of the perfect and mixed superlattices exhibit asymptotic flattening at longer MFPs.The cross-plane accumulation functions of all superlattices contain steplike jumps at longer MFPs, a consequence of the finite resolution of the Brillouin zone. 59These longer MFPs correspond to the low-frequency modes that follow the ω −2 lifetime scaling [see Fig. 5] and have a nonzero group velocity component in the cross-plane direction [see Fig. 2(a)].The significant contributions of longer MFP modes in the cross-plane direction suggest a mechanism for reducing the cross-plane thermal conductivity through the finite size of thin-film superlattices. 20he majority of the contribution (greater than 90%) to crossplane thermal conductivity for all superlattices, perfect and mixed, is from modes with MFPs greater than the superlattice period length.Furthermore, 40% (2 × 2) to 60% (14 × 14) of cross-plane and in-plane conductivity is from modes with MFPs greater than the system size in those directions (N x La for the cross-plane and N y a for the in-plane).For the bulk materials, the entire contribution to thermal conductivity comes from modes with MFPs greater than the lattice constant, with 80% of the contribution from modes with MFPs greater than the system length. 60In the perfect superlattices, MFPs greater than the period length cannot be interpreted any differently than MFPs greater than the lattice constant of a bulk crystalline structure.Based on the accumulation functions and good agreement with GK thermal conductivity predictions, our results indicate that, for the range of superlattices studied here, perfect superlattice phonons and bulk phonons can be treated in the same theoretical frameworks.In this case, the concept of coherent effects can be argued to be purely interpretational and is not required to model phonon transport in these superlattices.
Mixing shifts the in-plane thermal conductivity contribution to shorter MFPs for the 2 × 2 and 4 × 4 superlattices.Differences between perfect and mixed accumulations functions for cross-plane thermal conductivity manifest at intermediate and longer MFPs.These differences for the in-plane and crossplane, however, become increasingly smaller with increasing period length.Furthermore, mixing does not change the range of MFPs.While the perfect and mixed accumulation functions share similar trends, the discrepancy between the NMD and GK in-plane thermal conductivity predictions indicates that some of the mixed modes are fundamentally different than the those of the perfect superlattices.The modes in the mixed superlattices cannot be perfectly described by the modes of perfect superlattices and, therefore, the mixed in-plane accumulation functions are not complete representations of the thermal conductivity.

V. SUMMARY
We used NMD to predict phonon properties at a temperature of 20 K in perfect and mixed LJ superlattices with a mass ratio of three.Differences between in-plane and cross-plane components of group velocity are responsible for the differences between in-plane and cross-plane thermal conductivities [see Figs.2(a)-2(c) and Tables II and III].By fitting the mode power spectra to Lorentzian line shapes, we observe a ω −2 lifetime scaling for low-frequency modes in perfect and mixed superlattices in Fig. 5.In longer period-length perfect superlattices (8 × 8 and 14 × 14), two distinct lifetime trends that terminate at the maximum bulk frequencies were attributed to the separation in the participation of these modes in the layers of the superlattice (see Fig. 3).In mixed superlattices, these trends become smeared and a ω −4 lifetime scaling emerges at intermediate frequencies in the 2 × 2 and 4 × 4 superlattices, indicating elastic mass point defect scattering of Debye-like phonons.We find that interspecies mixing disrupts the secondary periodicity (see Fig. 4) and reduces phonon lifetimes, manifesting in a discrepancy between the in-plane thermal conductivity predictions for mixed superlattices from GK with NMD and Tamura theory [see Table III].The discrepancy is the consequence of the use of eigenvectors from perfect superlattices to describe the modes in mixed superlattices.
Our results demonstrate that further effort is required in order to improve our understanding of the effects of mixing on superlattice phonons.A crucial step is obtaining detailed experimental characterization of the quality of these interfaces, which can then be used as input into the modeling frameworks presented here.Due to the computational cost of DFT-based approaches, MD simulation will continue to be a valuable tool in establishing models of phonon transport in superlattices where the perturbative approach to mixing breaks down.a range of system sizes, plot the inverse of thermal conductivity versus the inverse of the system length, fit a line through the data, and then take the vertical axis intercept as the bulk value. 62This method was not used in previous superlattice studies 16,17,20 and is not used here.The complicated dispersion [see Figs.2(a)-2(c)] does not guarantee that this approach is valid and, as such, understanding size effects in superlattices warrants further work.
FIG. 1. (Color online) Atomic representation of a 4 × 4 superlattice for (a) perfect and (b) 80/20 interfacial mixing cases.Orange atoms have mass m and black atoms have mass 3m.

FIG. 4 .
FIG. 4. (Color online) Power spectra for selected modes of the 4 × 4 perfect and mixed superlattices [indicated by the labeled gray square markers in Figs.2(a)-2(c)].Dark blue corresponds to a perfect superlattice, red corresponds to mixing of 80/20, and light blue corresponds to mixing of 60/40.Reported lifetimes calculated from the fitting of the Lorentzian functions (not shown) are also included.By removing a single MD seed, the average uncertainty in the fitting was determined to be 7.5%.

TABLE II .
Cross-plane thermal conductivity predictions (W/m K).

TABLE III .
In-plane thermal conductivity predictions (W/m K).

TABLE IV .
Size-dependent cross-plane NMD predictions of thermal conductivity (W/m K).