Contributions of mass and bond energy difference and interface defects on thermal boundary conductance

The impact of mass and bond energy difference and interface defects on thermal boundary conductance (TBC) is investigated using non-equilibrium molecular dynamics (NEMD) with the Lennard-Jones (L-J) interatomic potential. Results show that the maximum TBC is achieved when the mass and bond energy of two dissimilar materials are matched, although the effective thermal conductivity is not necessarily a maximum due to the contributions of the thermal conductivity of the constituent materials. Mass and bond energy differences result in a mismatch between phonon dispersions, limiting high frequency phonon transport at the interface. This frequency mismatch is defined by a frequency ratio, which is a ratio of the characteristic frequencies of the two materials, presented in the discussion section, and is a reference of the level of phonon dispersion mismatch. Inelastic scattering may result at higher temperatures, especially when there exists a bond energy difference, resulting in strain in the lattice, which w...


I. INTRODUCTION
As the size of electronic, photonic, and phononic devices becomes smaller there is an increasing need to understand thermal transport across interfaces between two different materials.Interfaces are considered to be a thermal transport bottleneck in nanostructured materials and therefore play a critical role in the effective thermal conductivity.In order for next generation devices to perform at the necessary efficiency and desired performance, predicting thermal transport at the interface is crucial. 1,2When heat is transported across the interface between dissimilar materials there exists a resistance to conduction heat transfer creating a temperature discontinuity.This temperature discontinuity can be described by the thermal boundary resistance, or Kapitza resistance, 3 which is the inverse of the thermal boundary conductance (TBC).The TBC is given in equation (1).
where q is the heat flux and ∆T is the difference in temperature across the interface.The first theoretical model to predict the TBC was the acoustic mismatch model (AMM). 4he AMM is based on an assumption of a perfect interface, where incident phonons are treated as plane waves and assumes no scattering occurs at the interface, therefore phonons can only reflect, reflect and mode convert, refract or refract and mode convert.This model is valid only at relatively low temperatures (< 20 K) and perfect planar interfaces.In an attempt to capture scattering at the interface, the diffuse the mismatch model (DMM) was introduced. 5In the DMM phonons are diffusely and elastically scattered at the interface and are re-emitted to either side of the interface.The transmission probability predicted by the DMM depends on the density of states of materials on both sides of the interface and the selection rules that assume different scattering at the interface.To understand the role of scattering, Duda et al. 6 studied different scattering mechanisms, including complete diffuse scattering, partial diffuse scattering, elastic scattering, and inelastic scattering, and showed the dependency of transmission probability on the principle of detailed balance, and found that one of their specific cases produces the DMM developed by Swartz and Pohl. 5 Hopkins and Norris developed the joint frequency diffuse mismatch model (JFDMM) with the assumption of inelastic scattering at the interface, and found that inelastic scattering can contribute to thermal transport at high temperatures. 7Although these models have advanced the understanding of phonon transmission across interfaces by including scattering, a complete and accurate prediction of thermal boundary resistance at higher temperatures and in highly disordered interfaces is still lacking.In practical applications neither the AMM nor the DMM can exactly describe the phonon transport across the interface.
Over the last several decades, since Pickett et al. 8 first used molecular dynamics (MD) simulations, MD has become a widely used technique to study thermal transport across the interface of dissimilar materials, and a variety of different nanostructures, for example superlattices, [9][10][11][12][13][14][15] nanowires, 16,17 and carbon nanotubes. 18,191][22][23][24] Skye and Schelling explored scattering from mass disorder and found it to be a dominant factor in thermal transport, while the bond disorder was found to play a minor role in silicon-germanium alloys using Stillinger-Weber potentials. 21Hu et al. studied the interface of single-crystal solid and amorphous polyethylene and found a strong dependence of Kapitza conductance on interfacial bond energy. 22This result was attributed to the "graded" interface which is created by soft polymer-polymer bonds and stiff Si-Si bonds.In another study using a system of Au-self-assembled monolayer (SAM)-Si, it was found that the TBC at one interface is dependent on the bond energy of the other interface due to proximity of the 2 interfaces. 239][30][31] Research has shown that interface mixing between two solid materials has a significant effect on TBC 28,29 while vacancies have a negligible impact. 30Li and Yang studied the effect of lattice mismatch across the interface of Si and Ge-like materials with defected interfaces and showed that defects reduce the impact of lattice mismatch on TBC. 32lthough previous studies have been performed and found the impact of mass and bond energy between two materials on TBC, they have not distinguished the relative contributions of their impact on TBC.In this study, we perform non-equilibrium molecular dynamics (NEMD) simulations to investigate the relative and combined contributions of mass and bond energy differences on TBC at the interface.Mass and bond energy differences result in a mismatch between phonon dispersions, limiting high frequency phonon transport at the interface.This frequency mismatch is defined by a frequency ratio which is a ratio of the characteristic frequencies of the two materials, presented in the discussion section, and is a reference of the level of the phonon dispersion mismatch.In addition, the impact of defects, such as interdiffusion and vacancies, on TBC is explored.Finally, we propose that it is ultimately the phonon dispersion mismatch and the level of inelastic scattering at the interface that dictate the phonon contribution to TBC.

