Examining the Effects of Stiffness and Mass Difference on the Thermal Interface Conductance Between Lennard-Jones Solids

To date, the established methods that describe thermal interface conductance (TIC) and include mode-level dependence have not included anharmonicity. The current intuition is therefore based on the behavior in the harmonic limit, whereby the extent of overlap in the bulk phonon density of states (DoS) (e.g., frequency overlap) dictates the TIC and more frequency overlap leads to higher TIC. Here, we study over 2,000 interfaces described by the Lennard-Jones potential using equilibrium molecular dynamics simulations, whereby we systematically change the mass and stiffness of each side. We show that the trends in TIC do not generally follow that of the bulk phonon DoS overlap, but instead more closely follow the vibrational power spectrum overlap for the interfacial atoms. We then identify the frequency overlap in the interfacial power spectra as an improved descriptor for understanding the qualitative trends in TIC. Although improved, the results show that the basic intuition of frequency overlap is still insufficient to explain all of the features, as the remaining variations are shown to arise from anharmonicity, which is a critical effect to include in interface calculations above cryogenic temperatures.


Kiarash Gordiz 1 & Asegun Henry 1,2
To date, the established methods that describe thermal interface conductance (TIC) and include modelevel dependence have not included anharmonicity. The current intuition is therefore based on the behavior in the harmonic limit, whereby the extent of overlap in the bulk phonon density of states (DoS) (e.g., frequency overlap) dictates the TIC and more frequency overlap leads to higher TIC. Here, we study over 2,000 interfaces described by the Lennard-Jones potential using equilibrium molecular dynamics simulations, whereby we systematically change the mass and stiffness of each side. We show that the trends in TIC do not generally follow that of the bulk phonon DoS overlap, but instead more closely follow the vibrational power spectrum overlap for the interfacial atoms. We then identify the frequency overlap in the interfacial power spectra as an improved descriptor for understanding the qualitative trends in TIC. Although improved, the results show that the basic intuition of frequency overlap is still insufficient to explain all of the features, as the remaining variations are shown to arise from anharmonicity, which is a critical effect to include in interface calculations above cryogenic temperatures.
Thermal transport properties in nanostructures are often significantly impacted by the introduction of interfaces to the system [1][2][3] . Unlike the study of thermal conductivity in crystalline materials that has advanced tremendously in recent years 4 , many questions still remain regarding heat conduction through an interface between dissimilar materials. It is now possible to not only calculate thermal conductivity for existing crystalline materials from different methods 5 including first-principles 4,6 , but one can even predict the thermal conductivity of pure crystalline materials that have yet to be synthesized and experimentally measured 7 . For the study of heat transfer across interfaces, however, the situation is different, as no excellent agreement with experiment over a wide range of temperatures or for various material systems has ever been reported [8][9][10][11][12][13][14][15] .
The picture that is ubiquitously used for thermal interface conductance (TIC) is based on phonon gas model (PGM) 16 , which casts the problem as a transmission problem for phonons impinging on the interface: Here, C and v g are the heat capacity and group velocity of a phonon that impinges on an interface. The phonon energy can then either be transmitted or reflected with a specific transmission probability defined as τ. The PGM treatment of TIC lumps all of the physics associated with the nature of the other material, the quality of the interface and the nature of the chemical bonding/interactions (e.g., everything other than heat capacity and group velocity) into the transmission probability τ. Thus, different physical effects such as interface roughness, inter-diffusion, stress, and dislocations all manifest through τ. It is this single descriptor that encompasses all of the dynamics of the interface and it is the only place where the properties of the other material (e.g., its stiffness, chemical composition, atomic level topography) enter the calculations.
To date, the methods that have been developed to describe τ, or the mode level TIC contributions have been unable to include anharmonicity, with the exception of Hopkins and coworkers for the Diffuse Mismatch Model (DMM) 11 and Mingo's work developing the Atomistic Green's Function (AGF) method 15 . Nonetheless, the only guiding intuition we have is still based on the harmonic limit, where only elastic phonon interactions can occur, which restricts a mode with frequency ω to interactions with other modes of the same frequency. It is then intuitive to reason that if more modes of a given frequency are available on the other side, the probability of phonon transmission is higher and this general trend is reproduced by all of the previous methods that elucidate the mode level contributions 8,17,18 . This basic principle gives rise to the idea that high TIC occurs when a large degree of frequency overlap between the phonon density of states (DoS) exists between the two materials. However, this intuition is only rigorously correct as one approaches absolute zero. At higher temperatures, several new and emerging methods have shown that anharmonic effects (e.g., inelastic interactions) become non-negligible at temperatures as low as 10 K 19,20 . Nonetheless, even though anharmonicity is important, a remaining issue is the fact that the guiding intuition of whether or not frequency overlap is a truly useful descriptor for TIC, has yet to be verified independently i.e., by using an independent method to evaluate TIC that is not based on the same basic framework, such as an experiment or another type of calculation. Here, we used equilibrium molecular dynamics (EMD) to calculate TIC independent of models that yield results consistent with frequency overlap as a means of evaluating the conventional intuition. For simplicity, we used the Lennard-Jones (LJ) potential to systematically change the mass mismatch and stiffness mismatch at the interface, which then changes the degree of frequency overlap. The degree of overlap is then evaluated by calculating the phonon DoS for the atoms in each respective system and the TIC is calculated using DMM to assess if there is qualitative agreement with the EMD results.

