A formalism for calculating the modal contributions to thermal interface conductance

A formalism termed the interface conductance modal analysis (ICMA) method is presented, which allows for calculations of the modal contributions to thermal interface conductance within the context of molecular dynamics (MD) simulations, which inherently include anharmonicity to full order. The eigen modes of vibration in this study are calculated from harmonic lattice dynamics (LD) calculations, however the generality of ICMA formalism also allows for incorporation of anharmonic LD results into the calculations. ICMA can be implemented in both equilibrium MD (EMD) and non-equilibrium MD (NEMD) simulations and both methods show qualitative agreement as validated through study of a simple system of Lennard-Jones solids. The formalism itself is based on a modal decomposition of the heat flow across an interface, which is then substituted into expressions for the conductance either based on EMD or NEMD. As a MD based method, it not only includes anharmonicity, but also naturally includes the atomic level details of the interface quality. The ICMA method now enables more in-depth study of various effects such as temperature, anharmonicity, interdiffusion, roughness, imperfections, dislocations, stress, changes in crystal structure through a single unified model, as it can essentially treat any material or object where the atoms vibrate around equilibrium sites (e.g., ordered or disordered solids and molecules). This formalism therefore serves as an important step forward toward better understanding of heat flow at interfaces.


Introduction
When heat flows across an adjoining interface between two different materials there will be a temperature discontinuity at the interface. The interfacial heat flow can be written as the product of the interface conductance G, which is the inverse of the interfacial resistance, and the temperature difference across the interface T.
D The heat flow across the interface can be carried primarily by electrons in electrically conducting materials, but the contribution from the atomic motions is present in all materials. For solids and rigid molecules, these atomic motions correspond to vibrations around an equilibrium site, which can be decomposed into a series of eigen modes via the lattice dynamics (LD) formalism [1] where the actual anharmonic dynamics manifest as the modes having time varying amplitudes. At a given instant, by knowing the amplitudes of these eigen mode vibrations, one can sum the contributions of all the different eigen modes to recover the vibrations of each atom. The eigen modes are termed phonons in crystalline materials, since they generally correspond to propagating waves such as sound waves. However, in disordered/amorphous solids or molecules many eigen modes may not propagate or resemble the usual definition of a phonon. This has prompted our use of the term eigen mode in the ensuing discussion in order to retain generality, since the formalism presented herein is not restricted to crystalline solids. Within the last decade, techniques have been developed to accurately calculate the individual eigen mode contributions to thermal conductivity from first principles [2][3][4]. These techniques now allow for predictive calculations of the modal contributions to thermal conductivity for materials and nanostructures that have yet to be synthesized [5,6]. Knowing the contributions of specific eigen modes can enable rational design and selection of materials by crafting certain features that will target certain group of modes (i.e., acoustic, optical, longitudinal or transverse), to either inhibit or enhance their transport [7][8][9]. This quantitative capability has improved our ability to predict classical size effects [10] and other nanoscale phenomena [11], which are important effects that ultimately limit heat dissipation in applications such as microelectronics [12]. As feature dimensions in micro/nanoelectronics continue to shrink, interfacial resistance is now much more important and in some cases can become the bottleneck to heat removal [13,14]. Compared to the tremendous advancements in predicting thermal conductivity for crystalline materials over the last decade, predicting interface thermal conductance still faces many challenges. The central issue is that we lack quantitative understanding of the underlying processes that occur at interfaces, because we currently have no way of determining the modal contributions with full inclusion of anharmonicity.
The AMM and DMM take the limit of purely specular and diffuse scattering respectively. Many improvements have been made to these methods [28,[38][39][40][41][42], but neither can include the atomic level detail of the interface quality (e.g., roughness, interatomic diffusion, stress, imperfections etc) and therefore they are fundamentally limited. The development of the AGF method was a major step forward, as it incorporated the atomic level details and also accounts for quantum effects [30,43]. However, most applications of the AGF method have been limited to small system sizes and harmonic interactions, due to analytical complexity and computational expense [30]. Mingo nonetheless has shown that, in principle, anharmonicity can be included in the AGF [30]. To our knowledge, however, the anharmonic AGF has not yet been widely used. The WP method allows the calculation of a mode's transmissivity in the context of a MD simulation, however, it requires that all other modes have zero amplitude [30,[43][44][45]. This effectively corresponds to the T 0 K = limit and therefore is unable to examine the effects of temperature dependent anharmonicity. As a result, the WP method simply reproduces the same results as the harmonic AGF approach [13].
Although exceptions exist [46], the majority of the literature associated with using the DMM has only been able to evaluate elastic scattering interactions where the transmission of a mode's energy across the interface is purely governed by whether or not other modes with similar frequency exist on the other side of the interface [13,31]. It has been argued that inelastic scattering may not be important at room temperature, particularly in systems with nanoscale features [21], while others have argued that anharmonicity can have notable contributions to the thermal interface conductance at high temperatures [28,44,45,[47][48][49]. Until a method that includes anharmonicity/inelastic scattering is tested, it is difficult to conclude whether or not anharmonicity is important. The development and testing of such a method is the central focus of the ensuing discussion.
There have been two major breakthroughs over the past decade that have enabled our new approach. First, McGaughey and Kaviany [50,51] showed that by projecting the instantaneous positions of the atoms in a molecular dynamics (MD) simulation onto the eigen mode solutions (i.e., the mode shapes from LD), one can calculate the instantaneous mode amplitudes as follows where k is the wave vector, n is the mode frequency, p i is the polarization vector for atom i, m i is the mass for atom i, r i is the equilibrium position for atom i, x i is the displacement from the equilibrium for atom i, * denotes the complex conjugate, and N is the total number of atoms. From the mode amplitude X, one can track the equilibrium fluctuations in the mode energy and can directly calculate the relaxation times of individual modes [50,51]. This approach agrees with other independent methods and has contributed greatly to our ability to now calculate individual mode contributions to thermal conductivity from first principles [52][53][54].
The second major development that enabled our new formalism was introduced first by Barrat et al [55] and later by Domingues et al [56]. Both derived an expression for the conductance between any two groups of atoms, based on the fluctuation dissipation theorem. This approach has shown agreement with other methods such as NEMD [57] and now allows for calculation of thermal interface conductance without requiring an externally applied heat flow. It is also important to note that the method introduced by Barrat et al [55] and Domingues et al [56] is general and can be applied to any phase of matter, where the atoms are simply divided into two groups such that the interface can have any arbitrary geometry or condition. Our new formalism introduced in the next section is termed the interface conductance modal analysis (ICMA) method as it combines the essential features of the modal decomposition method introduced by McGaughey and Kaviany [50,51] with the equilibrium MD (EMD) conductance expression derived by Barrat et al [55] and Domingues et al [56]. With the ICMA method, we can determine the modal contributions to conductance at an interface and gauge the effects of anharmonicity/inelastic scattering as well as others (i.e., roughness, interdiffusion, stress, imperfections, dislocations etc) with fidelity.

ICMA formalism
Suppose that we construct a system of two materials where each atom vibrates around an equilibrium site. The two materials can be labeled A and B and consist of N A and N B atoms respectively that can move in threedimensions. When we bring these two systems into contact forming an interface, we can solve the equations of motion in the limit that the interactions are harmonic to obtain a set of 3N=3(N A +N B ) eigen solutions via the LD formalism [1]. We can then write the atomic displacements and velocities as follows, where the summations are over all the eigen modes n in the system, x , i x , i  m i are the displacement from equilibrium, velocity and mass of atom i respectively, and e n i , is the vector (containing 3 Cartesian components) associated with atom i that is contained within the N 3 1 eigen vector matrix that represents the nth solution to the equations of motion. The vector e n i , specifies the direction and magnitude of vibration attributed to the atom i as it participates in the collective motion associated with mode n. X n and X n  are the normal mode coordinates for the position and velocity of mode n, which can be calculated from the inverse of the operations in equations (2) and (3) as where the summations are over all the atoms i in the system and * denotes complex conjugate. For crystalline materials, where there is long range order and periodicity, the eigen vectors e n i , are usually expressed in terms of a wave vector k and branch n in the dispersion diagram as follows e k r e k exp i , 6 where r i,0 is the equilibrium position, and e k i n ( ) is the eigen mode displacement for atom i. This is the more common definition, which is consistent with equation (1). Here, however, we have used a single index n as a label for the eigen solutions instead of the wave vector k and polarization/frequency to identify each mode. This is done to retain generality, since it is not required that the system exhibit long range order/periodicity and therefore not all modes have to correspond to propagating wave solutions. Following the approach pioneered by Barrat et al [55] and Domingues et al [56], we can write the instantaneous heat flow across the interface by simply grouping the atoms into two groups, namely A and B. At each instant, in a microcanonical ensemble, the rate at which energy is transmitted across the boundaries of material A is equal to the instantaneous rate of change of the energy in material B. The Hamiltonian of a system having N atoms can be written as where i F is the potential energy associated with a single atom [52,58]. Using this definition for an individual atom's Hamiltonian, the instantaneous energy exchanged between material A and B is given by Equation (9) is general and can be applied to any model for the atomic interactions that can be represented as a sum of individual atom energies. If only pairwise interactions are present between material A and B, equation (9) can be reduced to where f ij is the pairwise interaction force between the two materials [48,56,59]. For pairwise interactions, it is a natural choice to partition half of the energy in the interaction with atom i and the other half with atom j. From this relation, Domingues et al [56] used the fluctuation-dissipation theorem to show that the correlation in the equilibrium fluctuations of the heat flow are related to the conductance via where G is the thermal conductance between the two materials, A is the interface contact area, k B is the Boltzmann constant, T is the equilibrium temperature of the system, and á ñ  represents the autocorrelation function. This result is in principle similar to the more widely used Green-Kubo expression for thermal conductivity calculations [52]. To simplify the nomenclature, we will use Q instead of Q A B  for interfacial heat flow in the ensuing discussion.
It follows from equation (11) that if one could obtain the modal contributions to the heat flow across the interface, such that, at every instant the sum of those contributions returned the total Q, then G can be rewritten as This would then give the individual contribution of each mode to G as Moreover, the definition of the modal heat flux in equation (12) would allow us to substitute for the total heat flux in both places in equation (11) yielding another expression for G as, with individual contributions from correlations between pairs of modes equal to, This format expresses the interface conductance as the summation of all the auto-correlations and crosscorrelations between eigen mode pairs n and n¢ in the system (G G n n n n , , = å ¢ ¢ ) and would provide new insight into which modes interact. This insight into the mode-mode interactions by ICMA may in some ways be related to the recent elucidation of interfacial scattering interactions using the non-equilibrium Green's function (NEGF) method by Ong and Zhang [60]. However, their application of the NEGF approach did not include anharmonicity and as a result, the ICMA method provides additional detail and potentially deeper insights into the physics than purely harmonic approaches.
The critical step is then to determine Q , n subject to the requirement that Q Q .
In using the substitution of the modal contributions to velocity (equation (3)), we postulate that any quantity that is a direct function of the atomic displacements and/or velocities can be decomposed into its modal contributions via equations (2) and (3). This includes quantities such as temperature [61], pressure [62], entropy [63], heat capacity [64], etc. The issue with equation (9), however, is that it includes the vibrations from both sides of the interface. Using such an expression is inconsistent with the prevailing paradigm used to understand heat flow at solid interfaces, which is based on the phonon gas model (PGM) and the Landauer formalism [65]. The inconsistency arises from the fact that in Landauer formalism although the transmission probabilities are calculated based on the interaction of the atoms at the two sides of interface, the final summation, as a consequence of the principle of detailed balance, however, is only cast in terms of the modes on one side of the interface, where the net heat flow is written as [28,65], In equation (20) the summations are over different polarizations (p) and allowed wave vectors (k x y z , , ) in either material A or B [65]. Here, it is critically important to recognize that current understanding of interfacial heat flow has been, in essence, entirely based on using the modes of an isolated bulk crystal. However, since the modes of an isolated bulk crystal do not involve the atoms of the other material that forms the interface, an expression for the interfacial heat flow that uses the velocities of both sides is conceptually inconsistent with the prevailing view of interfacial transport. Equation (20) is therefore inconsistent with equations (9) and (10), because one cannot project the vibrations of atoms from side B onto the modes of side A or vice versa. Such an operation is ill defined if only the modes associated with one side are used to describe the heat flow. In this sense, the fact that the velocities from both sides of the interface manifest in equation (9), raises interesting questions around whether or not a different set of vibrational modes should be used for describing interfacial heat flow; particularly since previous work has shown that the vibrations of the atoms near an interface can exhibit distinctly different frequency content than the rest of the bulk material [22,66]. For this reason, we perform a LD calculation on the entire combined system, so that equation (18) can be used directly without modification.
It is important to note that positive conductance is not strictly guaranteed by equation (11). Nonetheless, cases of negative conductance have not been observed, as negative thermal conductivity or negative conductance would violate the second law of thermodynamics. It should also be noted, however, that negative conductance contributions from individual modes are not prohibited by this formalism, as it is possible that the integral of the cross-correlations between a mode and all others could be more negative than the auto-correlation is positive. Consequently, in this formalism an individual mode can exhibit a negative conductance contribution. However, because the total conductance is positive in practice, such an occurrence does not violate the second law. Negative conductance contributions, on the other hand, are still inconsistent with the PGM. The PGM, however, is not a fundamental law, rather it is an accurate model for treating phonon-phonon interactions in crystalline bulk materials that has been widely used. There are other classes of materials, such as disordered materials (i.e., alloys and amorphous materials), where the validity of the PGM is questionable and the formalism presented herein merely offers a different and arguably more general viewpoint that does not require the definition of phonon velocities. As a result one can expect to observe negative contributions to the conductance from individual modes, but overall positive conductance from the entire system of modes, in accordance with the second law. It should also be noted that similar to any other property calculated with MD simulations (i.e., thermal conductivity [67]), the values of interface thermal conductance predicted by this formalism will depend on the interatomic potential used. This issue, however, does not affect the generality of the ICMA formalism as the interface conductance can in principle be calculated ab initio. Here, one could construct an accurate interatomic potential as a Taylor expansion of the energy with respect to the atom displacements as has been demonstrated by Esfarjani et al using the direct displacement method [68].

Simulation methodology
As an initial application of the ICMA method, we studied a simple interface between Lennard-Jones (LJ) crystals. The LJ potential is defined based on the following formula [69]: We select equal values of e and s for both materials A and B, which results in equal lattice constants for the two sides. An acoustic mismatch exists at the interface because the mass of the atoms on side B are four times the mass of atoms on side A (m m 4 .
B A = ) Both sides have face-centered-cubic (FCC) lattice structures. In LJ systems the simulations can be performed in LJ dimensionless units [57], but to have the results correspond to a physically meaningful system, we chose LJ parameters in our simulations to be equal to that of argon =´- [70]). Thus, side A represents solid argon (mass m) and side B represents a fictitiously heavier solid argon (mass m 4 . ) In our EMD simulations, the system consists of 3×3×60 fcc unit cells (each side 30 unit cells long), which includes 2160 atoms and 6480 eigen modes. A schematic of the structure can be seen in figure 1(a). We confirmed that increasing the size of cross section does not change the features observed in the results, which is in agreement with other reports [31,59]. We also checked the effect of system size on the calculated modal contributions, which is presented in the results section (figure 2), which shows that increasing the length of the system has a negligible impact on the general trend of the calculated modal contributions. Initially, the system is simulated for 2 ns to reach equilibrium. Then, modal heat flux contributions (Q n ) are recorded for 5 ns in the micro-canonical ensemble. The modal contributions to the heat flow, Q , n are then used in post processing (see equation (14)), which leads to the calculation of modal thermal conductance (G . n ) Examples of calculated correlation functions and the respective evolution of interface thermal conductance are shown in figures S1 in the supplementary materials.
The modal heat flux in equation (19) can also be readily used in NEMD simulations to obtain modal contributions to interface conductance. Here we used G , Q T = D where Q is the time-averaged heat flux in the NEMD simulation, and T D is the temperature difference at the interface, determined by extrapolation of the temperature gradients in each respective material [23,[71][72][73]. The modal contribution to the conductance (G n ) then becomes where the same T D is used for every mode. We used periodic boundary conditions along all Cartesian directions. Hot and cold heat baths were placed in the middle of the left side (mass m) and right side (mass 4 m) of the interface, respectively (see figure 1(b)). A thermal power equal to 6.5 nW is input to the system from the hot bath and extracted from the system from the cold bath. The system is simulated for 2 ns to reach steady state, after which the temperature profile remains constant throughout the structure, which can be seen in figure 1(c). We then averaged the modal heat flux contributions for 2 ns and used equation (22) to calculate the modal contributions to the interface conductance of the interface of two LJ materials. For both EMD and NEMD simulations a time step of 1 fs was considered and the standard deviation in G was less than 5% for five independent initial conditions [74].

Results and discussion
For this system, using EMD simulations, we calculated the thermal interface conductance accumulation function for three different system lengths, which is shown in figure 2. The results show that choosing different lengths for the structure does not change the calculated modal contributions to interface thermal conductance (e.g., <6%). In figure 2 we also identify the frequency above which the heavier side does not have any modes available that can be supported in the bulk of the material. Above this frequency, all of the conductance can be attributed to inelastic/anharmonic processes, since purely harmonic interactions would cause the accumulation to reach 100% by .

max ,B
w Thermal interface conductance accumulation functions for the LJ interface obtained from NEMD simulations are also shown in figure 3 compared with the EMD result. Although small differences exist, the result shows that ICMA yields a similar trend in the normalized modal contributions between EMD and NEMD approaches. Here, we have used normalized thermal interface conductance accumulations to compare EMD and NEMD approaches, because for this case the total values of conductance differ for the two approaches, as has been observed previously [75,76].
The effect of temperature on the calculated modal contributions to thermal interface conductance from EMD simulations is also presented in figure 4. The observed decrease in the values of total thermal interface conductance as the temperature decreases is similar to the trend reported by Duda et al [77].
The temperature dependent calculations presented in figure 4 clearly show the importance of anharmonicity. Unlike thermal conductivity in crystals, which is well known to decrease due to anharmonicity, at interfaces, anharmonicity opens up channels for heat conduction, by allowing modes with different frequencies to interact. This has the net effect of increasing the conductance in this system, which is evidenced by the higher conductance values at higher temperatures and more specifically the increased contribution to the total conductance from modes with frequencies above max ,B w .

Conclusions
A new approach, termed the ICMA method has been presented, which allows one to calculate the modal contributions to thermal interface conductance. The approach is based on the modal decomposition of the instantaneous heat flow across an interface, which can be implemented in either EMD or NEMD simulations. Results from both EMD and NEMD simulations of a LJ solid interface with ICMA show agreement for the normalized thermal interface conductance accumulation versus frequency. This new formalism provides an alternative and arguably more complete picture of the modal contributions to conductance at interfaces since it includes anharmonicity to full order. Many additional comparative studies to understand the effects of temperature, anharmonicity, interdiffusion, roughness, imperfections, dislocations, stress, changes in crystal structure etc are needed. The new ICMA formalism presented herein serves as a critical step forward, which now   [77]. The values are calculated from EMD simulations and are for the 60 unit cell structure. Lattice constants for the simulated structures at 60 K, 30 K, 10 K, 1 K are 5.278 Å, 5.272 Å, 5.267 Å, and 5.264 Å, respectively. enables a high degree of fidelity in studying such effects and is expected to lead to even more new insights on how to rationally design and select materials for applications involving heat flow across interfaces.