II. MOLECULAR DYNAMICS METHOD
4][35][36][37][38] NEMD is especially useful for the study of thermal boundary conductance across interfaces between two materials.In this study we use NEMD simulations of argon and argon-like systems using the Lennard-Jones (L-J) interatomic potential to calculate the effective thermal conductivity of the system and the TBC at the interface.All simulations performed in this work use the LAMMPS molecular dynamics package. 39An illustration of the computational domain is shown in figure 1; x is the direction of transport, and only y and z directions have periodic boundary conditions.The computational domain does not impose a periodic boundary in the direction of transport, but rather has a single interface between two baths.The computational domain consists of 36352 atoms where the system size of 128 × 8 × 8 unit cells (UCs) in x, y, and z directions, respectively, with an additional 14 UCs (7 on each side) in the direction of heat flow for fixed position walls, heat source, and sink.6 UCs (3 on each side) are assigned to fixed walls to eliminate evaporation and enforce a known constant heat flux across the system and interface, and 8 UCs (4 on each side) are used as heat baths.
A L-J potential is used to describe the interatomic interactions because it provides a realistic description of interactions between two FCC structures.Only the phonon, or lattice vibration, contribution is considered for thermal energy transport within the solid system. 40The L-J potential is given in equation (2).
where ϵ and σ are parameters that define the energy and length scales of the potential function, respectively, and r is the interatomic separation distance.Parameters for argon are taken from Chen et al., 41 and the interactions between the argon and argon-like materials are defined using the Lorentz-Berthelot mixing rules, ϵ 12 = √ ϵ 1 ϵ 2 and σ 12 = (σ 1 + σ 2 )/2. 42The argon-like material properties are chosen relative to argon and are defined in detail later.
Prior to the NEMD simulations, additional simulations are performed in order to find an appropriate time step.Both pure argon and argon-like materials are simulated to check the validity of the time step for the high frequency materials.The maximum time step is chosen where the stable energy distribution for both cases is observed during the simulation.As a result, a time step of 10 fs is used for all NEMD simulations.The zero pressure NPT (constant number of atoms, pressure, and temperature) ensemble is performed for 200 ps under periodic boundary conditions in order to account for temperature dependent equilibrium atom spacing within the system.In the main NEMD simulation, periodic boundary conditions are used in y and z directions (see figure 1) with the previously calculated lattice constant from the NPT relaxation.The NEMD simulations are run for 200 ps within the NVT (constant number of atoms, volume, and temperature) ensemble to reach an equilibrium energy distribution state at the desired temperature, 10 K, 25 K, and 50 K, followed by 15 ns of the NVE (constant number of atoms, volume, and energy) ensemble while the constant heat flux, between 9.0062 × 10 7 and 1.2609 × 10 8 W/m 2 for different simulations, is applied across the system.The heat flux is imposed across the material by adding and subtracting an equal and prescribed amount of kinetic energy from the bath regions.The imposed heat flux results in a temperature gradient and the thermal conductivity and TBC are calculated from Fourier's law and equation (1), respectively.As the simulations run, the temperature distribution across the material is calculated for each unit cell length from the kinetic energy based on kinetic theory.