Model Details
The systems are constructed from two face-centered-cubic (FCC) LJ solids that are lattice matched, which is a structure that has been used extensively in other interfacial heat transfer studies 19,[21][22][23] . The LJ interatomic potential is taken as, 12 6 where i and j are atom indices and r ij is the distance between atoms i and j. We have three tunable parameters on each side of the interface: the depth of the potential well (ε) (e.g., the stiffness), the location of the potential well minimum σ ( ) / 2 1 6 , and the mass (m) of the atoms. For any set of LJ parameters for the materials 1 and 2 forming the interface (e.g., σ , 1 2 and ε , 1 2 ), we used mixing rules to define the LJ parameters for the cross-species interactions: σ σ σ = ( + )/2 12 1 2 , ε ε ε = ⋅ 12 1 2 24 . As a result, the quantitative values for TIC and qualitative trends for TIC obtained herein are specific to these choices. However, as will be shown later, the conclusions are interesting and instructive for our intuition, and it is straightforward to reason that the conclusions are likely to still have meaning and relevance for more technologically relevant interfaces. Hence, the interest here in such a simplified system is to enable a systematic and independent evaluation of the effects of mass and stiffness mismatch. For other, more realistic interatomic potential choices, the interactions are typically more complicated and most often the stiffness and lattice parameter cannot be changed independently. This can make it difficult to isolate the effect of stiffness itself (e.g., over a large range of stiffness values), without inducing other effects such as lattice mismatch and its associated stress and defects. Here, however, by using LJ solids, we can isolate the effect of stiffness, by keeping the two sides lattice-matched. This is accomplished by simply selecting equal values of σ for the two sides of the interface 22 , which then allows us to independently explore the effects of ε and m on the TIC.
A schematic of the model structures is shown in Fig. 1. We start from the case of having solid argon on both sides of the interface. LJ parameters and mass of argon atom (σ Ar , ε Ar , and m Ar ) are chosen from Ref. 24. Then, keeping all the parameters constant on one side of the interface (e.g., solid Ar), we independently increase ε and m on the other side of the interface from ε Ar and m Ar to ε 10 Ar and m 10 Ar . In this way, we examine the impact of a 10X difference in stiffness and mass. This allows us to generate many interfaces all having solid argon on one side and another solid with different properties on the other. Henceforth, we will refer to these two sides as the "constant" side and the "varying" side respectively. This allows us to study a broad range of distinct interfaces instead of evaluating interfaces by modifying only one parameter as has been done previously 19,[21][22][23] . The results can then be examined as two-dimensional maps of m and ε ( ε − m map) with x-and y-axis representing the changes in m and ε on the varying side and the height (z-axis) equal to the property of interest (e.g., TIC or frequency overlap). Each point on the ε − m map therefore represents a distinct structure and interface. In this study, we changed the values of m and ε along the m and ε axes with a resolution of . m 0 2 Ar and ε . 0 2 Ar , respectively, which resulted in a total of 2116 distinct interfaces. We considered 3 × 3 × 20 (x, y, z) FCC unit cell systems, for both the constant and varying sides. Larger cross sections up to 5 × 5 unit cells and longer systems up to 40 unit cells exhibited less than 5% difference from the 3 × 3 × 20 structures suggesting the 3 × 3 × 20 structures are representative of the limiting behavior of contact between two bulk crystals. The z-direction is chosen along the [100] crystallographic direction and is perpendicular to the planar interface. Periodic boundary conditions are applied in all three dimensions. The LJ cutoff radius is chosen to be equal to σ 3 Ar and a small time step of 0.15 fs is used for all EMD simulations to properly model the dynamics of the highest frequency systems on the ε − m map. To reach equilibrium, we simulated each structure successively in the isothermal-isobaric (NPT) ensemble for 0.5 ns and in canonical (NVT) ensemble for another 0.5 ns. Then, using Eq. 3, we outputted the values of interfacial heat flow under the microcanonical (NVE) ensemble for 3 ns and used Eq. 2 to calculate the TIC value. We averaged the results over 20 ensembles 25 with different initial random velocities, which reduced the statistical variation to less than 4% for all the interfaces on the ε − m map (a total of 42,320 independent EMD simulations). All simulations were conducted using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package 26 and the force routine was modified to include the interfacial heat flux calculation concurrently with the trajectory evolution.
In EMD, the TIC is calculated from the autocorrelation function of the interfacial heat flow, using the fluctuation-dissipation theorem 27 and can be written as 28,29 , where A, T , k B , and Q are the cross-section, simulation temperature, Boltzmann constant, and instantaneous interfacial heat flow across the interface, respectively, and ⋅ ⋅ ⋅ represents the ensemble average. For a pair-wise inter-atomic potential such as LJ, Q can be defined as 28 , where i and j represent atomic indices, A and B represent the two sides of the interface, f ij is the force from atom j on atom i, and  v i and  v j are the velocities of the atoms i and j.

