Machine learned interatomic potentials for modeling interfacial heat transport in Ge/GaAs

Abstract

Molecular dynamics simulations provide a versatile framework to study interfacial heat transport, but their accuracy remains limited by the accuracy of available interatomic potentials. Past researchers relied on simple analytic potentials and the use of mixing rules to model interface systems, but with minimal justification for their use. Contemporary researchers have access to highly flexible machine learned interatomic potentials (MLIPs), yet these remain understudied within this problem domain. In either case, efforts are needed to rigorously assess and validate interatomic potentials prior to their use in MD simulations of interfacial heat transport. Here, we pursue these efforts while developing interatomic potentials for a model Ge/GaAs system. We first assess the quality of ab initio harmonic force constants (IFC2s) obtained from the interface structure that is used to generate fitting data for MLIPs. We show that these force constants can converge to near bulk-like values within 1 nm away from the interface while also exhibiting a complex relationship across the interface, which likely precludes any successful application of mixing rules. Subsequently, we develop two different MLIPs to model Ge/GaAs using the linear spectral neighborhood analysis potential (SNAP) functional form: one standard SNAP fit to total forces, and another hybrid SNAP fit to anharmonic force components and combined with a harmonic Taylor expansion potential. Each potential is evaluated with respect to bulk thermal properties, interface IFC2s, and stability considerations, providing a valuable comparative analysis that can guide future work.

Introduction

The physics of thermal boundary conductance (TBC) is not well understood, yet it plays an increasingly important role in the thermal management of nanoscale devices [1], [2], [3], [4] and can be leveraged to enhance the performance of applications like thermoelectrics [5], [6] and thermal barrier coatings [7], [8]. Experimental efforts to accurately measure the TBC of interfaces have expanded tremendously over the past couple decades [9], [10], [11]. However, challenges remain in developing sufficiently sensitive metrology for some material systems and in de-convoluting the effects of structural and compositional disorder, strain, and defects near interfaces.

Theoretical efforts are an essential complement to experiments, with the ability to clarify the underlying physics and to enable rational design by systematically exploring the effects of different interfaces and interfacial heterogeneities. One class of theoretical approaches that has been extensively used to try to understand TBC is based on the phonon gas model (PGM) and the Landauer formula, which expresses conductance in terms of phonon particles carrying energy through an interface with some probability, termed the transmission probability [12]. A couple of simple approximations have historically been used to model this transmission probability, namely the acoustic mismatch model [13], [14] and diffuse mismatch model [15], [16], and recently more sophisticated atomistic-level approaches like the atomistic green's function (AGF) [17], [18], [19] and wave packet method [20], [21] have been employed.

Aside from AGF, it is also possible to directly simulate the atomistic dynamics of a system using molecular dynamics. Such an approach is valuable in its generality, since it is capable of studying arbitrarily complicated systems, with full inclusion of temperature-dependent anharmonicity. Additionally, MD enables the study of both statistically averaged properties and detailed temporal dynamics, albeit in a classical limit. With respect to predicting TBC values, interfaces can be simulated under an equilibrium or nonequilibrium framework [22], [23], and methods like interface conductance modal analysis (ICMA) [24], [25], [26] can be used to determine TBC and provide a detailed analysis of the modal contributions.

At the heart of an MD simulation is the choice of the interatomic potential, which models the potential energy surface (PES) of the system and the resulting forces experienced by atoms as a function of their respective positions. Aside from the limitations associated with simulating classical dynamics (i.e., not accounting for quantum dynamics or statistics), the accuracy of an interatomic potential is the primary determinant of the accuracy of an MD simulation and any properties extracted from it. Specifically, in regards to thermal properties, the key values of interest are the interatomic force constants (IFCs), which can be obtained by taking a Taylor expansion (TE) of the PES at a given configurational minima. The lattice dynamics formalism, used to define and compute many thermal properties, generally takes as a primary input the second order (harmonic) and higher order (anharmonic) IFCs of a system [27]. Harmonic IFCs (IFC2s) determine the system's eigenmodes while the anharmonic IFCs govern the interactions between the modes. Developing interatomic potentials that accurately replicate the IFCs relative to some trusted reference is therefore essential for the reliable prediction of thermal properties using MD simulations.

In the past decade, a number of researchers have developed potentials focused on accurately modeling the thermal properties of bulk systems. Generally, ab initio data has been used as a reference, since methods like density functional theory (DFT) have been highly successful in predicting the phonon dispersions and thermal conductivity of bulk systems, the latter via a Boltzmann Transport Equation (BTE) framework [28]. One approach has been to take common analytic potentials and fit them to ab initio forces and IFCs [29], [30], [31], [32], [33], although high accuracy with respect to thermal properties can be difficult to achieve due to the simple functional form of these potentials [34]. In contrast, machine learned interatomic potentials (MLIPs) utilize highly flexible functional forms that can reproduce an ab initio PES with high fidelity, though typically requiring an appreciable amount of training data. The bulk thermal properties of a wide variety of materials systems have been successfully modelled with MLIPs based on linear models [35], [36], [37], Gaussian processes [38], [39], [40], and neural networks [41], [42], [43], [44], [45], [46].