A. Impact of mass and bond energy difference
In this work the mass and bond energy ratio are varied to study the role of the individual and combined impact on TBC.The mass of the material, m, is varied first while all of the other parameters are fixed.In the second set of simulations the bond energy, defined by the interatomic potential well depth, ϵ, is the only parameter that is varied.Finally the combined impact of m and ϵ on TBC is studied by varying both parameters at the same time.The TBC from the mass variation is shown in figure 2(a), and the thermal conductivity is shown as an inset.In this simulation the mass of solid argon is fixed on one side of the computational domain, and the mass at the other side is varied.Figure 2(b) represents the temperature profiles, showing two different cases at 10 K across the computational domain during the NEMD simulations.The temperature is linearly distributed from the hot bath to the cold bath and the discontinuity is created at the interface of two dissimilar materials.The thermal conductivity presented here represents an effective thermal conductivity because it is calculated with the full length of the system, however it will be shown as thermal conductivity hereafter.As expected, the maximum TBC and thermal conductivity are obtained when the mass of two dissimilar materials are matched.The TBC decreases and saturates as the mass ratio is varied from unity.
The individual impact of the bond energy on TBC and thermal conductivity is shown in figure 3. The maximum TBC is achieved where the bond energy is the same for two materials at low temperatures, which is similar to the case of mass variation, however, the thermal conductivity is not necessarily maximum at that ratio because the contribution of the individual thermal conductivity of the constituent materials (as is demonstrated by the different temperature gradient on the right hand side of the interface in figure 2(b) for the bond energy ratio case).Additionally, the maximum TBC at 50 K is not observed at the bond energy ratio of 1, as shown in figure 3. The authors believe this behavior indicates the onset of inelastic scattering of high frequency phonons from a combination of strain at the interface and increased vibrational amplitudes at 50 K vs. 10 K.An increase in TBC could be accounted for by inelastic scattering by allowing high frequency phonons (above the maximum acoustic frequencies of the material on the opposite side of the interface) to contribute to transport by scattering into multiple lower frequency phonons.Phonon modal or spectral based techniques are required to confirm this belief.This effect has been observed in another study, and can also result in TBC predictions greater than the DMM and other effects like thermal rectification. 43,44he combined impact of m and ϵ on TBC and thermal conductivity at 10 K is studied and presented in figure 4. Results of the combined impact show that the decrease in TBC is almost equivalent to that from the mass ratio alone.Even though the mass and bond energy differences have their individual contribution to TBC as it is presented above, we observe that the mass difference is the more dominant factor in TBC at the interface of two dissimilar materials.Figure 4 is plotted as a function of maximum phonon frequency ratio because the systems in these simulations are treated as simplified mass-spring systems.We have assumed a simplistic mass-spring calculation to represent the system, but only do so to demonstrate the impact of mass and bond energy on vibrational properties.A characteristic frequency in the mass-spring system is estimated as ω = √ k/m, where k in this system is a function of the potential function near the equilibrium position and will be presented later.The characteristic frequency ratio of dissimilar materials is plotted in figure 5 as a function of mass and bond energy difference showing that the differences in mass and bond energy result in different characteristic frequencies in the two materials.Based on this result, the authors hypothesize that transport across an interface between dissimilar materials is a function of the mismatch in phonon vibrational frequencies and the level of inelastic scattering, though further investigation is needed to confirm the rate of inelastic scattering.

B. Impact of defects
In order to study the impact of defects across the interface on TBC, interdiffusion regions of widths of 4, 10, 50, and 90 UCs, and vacancies are created at the interface.Similar to the previous set of simulations, these simulations focusing on the impact of defects are run at three different temperatures.
As expected the TBC and thermal conductivity at 10 K, shown in figure 6, decrease as the width of the interdiffusion region increases.In figure 6(a), a small change in the characteristic frequency ratio results in a much more significant change in TBC versus the introduction of the interdiffusion region, although a clear trend of decreasing TBC is still observed with increasing the interdiffusion region.Interestingly the thermal conductivity of the alloy, shown in figure 6(c), is greater than that of the system with 90 UCs of interdiffusion region.Presumably this is due to additional resistances in the form of a thermal boundary resistance between the alloy region and the homogeneous regions on each side of the interdiffusion region.Figure 7 shows that the TBC decreases with increasing the interdiffusion region.However as the width of the interdiffusion increases to 50 UC and greater, the TBC appears to become temperature independent.This apparent temperature independence could point to alloy scattering beginning to dominate over multiphonon scattering.This could be similar to what is seen in superlattices and nanowires where weaker temperature dependent thermal conductivity is observed due to increased interface and boundary scattering, respectively. 45,46s with interdiffusion, vacancies are expected to be present and contribute to thermal transport behavior, even in the highest quality single crystalline materials.Thus the impact of vacancy density at the interface is considered by simulating pure argon systems with a variable level of vacancy density.Since only one material is used in these simulations, no interfaces exist, but an impact on thermal transport behavior is established.Figure 8 shows that the thermal conductivity is independent of vacancy density and is, therefore, not expected to contribute significantly to TBC, which agrees with the result from Zhou et al.