Results and Discussion
We first determined the frequency content of each structure by lattice dynamics (LD) calculations 30  interval under consideration is shown in Fig. 2a. It can be seen that the maximum values of DoS overlap between the two sides are present along the diagonal, which is expected since individual oscillator frequencies scale as ω = / K m, where for a monatomic LJ solid m is the mass of the atoms in the system and ε σ ≈ . / K 57 15 2 is the approximate characteristic spring constant between two first nearest neighbor atoms in the FCC lattice. The systems corresponding to the diagonal in Fig. 2a have the same ε/m ratio as the constant side therefore these systems have the same characteristic frequencies and identical dispersion curves, which leads to the maximum DoS overlap. The ratio of the characteristic frequency of the varying side to the characteristic frequency of the constant side is also shown in Fig. 2b and is another way of visualizing the trends in DoS overlap (Fig. 2a). Figure 2a,b suggest that the highest contrast in vibrational properties between the bulk of the two sides occurs when the stiffness is the same with maximum mass difference (bottom right) and when the mass is the same and but the stiffness difference is maximum (top left). The guiding intuition of frequency overlap would then suggest that these systems should exhibit the lowest values of TIC for the interval considered.
We then used the DMM to calculate the values of TIC for the interfaces on the ε − m map (see Fig. 3). The details associated with using the DMM are available in the original work of Swartz and Pohl 9 and works by Hopkins et al. 10,11,32 . Here, we used the DMM with the following assumptions: (1) isotropic properties for both sides of the interface, (2) Debye model for the phonon dispersion and DoS, and (3) elastic interactions between phonons across the interface. In addition, the TIC integral in the DMM should be performed on the softer side of the interface for all the structures on the ε − m map and we evaluated the TIC values at T = 40 K. It should be mentioned that solid argon has a melting temperature of around 86 K and it is the softest structure among all of the solids considered. Therefore, none of the structures melt at T = 40 K. It should be noted that, even the approaches that can incorporate anharmonic interactions into DMM calculations 11 would not yield qualitatively different predictions, since all of them are based on the vibrational modes calculated from the bulk of the materials. Even if we relax all of our earlier assumptions associated with performing the DMM calculations, the TIC for systems along the diagonal will not change. TIC remains constant because DMM fundamentally utilizes the vibrational information of the bulk of the materials, which are indeed the same for both sides in systems located on the diagonal. Therefore, the DMM in any of its current incarnations does not yield different predictions for drastically different materials that ultimately have the same phonon dispersion.
It should be noted that one should also consider the Acoustic Mismatch Model (AMM) since the interfaces are smooth. Fundamentally, the AMM is arguably more applicable than the DMM, because for such smooth interfaces one also has to take into account momentum conservation in the directions parallel to the interface 33 . Thus we also examine the AMM, as we have calculated the TIC values it predicts on the ε − m map in Fig. 4. The details associated with usage of the AMM can be found in Ref. 34. In Fig. 4 the AMM trend is drastically different from the DMM calculations (Fig. 3), which signifies the fundamental contrast between the two approaches. Interestingly, although the AMM technique is still based on the bulk vibrational information 34 , unlike the DMM approach, it can distinguish the different interfaces located on the diagonal of the ε − m map. The reason is that the definition of transmission coefficient in AMM is based on the acoustic impedance (i.e., ρc, where ρ is the density and c is the speed of sound in the material), and because of the large difference in mass on the ε − m map, there is a large  difference in density. As a result the AMM predicts that the acoustic impedance and therefore the TIC should vary along the diagonal of the ε − m map. The EMD values of TIC are shown in Fig. 5a and they are in qualitative disagreement with the guiding intuition obtained from both frequency overlap analysis (Figs 2 and 3) and the AMM calculations (Fig. 4). This disagreement  raises significant questions around the applicability of the harmonic analyses based on the bulk vibrational characteristics of materials. Clearly, utilizing the bulk vibrational information from sides of the interface in DMM and AMM fails to properly describe the trends in TIC for the systems on the ε − m map. One interesting insight that can in general be obtained from the observed trends in Fig. 5a is that the mass-mismatch between the two sides of the interface seems to have a much stronger effect on the TIC than the stiffness mismatch. The significance of mass mismatch alone on TIC has also been noted in other reports 22 , but the ε − m map also reveals a deviation from this trend at the upper ridge of the map. Instead of using system properties that reflect the bulk vibrations, it was postulated that what may serve as more useful descriptors are the properties that are specifically related to the interfacial region's vibrations. As a first hypothesis, we examined the spectrum of the vibrations near the interface (i.e., interfacial power spectra (IPS)), which have been previously observed to differ from the bulk in some structures 35,36 . We determined the power spectrum of the interfacial atoms at each side of the interface by calculating the averaged power spectrum of velocities 37,38 . The atomic velocities ( ) α= , , v x y z are extracted every 10 time steps (e.g., 1.5 fs), and the power spectrum is calculated using 37,38 , Then, using Eq. 4, the degree of overlap between the interfacial DoS at the two sides of the interface is calculated for each structure and the results are shown in Fig 6. It can be seen that if we use the frequency content of the interfacial atoms, the frequency overlap ( Fig. 6) differs significantly from the bulk vibrational information (Fig. 2a) and the trends are much closer to the TIC values on the ε − m map (Fig. 5a). The fact that IPS presents a better description for heat transfer at the interface of two solids has also been reported by Chen et al. for the heat transport across SiGe interface, where more complex many-body potential is used in the MD simulation 36 . However, a recent study by Alexeev et al. 39 showed that for the solid-liquid interface, the liquid layering effect near the interface, rather than the frequency overlap, plays an important role in determining the interfacial heat transport. Figures 5a and 6 show that the mass ratio plays the dominant role in determining the TIC and degree of interfacial region DoS overlap on the ε − m map. It was then realized that this is because even though the mass and stiffness affect the characteristic frequency in the same way, they do not affect the interfacial interactions in the same way. Most notably, as one moves from one side of the interface to the other, there is a continuous transition in the stiffness, but the transition in mass is discontinuous. Consider for example, when the stiffness on the other side of the interface is 2X larger, the atoms on the softer side, due to the use of mixing rules, experience stiffer long ranged interactions with atoms from the stiffer side that become increasingly strong as one moves toward the interface. At the interface itself, the last plane of softer atoms experiences weaker bond strength on one side, and a much stiffer bond strength on the other. Thus, as one increases the stiffness of the stiffer side of the interface, the increased stiffness also affects the atoms in the softer material and consequently makes them vibrate at higher frequencies. The net effect is that it keeps the frequency content of the atoms in the interfacial region matched for any choice of stiffness. This then leads to a fixed amount of interfacial DoS overlap, regardless of the difference in stiffness. The mass difference, on the other hand, is fundamentally different, because it does not penetrate into the other material the way the strength of interaction does. Since changing the mass difference is independent, one can decrease the frequencies of vibration for one side of the interface without changing the frequencies on the other side, thereby explaining why Figs. 5a and 6 seem to only show trends with respect to mass, but nearly constant TIC vs. stiffness ratio. Prior to this investigation, this behavior would have been non-intuitive and thus the results are highly corrective for our intuition. properties play a key role in the interfacial heat transfer, we then showed that incorporating the vibrational power spectra of atoms at the interface into the calculation of the transmission coefficient in the DMM (i.e., interfacial power spectra approximation (IPSA) approach) significantly improved the predicted TIC values. The remaining differences between the EMD and IPSA predictions are attributed to anharmonic interactions at the interface. Incorporating such effects demands more involved techniques, such as MD based approaches 20 , that allow the inclusion of inelastic coupling between different modes of vibration. 1 56 , σ . 1 9 , and σ . 2 2 , respectively 44 . The proportionality of forces at these separation distances determine the factors used to weight the power spectra of atoms at the interface in the calculation of the total averaged interfacial power spectra. For the LJ interatomic potential and FCC lattice, the weighting factors for the first, second, third, and fourth nearest neighbors were approximately 0.64, 0.31, 0.04, and 0.01, respectively.