Another approach to modeling bulk systems with high thermal accuracy is to leverage the definition of IFCs and use a Taylor expansion directly as the interatomic potential. Such Taylor expansion potentials (TEPs) by definition have exact force constants if parameterized by the ab initio values and as a result can achieve extremely high force accuracy when incorporating 3rd, 4th and higher order terms [47], [48], [49], [50]. Nevertheless, their widespread adoption among the thermal transport community has been hindered by frequent instability problems [34], [50]. To overcome this, Rohskopf et al. recently developed a method combining a harmonic translationally-invariant TEP (TITEP) with an additional analytic potential or MLIP [51], obtaining high force accuracy, exact harmonic IFCs, and good 3rd order IFC accuracy.

While the development of “vibrationally accurate” interatomic potentials to model bulk systems has seen considerable progress, the situation with interfaces remains relatively unexplored, with unique challenges compared to modeling bulk systems. Such an interatomic potential needs to accurately model three distinct regions -- the two bulk regions away from the interface and the interfacial region itself. The interface is especially challenging, as it exhibits a unique bonding environment where atoms from the respective bulk constituents interact with each other in a chemically distinct way. Here, a fundamental choice needs to be made as to how to handle these bonds. Should the interaction be modeled in a way that explicitly incorporates interface data during fitting, or can a simpler approach be taken, whereby interactions are derived purely from constituent bulk interactions?

Historically, researchers studying the thermal properties of interfaces using MD methods have adopted this latter approach, often out of a lack of better/easier options. The simplest of such approaches is to use a single potential to model both bulk regions and only vary the mass between the two sides [52], [53], [54]. A somewhat more versatile approach that has seen widespread adoption is the use of mixing rules to model interface bonds, in which the parameters of two bulk potentials of the same functional form are mixed according to specific rules. Despite being widely used in the past for many interface systems [55], [56], [57], [58], [59], [60], [61], it is unclear how well potentials with mixed parameters can represent the bonds across an interface, and specifically, their associated force constants.

If a mixing rule strategy cannot capture interface bonds well, it is necessary to develop a potential fit explicitly to interface data (in addition to bulk data). Here one of the most viable strategies is to use a MLIP, leveraging the flexibility of their functional forms to model the different bonding environments of bulk and interface regions with a single potential. Another interesting route would be the development of a TE-style potential, promising near exact IFCs if a method is developed to incorporate bulk and interface IFC sets into a single potential. While a number of prior studies have used MLIPs to model interfaces [41], [62], [63], [64], [65], [66], only one was explicitly focused on interfacial thermal transport. Likewise, to the best of our knowledge, only a single paper has used a TEP to model an interface system, and it did not incorporate explicit ab initio interface IFCs [67].

However, a couple key questions need to be addressed before we can confidently use MLIPs and TE-style potentials fit to ab initio data to study thermal transport at interfaces. For one, it is not clear how large the DFT supercell needs to be in order to model an interface, so that its underlying force constants transition back towards their bulk values. Secondly, it is rarely, if ever, assessed as to how well a given MLIP can approximate the force constants at an interface region, and even then, the IFCs of bulk regions are usually only evaluated indirectly via dispersion and thermal conductivity calculations. Given their importance to thermal transport, the ability of MLIPs and related potentials to reproduce interface IFCs needs to be explicitly evaluated.

In this study, we answer the above questions through an investigation of a Ge/GaAs interface as a model system. Ge/GaAs is specifically of interest because it is nearly lattice matched, rendering it straightforward to study with MD (i.e., minimal defects/dislocations). In addition, the TBC of epitaxial Ge/GaAs samples can be measured using a novel transducerless thermoreflectance techniques [68].

To start, we obtain ab initio harmonic IFCs for the Ge/GaAs interface system and examine the convergence of the IFC2s with distance from the interface. These interface force constants are then used to assess whether mixing rules can provide an effective model for the interfacial bonds. We then turn our focus to methods that explicitly model interface interactions and evaluate their performance with respect to bulk and interface thermal properties. We study a simple MLIP, a spectral neighborhood analysis potential (SNAP) [69], and a hybrid potential inspired by the work of Rohskopf et al. [51], which combines bulk and interface TITEPs with a single anharmonic SNAP potential, directly incorporating IFC2s into a potential that models an interface system.