IV. DISCUSSION AND CONCLUSION
In this study the stiffness of the "spring" in the mass-spring system, k, is obtained by linearization of the L-J interatomic potential.The linearization is performed using a Taylor series expansion of the L-J interatomic potential at the equilibrium separation distance.The result of the linearization of the L-J potential at the equilibrium separation distance is shown in figure 9, and the expression for k is given in equation ( 3). ) , where ϵ and σ are again the energy and length scales from the L-J potential function, respectively.
In order to study how the mass and bond energy of the material influence the frequency, the phonon dispersion curve is obtained for different variations of argon with changes in mass and bond energy.Figure 10 represents the longitudinal acoustic (LA) phonon branch for these different cases.For the purpose of validation, the phonon dispersion of solid argon is calculated directly from MD simulations using an auxiliary post-process code developed by Kong. 47 This result agrees well with the phonon dispersion of solid argon calculated both analytically 48 and computationally. 49In the case of using a doubled bond energy, the frequencies above the normalized wave vector of 0.5 are greater than the maximum frequency of the solid argon for the same LA branch.Since these frequencies do not exist in the solid argon, we assume that they cannot contribute to TBC, unless an inelastic scattering event has occurred, and is the reason for large reductions in TBC when even small differences in the mass and/or bond energy exist.We hypothesize that at higher temperatures and in disordered interfaces inelastic scattering begins to occur.If this is the case, higher frequency phonons are disproportionately affected resulting in a deviation in the maximum TBC away from matched materials in some cases, as it is observed in figure 3.This study has demonstrated that modifying the mass and bond energy ratios of the constituent interface materials results in changes to the phonon dispersion, and thus a phonon dispersion mismatch between the two different materials and therefore changes in TBC across the interface.Using the NEMD method with L-J interatomic potential the individual and combined impact of mass and bond energy difference as well as the impact of defects on TBC are presented.Results from the study of mass and bond energy difference show that the maximum TBC is achieved when there is no difference in two materials.Some cases deviate from this behavior and the authors believe it is a result of inelastic scattering.The mass and bond energy ratios of 2 result in 97 % and 83 % reductions in TBC at 10 K, respectively, as a result of increased phonon dispersion mismatch.
Defects at the interface also play a significant role on TBC.TBC is dependent on the length of the interdiffusion region as well as temperature.The temperature dependence becomes less significant with increases in the length of interdiffusion region, presumably due to increased alloy scattering.Figure 11 shows that the larger the interdiffusion region is the greater reduction of TBC is observed with small differences in the characteristic frequency.This result differs from that of Stevens et al. 29 due to the difference in the temperature range of the study, resulting in different phonon scattering behavior at the interface.This result also shows some agreement with recent results from Tian et al., 31 although the current study does not include interdiffusion thicknesses thin enough (below 4 UCs) to demonstrate an increased TBC.The impact of vacancies appears to play an insignificant role in thermal transport.
Length Scale of the Potential Function

FIG. 1 .
FIG. 1. Computational domain.Two fixed walls are shown in black, hot and cold baths are shown in red and blue, respectively, and light gray and dark gray represent two different materials, argon and argon-like materials, respectively.

097160- 5 FIG. 2 .FIG. 4 .
FIG. 2. (a) Individual impact of m on TBC.The inset shows the impact on thermal conductivity.(b) Temperature distribution of the domain for the case of mass ratio of 2 and the bond energy ratio of 1.25 at 10 K. Dashed lines do not present the trend but show guide lines to help the eye.

FIG. 5 .
FIG. 5.Frequency ratio as a function of m and ϵ ratio.Dashed lines do not present the trend but show guide lines to help the eye.

FIG. 6 .FIG. 7 .
FIG. 6.(a) Impact of interdiffusion on TBC.(b) Magnification of the result of TBC with interdiffusion.(c) Impact of interdiffusion on thermal conductivity.Dashed lines do not present the trend but show guide lines to help the eye.

FIG. 8 .
FIG.8.Vacancy effects on thermal conductivity.Dashed lines do not present the trend but show guide lines to help the eye.

k = 2
. 11. Impact of interdiffusion on TBC.The inset shows the impact on thermal conductivity.Dashed lines do not present the trend but show guide lines to help the eye.
FIG. 10.Phonon dispersion curves for different cases.Dashed lines do not present the trend but show guide lines to help the eye.All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported license.See: http://creativecommons.org/licenses/by/3.0/Downloaded to IP: 129.123.126.110On:Thu, 08 Oct 2015 22:10:25FIG