{'fraction_non_alphanumeric': 0.023859404096834264, 'fraction_numerical': 0.004480912476722532, 'mean_word_length': 4.979471120389701, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Probing the validity of the diffuse mismatch model for phonons using atomistic simulations

Due to it's simplicity the diffuse mismatch model (DMM) remains a popular description of phonon transmission across solid-solid boundaries. However, it remains unclear in which situations the DMM should be expected to be a valid model of the underlying physics. Here, its validity is investigated mode-by-mode using a 3-dimensional extension of the frequency domain, perfectly matched layer (FD-PML) method, to study the interface between face-centered cubic solids with interdiffused atoms. While submonolayer levels of interdiffusion are found to increase the number of available modes for transmission, consistent qualitatively with the DMM, we do not find quantitative or qualitative convergence toward the DMM at higher levels of interdiffusion. In particular, contrary to the fundamental assumption of the DMM, modes are not found to lose memory of their initial polarization and wavevector. The transmission coefficients of randomly interdiffused and smoothly-graded interfaces are also compared. While smoothly graded interfaces show strong anti-reflection properties, selection rules still prohibit transmission of many modes, whereas interdiffused interfaces are not subject to such rules and achieve similar thermal interface conductance by transmitting with lower probability but using a wider range of modes.


I. INTRODUCTION
At the interface between two atomically conformal solids, discontinuities in the vibrational properties result in a finite thermal boundary conductance (TBC) for phononic heat transport. The TBC becomes especially important whenever multiple interfaces are closely spaced together or when heat originates from sources which are highly confined by interfaces.
Examples of the former include superlattices, thin films, and high density nanocomposites.
Examples of the latter include finFET transistors, phase-change memory, and the plasmonic transducer of a heat-assisted magnetic recording head 1 . Despite the technological importance of engineering TBC, there are still no widely accepted theories for predicting TBC a priori that are in excellent agreement with experiment, though modern theories based on the diffuse mismatch model (DMM) do generally predict the correct order of magnitude over the entire range of temperature. There are generally two well-developed theoretical frameworks that have been established for modelling TBC: (1) direct non-equilibrium computational simulations in the time-domain by molecular dynamics and (2) phonon-gas treatments which attempt to determine the transport contribution from each phonon mode.
Landry et al 2 have compared TBC predictions obtained using Molecular Dynamics in the harmonic limit with those obtained using both equilibrium and non-equilibrium phonon gas approaches using transmission coefficients (τ ) determined through a Lattice Dynamics based scattering boundary method. Their results only show good agreement between the near-equilibrium phonon gas approach and the MD approach at mass ratios, m r > 2, while for smaller mass ratios they show that the non-equilibrium correction can be applied to the phonon gas treatment. However, in conventional MD approaches, mode-by-mode information is not obtained, which limits the theoretical insights obtainable by such models.
However, recent developments combining modal analysis with MD may reverse this trend 3-6 .
One potential advantage of the MD approach is the natural inclusion of inelastic phonon interactions whenever large displacements are present. This manuscript focuses on the latter approach.
Using the phonon gas approach, the thermal interface conductance can be calculated by performing integration over k-space in the first Brillouin zone as 7 Here, ω is the phonon frequency, v g is the group velocity,n is the unit normal to the interface plane, τ 1→2 is the energy transmission coefficient, and f is the Bose-Einstein statistical distribution. Due to the advent of density functional theory, obtaining the phonon dispersions and group velocities on either side of the interface is typically trivial, but methods of modelling transmission coefficients at solid-solid interfaces remain controversial.
Initial attempts to theoretically describe transmission coefficient included the Acoustic In the DMM, it is assumed that the interface is disordered such that, whatever the details of the interfacial structure and bonding, incident phonons are scattered so intensely that they lose memory of the polarization and direction from which they originated. Using the principle of detailed balance and conservation of energy, the transmission coefficient of an incident phonon depends only on its incident energy and is given by where A numerically discretized version of this equation was given by Reddy 7 , previously. The integration in Eq. 3 is to be performed over k-space surfaces with constant frequency and polarization. For isotropic phonon dispersions in the low frequency limit, where ω = c p ||k|| , the isofrequency surfaces are spheres, and the simpler expressions originally given by Swartz  with probability, where z 0 is the location of ideal interface and where P 1 is the probability of finding an atom of type 1 at location z; is the characteristic interdiffusion length and erf is the error function.
By definition, if an atom type is not 1, then it must be 2. To maintain simplicity, we study interfaces between FCC crystals. Our base case is an idealization of a Si/Ge interface. True Si/Ge takes the diamond crystal structure, but we adopt the approach of Tamura  individual phonon modes were calculated using a uniform k-grid in the Brillouin zone. In order to solve such a large system using minimal computing resources, we used sparse, iterative solvers. The biconjugate gradient method and the generalized minimum residual methods were both found to be effective without preconditioners. Note that because of the roughness of the interface, there were no guaranteed symmetries in the k x -k y plane, though as we will show one does nearly recover cubic symmetries for statistically large simulation domains. Once the transmission coefficient, frequency, and group velocity of each mode are known, the thermal interface conductance can be easily calculated using Eq. 1. For the case of a flat interface, the Lattice Dynamics technique developed by Young and Maris 9 was used to validate our transmission coefficients in a mode-by-mode manner. That technique is specifically designed to study the phonon transmission across the interface between FCC materials with a one atom basis, bonded by nearest neighbor forces, which exactly conforms to our computational case. For those k-points where the lattice dynamics approach is numerically stable, the predicted transmission coefficients were identical to better than 0.01.
Interface conductance for a flat interface obtained for a similar test case was reported by Here G is the thermal interface conductance obtained from simulations and G DM M is the TBC obtained using DMM, a is the lattice parameter of FCC-Si.

III. RESULTS AND DISCUSSIONS
It is instructive to observe the collective behaviour of thermal boundary conductance prior to discussing a mode-by-mode analysis. Figure 3 shows the thermal boundary conductance as a function of the interdiffusion distance, . If the Diffuse mismatch model (DMM) is a valid description of transport for a disordered interface, then one expects that as disorder is increased, G should tend toward the DMM value. Interestingly, Fig 3 indicates that this is not the case.
In the large interdiffusion limit ( = 1.6a), we see that the thermal interface conductance is not approaching the DMM value. For the model Si/Ge interface here, the DMM result for G is greater than that for a perfect interface. In Fig 3, it can be seen that small levels of interdiffusion do increase G to 85% of the DMM value. However, the maximum corresponds to submonolayer interdiffusion distances, . Beyond this level of interdiffusion, the thermal interface conductance is observed to fall monotonically and thus does not converge toward the DMM. In the thick interface limit ( → ∞), one might expect G ∼ −1 (i.e. rather than an interface conductance, one might think of the region as having a bulk thermal conductivity describing the resistance, /k ∼ 1/G). However, Fig. 3  this was attributed to an effective increase in surface area of the interface.
In order to understand the physics behind this behaviour, we analysed mode resolved transmission coefficient data over the entire Brillouin Zone. The phonon dispersion for each material is given in Fig 4 along  several important effects can be seen.

FIG. 4. Dispersion relation for FCC-Si and FCC-Ge along high symmetry lines
One of the most distinctive effects of disorder is to make certain types of mode conversion possible, and consequently to increase the number of available pathways for transmission. In Figs. 5b and 4c, an example of this can be seen in the transverse modes traveling along Γ-X: these have identically zero transmission coefficient for a perfectly smooth interface once the frequency is higher than the maximum frequency of TA phonons on the Ge side. However, submonolayer levels of disorder change this, enabling the high frequency/wavevector TA  the other hand, further increases to the interdiffusion distance (Fig 6c) results in the general decrease in all transmission coefficients with shorter wavelength phonons being disproportionately affected. One might attribute this to scattering, though as we've noted above the interface conductance in Fig. 3 does not appear to be converging to a diffusive limit, G ∼ −1 . Even at the large level of interdiffusion in Fig. 6c, the transmission coefficients to not follow the isofrequency lines. This shows quite clearly that a fundamental assumption of the DMM is incorrect at the mode-by-mode level: a phonon mode's transmission coefficient is not independent of its incident wavevector and polarization, and dependent only on frequency (i.e. memory-less). While one should not expect long wavelength modes to act diffusely since they do not react to local disorder as strongly, Fig. 6c shows that even many short wavelength modes maintain memory of their incident direction.
Finally, we would note that an interdiffused interface is one way to achieve a graded transition between interfaces. Several groups have proposed that a gradation between lay- we constructed a simulation where, rather than having statistically placed masses, a smooth transition in the masses (Fig. 1d) is assumed between materials, over a characteristic distance, . At any atomic site at coordinate, z, the mass of the atom is given deterministically by This gives the same average mass profile as the the interdiffused case, but without the presence of random disorder. Fig. 1 shows the schematic of different interfaces studied in this article. Such a gradation might be impossible to achieve practically since it would require a different homogeneous mass within each layer, but it serves as a useful point for comparison. Figure 7 shows the resulting phonon transmission coefficient throughout the k x -k z plane of the Brillouin zone. It can be seen that for a continuously graded interface, the transmission coefficients are near unity for a wide range of incident phonon wavevectors. Thus, the smooth mass-gradient acts as an effective phonon anti-reflection coating for phonons with sufficiently short wavelength. However, comparing Fig. 7 with Fig. 6c shows that if the mass gradient is instead achieved by interdiffusion over the same characteristic distance and with the same average mass profile, the anti-reflection effect is not reproduced and the transmission coefficients are generally smaller. Thus, these results show that a phononic anti-reflection coatings cannot be practically realised by mere interdiffusion of materials across an extended distance. The resulting thermal interface conductance from an interdiffused interface may even be lower than what could be expected for either an ideal interface or a perfectly diffuse interface, as shown in in Fig. 3 (see = 1.6a ). In addition, while a smooth mass gradient produces high transmission coefficients, the predicted thermal interface conductance (145 MW/m 2 -K) is actually ≈9% lower than the DMM (163 MW/m 2 -K) and only ≈5% higher than for the slightly interdiffused interfaces ( = 0.2a and = 0.4a in Fig. 3). This is because the smoothly graded interface with its perfect in-plane ordering appears to maintain inplane scattering selection rules, preventing mode conversion, and effectively blocking higher wavenumber modes from participating in transport. The interdiffused case = 0.4a (Fig. 6c) has no such restriction and achieves similar thermal interface conductance using a wider range of transmitting modes with lower average transmission coefficient. Thus, the effect of a high transmission coefficient obtained with a graded interface does not significantly outweigh the effect of increased pathways for transmission brought about by interface disorder.

IV. CONCLUSIONS
In summary, the validity of the diffuse mismatch model ( The FD-PML method for phonons has previously been mathematically described and validated for one-dimension chains as well as in two-dimensions for the specific case of determining the scattering cross section of embedded nanoparticles with simple cubic lattice structure and second nearest neighbor interactions 27 . Here a more general mathematical description is given, valid for arbitrary 3D arrangements of atoms. The dynamic equations governing the displacements (ũ) of each atom in the simulation cell (for example, box 2 and 3 in Fig. 2) can be written in the frequency domain for each atom under the harmonic approximation as where M is a diagonal mass matrix given by M = m i δ ij (summation not implied) where m i is the mass associate with the i th degree of freedom (each atom has 3-degrees of freedom), and K are the interatomic force constants given by where Φ is the interatomic potential and u i is the real-space displacement of the i th degree of freedom relative to the equilibrium position. These can be obtained from ab-initio calculation or based on some other description of the interatomic potentials; in the current paper, we used the simplest possible interaction, based on nearest-neighbor bonding via linear spring interactions. The diagonal elements for K must obey k ii = − i =j k ij .
The frequency domain displacements,ũ, are then divided into contributions from an incident wave (which is to be specified before the simulation) and a scattered wave (which is to be found as part of the solution procedure):ũ =ũ inc +ũ scat . Typically, the context for specifying the incident mode is self-evident, but the method provides tremendous flexibility as to how this is accomplished, and several reasonable variations are possible. For example, if the objective of the computation is to simulate the transmission coefficient from one bulk material to another, the incident wave displacements would be specified as a plane wave obeying the dispersion relation in one of the materials, and would be set to zero everywhere the relation isn't expected to hold; so for a Si-Ge interdiffused interface, one could specify an incident wave within some portion of the purely Si region, and zero everywhere else (including the interdiffused zone).
As discussed thus far, the equations for the system would be given as However, the issue of how to handle the outer boundary conditions is also critical. A typical incident plane wave carries energy into a system, and in most physical situations this energy is scattered or reflected by inhomogeneities. This scattered energy must be allowed to escape the system. For reflections from interfaces, the scattered modes are also expected to be infinitely extended in space, which is not possible to simulate numerically.
For scattered modes that do decay radially, the displacements are only zero infinitely far from the scattering center; numerically setting the scattered wave displacements to zero at any finite distance is equivalent to preventing energy from leaving the system. σ is a diagonal matrix of damping coefficients for each degree of freedom. The value of σ can be chosen with considerable freedom, and need not be homogeneous. We use an inhomogeneous σ , whose value for atoms inside the PMLD is given by where σ max is the maximum value of the damping coefficient, L P M L is the total length of the PML layer and the term r i represents the shortest distance of the i th atom from PD-PMLD boundary. For atoms inside the PD the damping coefficient is simply set to zero. Here we have chosen a parabolic distribution of the damping coefficient and a more comprehensive study for obtaining the PML parameters (σ max , L P M L ) is given in Kakodkar et al. 27 . The following empirical relations were used throughout this manuscript to set the PML parameters and were found to result in artificial PML energy reflection coefficients of