Introduction

Van der Waals heterostructures, formed by vertically stacking various monolayer two-dimensional (2D) materials, have emerged as a versatile platform for exploring fundamental physical phenomena and offering great opportunities for fabricating advanced electronic devices1,2. Beyond their remarkable functional capabilities, the intrinsic thermal transport properties of 2D-based heterostructures have attracted considerable attention, as efficient heat dissipation is paramount to ensuring reliable device performance and long-term operational stability3,4. Particularly, given that “the interface is the devil’’, yet it is omnipresent in devices, a pivotal and highly compelling question arises regarding the ideal upper limit of interfacial thermal conductance (ITC) in a pristine, atomically perfect 2D van der Waals heterostructure3,5,6,7,8,9,10,11.

Among the various members of this family, the graphene/hexagonal boron nitride (Gr/h-BN) heterostructure stands out due to its vanishingly slight lattice mismatch and the close structural and mass similarities between Gr and h-BN12,13,14, making it a prominent candidate for achieving the highest ITC among 2D heterostructures. Therefore, the challenge can be reframed as exploring the ideal limit of ITC in Gr/h-BN heterostructure.

However, considerable discrepancies remain between theoretical predictions and experimental measurements of the ITC in Gr/h-BN heterostructure. The experimental values range from approximately 7 to 52 MW m−2 K−1 3,6,7,15, while theoretical predictions span a much broader range, from 5 to 442 MW m−2 K−1 16,17,18,19,20. Due to the limited resolution of thermal measurement techniques such as time-domain thermoreflectance (TDTR) and photothermal Raman3,6,7,15, experimentally isolating the ITC of 2D heterostructures is particularly challenging. In light of this, molecular dynamics (MD) simulations have emerged as the predominant theoretical approach for predicting the ITC values and providing atomic-scale mechanistic insights. This method naturally incorporates full anharmonicity and scales computationally linearly with system size, enabling large-scale and long-duration simulations, making it an ideal approach for considering the slight lattice mismatch between Gr and h-BN and the accompanying surface reconstruction-driven moiré patterns12,13,21,22.

The accuracy of MD simulations, however, depends heavily on the reliability of the interatomic potentials. For the Gr/h-BN heterostructure system, the commonly employed combination of Tersoff23 and Lennard-Jones (LJ)24 potentials has proven inadequate in accurately capturing intralayer and interlayer interactions25,26. While the anisotropic interlayer potential (ILP)26,27,28,29,30, parameterized using state-of-the-art density functional theory (DFT), can simultaneously capture the adhesion and anisotropic repulsive characteristics of the interlayer interaction, the intralayer interactions still rely on empirical potentials such as Tersoff23 and adaptive intermolecular reactive empirical bond order (AIREBO)31. These limitations have prevented the ideal limit of ITC in the Gr/h-BN heterostructure from being fully explored, and the underlying moiré-pattern-related thermal transport mechanisms remain elusive–all of which are essential to provide calibrations for experimental measurements.

In this work, based on the neuroevolution potential (NEP) approach32,33,34, we develop an accurate yet highly efficient machine-learned potential (MLP) specifically tailored for Gr/h-BN heterostructure systems and use it to reconstruct the experimentally observed moiré patterns and periods accurately. Employing NEP-driven nonequilibrium MD (NEMD) simulations on a large Gr/h-BN system comprising over 300,000 atoms, augmented with quantum-statistical corrections to the classical MD outcomes, the ITC values in quantitative agreement with the experimental ones are achieved. Such an agreement, combined with the phonon density of state (PDOS) overlap criterion, enables us to report the ideal limit range of ITC for 2D heterostructures at room temperature. The quantum-statistical corrections derived from spectral thermal conductance also elucidate a moderate temperature dependence of ITC in the Gr/h-BN heterostructure.

Further, in the Gr/h-BN heterostructure, we identify a novel physical picture of stacking-sequence-dependent ITC (ABp  > DW  > AB  > AA). To elucidate the atomic-level mechanisms, van der Waals potential energy, interlayer distances, and energy-weighted phonon dispersion are systematically calculated. Their results collectively reveal that stacking sequences with larger interlayer distances are unfavorable for out-of-plane phonon energy transmission, manifested by smaller group velocities and phonon lifetimes, thereby driving the stacking-sequence-dependent ITC. This phenomenon can further be associated with moiré patterns, likely leading to moiré-pattern-dependent ITC.

Results

An accurate and efficient NEP model for Gr/h-BN heterostructure

MLP represents an emerging and powerful technology that leverages the robust fitting capabilities of machine learning to accurately describe potential energy surfaces (PESs) across the entire chemical and structural space provided by quantum-mechanical calculations, such as DFT35. A well-trained MLP can predict physical quantities and drive atomistic simulations with computational efficiency approaching that of an empirical potential, while maintaining the high fidelity of DFT and far beyond its spatial and temporal scales. Due to its high computational efficiency, the NEP framework has become a popular MLP method widely employed in thermal transport studies36,37,38,39,40,41, and is therefore utilized in this work.

Constructing a high-quality and diverse training dataset is a pivotal step in training the NEP model. As shown in the sketch-map representation in Fig. 1a, the dataset encompasses various configurations essential for building the NEP model, such as Gr and h-BN monomers with different compositions, and the Gr/h-BN heterostructure with varying numbers of layers. Each point in the sketch map represents a distinct structure, with similar structures clearly grouped and dissimilar ones well separated. The spatial positions are determined through principal component analysis (PCA) of learned descriptors of the local atomic environments42,43, highlighting the breadth of the configurational space covered by the dataset. Detailed information on the creation and composition of the dataset can be found in the Methods section.

Fig. 1: Training dataset and results for constructing the NEP model.
Fig. 1: Training dataset and results for constructing the NEP model.The alternative text for this image may have been generated using AI.
Full size image

a A sketch map consisting of 4099 structures from the Gr/h-BN training dataset. Each colored dot represents an individual atomic configuration, with color indicating the energy associated with each structure. Snapshots of key configurations are visualized to highlight the diversity of the dataset. b Comparison between the NEP predictions and DFT reference data of energy, force, and stress for the training dataset. In the parity plot, the color intensity visualizes the distribution and density of the training data.

The initial structures were subjected to single-point DFT calculations to label the total energy, force, and stress. We employed the Perdew-Burke-Ernzerhof (PBE) functional44 to describe the exchange-correlation interactions in the structure, and utilized the advanced many-body dispersion (MBD) framework45,46 to capture long-range dispersion interactions. Thus, the DFT reference data are labeled as PBE+MBD (for details, see Methods section). Figure 1b shows a well-trained NEP model for the Gr/h-BN heterostructure (see Methods section for details on training and hyperparameters), which achieves relatively high accuracy. The root mean square errors (RMSEs) for energy, forces, and stress relative to the DFT references are 0.8 meV/atom, 51.2 meV/Å, and 0.1 GPa, respectively.

To comprehensively validate the performance of the NEP model, we conducted a series of benchmark tests. As depicted in Fig. 2a, the NEP model accurately captures the PBE+MBD binding energy curves for different stacking types in the Gr/h-BN heterostructure, with ABp emerging as the most energetically favorable stacking configuration. Notably, the NEP model demonstrates a reliable prediction of domain wall (DW) stacking configurations with varying layer spacings without including them in the training set, highlighting its outstanding interpolation capability. The NEP model also accurately predicts the binding energy curves for different stacking types in both bulk and bilayer Gr and h-BN systems (see Fig. S1), owing to its nearly error-free fitting of DFT energies (see Fig. 1b). For the shallow sliding PES of the Gr/h-BN heterostructure, the difference compared to PBE+MBD is less than 0.3 meV/atom (Fig. 2b). Furthermore, the phonon dispersion of bulk graphite and h-BN (Fig. 2c, d) calculated using the NEP model aligns well with experimental measurements47,48, despite minor deviations along certain high-symmetry paths. A detailed comparison of the sliding PESs and phonon dispersions for other systems against DFT references, provided in Figs. S2 and S3, reveals minimal discrepancies, underscoring the ability of the NEP model to capture subtle differences in energy and force.

Fig. 2: Various benchmarks for validating the NEP model.
Fig. 2: Various benchmarks for validating the NEP model.The alternative text for this image may have been generated using AI.
Full size image

a Comparison of binding energy curves for different stacked Gr/h-BN systems, calculated using DFT and NEP model. Atomic snapshots corresponding to different stacking types are visualized. b Comparison of sliding PES predicted by DFT and NEP model for the Gr/h-BN heterostructure. The AA, AB, DW, and ABp stacking modes are marked at the corresponding positions on the PES map. c, d Comparison of the phonon dispersion between experimental measurements and NEP model predictions, with experimental data for bulk graphite and h-BN extracted from refs. 47 and 48, respectively.

Utilizing the equilibrium lattice constants obtained from the NEP model, which are very close to those given by DFT calculations and experiments (see Supplemental Table S1 for the equilibrium lattice constants and interlayer spacing), we constructed a bilayer Gr/h-BN heterostructure unit cell with minimal lattice mismatch and in-plane dimensions of 13.8 × 23.9 nm (see Methods section for details) to generate moiré patterns12,13,21,49. The inherent lattice mismatch in the heterostructure causes interlayer atoms to form stacking sequences that are less energetically favorable. To accommodate the lattice disparity, during geometric optimization, atoms in each layer shift to minimize the system’s energy, resulting in reconstructed surface regions where flat, near-optimal stacking areas are separated by narrow, elevated DWs36,50.

Figure 3a, b illustrate the unit-cell moiré pattern and its corresponding supercell of the Gr/h-BN heterostructure optimized using the NEP model. Different regions display various stacking configurations that align well with the moiré patterns observed experimentally with a scanning tunneling microscope (STM). The atomic reconstruction of the heterostructure leads to a loss of planarity, forming a strongly corrugated surface (Fig. 3c). The moiré pattern periodicity from the NEP model (about 13.7 nm, see Fig. 3d) is well consistent with the results of the theoretical model at a zero-twist angle12 and atomic force microscopy (AFM) measurements21, further validating that the NEP model accurately captures the sliding energy landscape corrugation in the Gr/h-BN heterostructure (see Fig. 2a, b). Additionally, the empirical potentials Tersoff+ILP and Tersoff+LJ can reproduce the experimentally observed moiré patterns and their periodicities (refer to Fig. S4). The LJ term (LJ parameters can be found in Ref. 51), however, significantly underestimates the overall sliding energy landscape corrugation in Gr and h-BN25, leading to substantial suppression of the height fluctuations in the moiré pattern obtained with the Tersoff+LJ potential (see Fig. S4b).

Fig. 3: Moiré patterns and periodicities of the Gr/h-BN heterostructure.
Fig. 3: Moiré patterns and periodicities of the Gr/h-BN heterostructure.The alternative text for this image may have been generated using AI.
Full size image

a The unit-cell moiré pattern obtained from geometric optimization using the NEP model, containing AA, AB, DW, and ABp stacking types. b The supercell of the moiré patterns. The color bar represents the z-coordinate of the atoms, showing the undulating ripples. The upper right corner shows the moiré patterns observed experimentally using the STM13. c Schematic diagram of the free-standing bilayer Gr/h-BN heterostructure. d Comparison of the periodicity (wavelength, λ) of moiré patterns obtained through different methods. The modeling result12 corresponds to the Gr/h-BN heterostructure with a zero-degree twist angle, while the experimental data are derived from AFM measurements21.

Further, the NEP model is employed to perform geometric optimization of a bilayer h-BN with a 0.385∘ twist angle52, successfully replicating the shape of the experimentally observed moiré pattern49 and highlighting the stacking types in different regions (see Fig. S5). This achievement demonstrates the strong generalization capability of the NEP model, even without prior training on twisted bilayer h-BN configurations.

For MD simulations, our NEP model, capable of capturing long-range dispersion interactions, achieves an impressive computational speed of 9.4 × 106 atom-step per second on four Nvidia GeForce RTX 4090 GPU cards (see Fig. S6), closely rivaling the performance of the Tersoff+LJ empirical potential, which attains 1.3 × 107 atom-step per second using 72 Intel Xeon Gold 6420 CPU cores. In contrast, the Tersoff+ILP potential, running on 108 Intel Xeon Gold 6420 CPU cores, achieves only 1.7 × 105 atom-step per second, as ILP requires frequent updates of neighboring atoms and employs a large radial cutoff. The fast and accurate NEP model provides foundations for a comprehensive study of thermal transport in Gr/h-BN heterostructure.

Ideal limit of interfacial thermal conductance in 2D heterostructures

By employing the NEMD method (refer to the Methods section for details) on a realistically large system containing over 300,000 atoms, we calculated the ITC of the Gr/h-BN heterostructure, as presented in Table 1 and Fig. 4a. At around 300 K, the ITC values calculated from the Tersoff+LJ and Tersoff+ILP potentials deviate significantly and moderately, respectively, from photothermal Raman and TDTR experimental measurements, whereas the NEP model yields results that closely match the experimental ones. The distribution of ITC data for different potentials from 20 independent calculations is shown in Fig. S9. Classical MD simulations, however, cannot inherently account for quantum Bose-Einstein statistics38,39,40. To address this limitation, we performed a proper quantum-statistical correction to the classical ITC based on spectral thermal conductance (see Methods section), achieving a closer match with experimental values at 300 K (see Table 1).

Table 1 Comparison of the ITC of the Gr/h-BN heterostructure at 300 K, as calculated using NEMD method driven by various potentials, with experimental results
Full size table
Fig. 4: ITC data for various heterostructure systems.
Fig. 4: ITC data for various heterostructure systems.The alternative text for this image may have been generated using AI.
Full size image

a The correlation between ITC and PDOS overlap at room temperature for different 2D heterostructures. The light orange band emphasizes the ideal limit range of ITC for 2D heterostructures. The values of PDOS overlap (see Fig. S12 for details) are marked below the labels of different systems, decreasing from left to right. The ITC data measured by Raman and TDTR experiments are highlighted with light blue shading, including Gr/h-BN3,6,7,15, Gr/MoS28,9, h-BN/MoS26, MoS2/WSe28, h-BN/WSe211, h-BN/WS211, and Gr/WSe28 heterostructures. The ITC for the Gr/h-BN heterostructure obtained from NEMD19,20, and AEMD17 simulations are further compared with experimental ones. For clarity, the essential data are presented in Table 1. b Temperature-dependent ITC for different heterostructure systems. For the Gr/h-BN heterostructure, the ITC includes the classical and quantum-corrected results from NEP calculations (with no available experimental data). For 2D/SiO2 heterostructures, including Gr/SiO272, h-BN/SiO273, MoS2/SiO274,75, and MoSe2/SiO275, the ITC values are based on experimental measurements.

We also used the approach-to-equilibrium MD (AEMD) technique (see Methods section for details) to calculate the ITC of Gr/h-BN heterostructure, revealing a significant underestimation compared to the experimental ones (see Supplemental Note S4 and Fig. S8), likely because some low-frequency phonons with long lifetimes do not contribute to thermal conductance in AEMD53. Additionally, we confirmed that the Gr/h-BN heterostructure system exhibits almost no thermal rectification (Fig. S10).

Among the various 2D heterostructures formed by Gr, h-BN, MoS2, WS2, and WSe2, the Gr/h-BN heterointerface achieves the highest ITC (Fig. 4a), because Gr and h-BN have the highest similarity in structure and mass14. Conversely, the Gr/WSe2 interface exhibits the lowest ITC herein, reflecting a pronounced disparity between the constituent materials, characterized by a substantial unit-cell mass mismatch (as shown in Fig. S11) and differences in thermal properties. The relative “similarity” between a pair of 2D materials can be simplified to two key factors: the closeness of their average atomic masses (Fig. S11), determining the energy transmission probability, and the degree of overlap in their PDOSs (Fig. S12), indicating the phonon modes available for energy transmission5,54. Figure 4a and Supplemental Table S2 highlight the correlation between PDOS overlap and ITC across different 2D heterostructures, revealing a clear trend that greater overlap corresponds to higher ITC. To provide a quantitative metric, we perform a Spearman correlation analysis between the PDOS overlap and the averaged ITC values (Fig. S13), which yields a Spearman coefficient of ρ = 0.78 with a p-value of 5.76 × 10−4, indicating a statistically significant and strong positive monotonic relationship. These results suggest that the PDOS overlap can serve as a straightforward positive predictor for determining ITC.

In general, the intrinsic ITC of 2D heterostructures through van der Waals interactions is expected to have an upper limit, which we refer to here as the ideal limit of ITC. Approaching this ideal limit experimentally is fraught with challenges. The fabrication of a flawless interface is inherently difficult due to the susceptibility of 2D materials to impurity contamination and wrinkling during exfoliation and transfer. Moreover, the strong dependence of experimental outcomes on the “cleanliness" of the interface, coupled with the low measurement sensitivity of ITC3,6,7,15, further complicates the accurate probing of this upper bound. Using the Gr/h-BN system, which demonstrates the highest ITC among 2D heterostructures (see Fig. 4a and Supplemental Table S2), as a prototype, and employing NEP-driven NEMD simulations augmented with the PDOS overlap criterion, we determined the ideal limit of ITC for 2D heterostructures at room temperature to be 58.4 ± 11.0 MW m−2 K−1. In Fig. 4a, we highlight this ideal limit range, anticipating it to serve as a reference standard.

Figure 4b compares the temperature-dependent ITC values for the Gr/h-BN heterostructure, as predicted by the NEP model, with experimental data for various 2D/SiO2 heterostructures. Fundamentally, the temperature dependence of ITC arises from the temperature dependence of the heat capacity. Below the Debye temperature (ΘD), the population of phonon modes available for heat transfer across the interface increases with temperature. Similar to heat capacity, however, the ITC saturates and remains constant for T > ΘD. The “softer” material (with weaker bond strength and larger atomic mass) at the interface will ultimately limit the ITC. At temperatures above ΘD,soft, the high-energy phonon modes of the “harder” material remain populated but fail to transmit efficiently across the interface due to the absence of matching phonon modes in the “softer” material. This mechanism likely accounts for the weak temperature dependence of ITC in MoS2/SiO2 and MoSe2/SiO2, as the Debye temperatures of transition metal dichalcogenides are significantly lower than those of graphite and h-BN (see Fig. S14).

As mentioned earlier, classical MD cannot account for quantum statistics, resulting in the ITC of the Gr/h-BN heterostructure, calculated using the NEP model (depicted as circles in Fig. 4b), exhibiting an almost temperature-independent behavior. By incorporating quantum-statistical corrections based on spectral thermal conductance (see Methods section and Fig. S15), we reveal a moderate temperature dependence of the ITC (denoted by pentagrams), comparable to that of the h-BN/SiO2 system. The temperature dependence of ITC for a broader range of heterostructures is comprehensively compared in Fig. S16.

Stacking-sequence-dependent interfacial thermal transport

It is now timely to delve into the quantum-statistical corrections based on spectral thermal conductance. As shown in Fig. 5a, the quantum-corrected spectral thermal conductance is slightly smaller than the classical one at 300 K, with substantial quantum-correction effects becoming apparent only at low temperatures (50 K, see Fig. S15). This aligns with the moderate temperature dependence of the ITC demonstrated in Fig. 4b. The ITC of the Gr/h-BN heterostructure is predominantly driven by contributions from long-wavelength, low-frequency (ω/2π < 10 THz) phonons. This finding is further corroborated by the spectral thermal conductance calculated using the empirical potentials Tersoff+ILP and Tersoff+LJ (Fig. S17), which consistently highlight the dominant role of low-frequency phonons.

Fig. 5: Spectral thermal conductance, interlayer energy and distance of Gr/h-BN heterostructure.
Fig. 5: Spectral thermal conductance, interlayer energy and distance of Gr/h-BN heterostructure.The alternative text for this image may have been generated using AI.
Full size image

a Classical and quantum-corrected spectral thermal conductance of Gr/h-BN heterostructure at 300 K calculated using NEP. The inset shows a magnified view covering the 1–5 THz to provide a clearer comparison between the classical spectral thermal conductance and the quantum-corrected one. The gray dashed line emphasizes the central frequency of the peak for the spectral thermal conductance. b, c Quantum-corrected spectral thermal conductance and the corresponding accumulated thermal conductance for different stacking types. d, e Interlayer van der Waals potential energy (characterizing long-range dispersion interactions) and interlayer distance of the Gr/h-BN heterostructure. The distinct stacking types corresponding to different regions are highlighted, with the gray dashed lines representing the unit cell.

Furthermore, we evaluate the quantum-corrected spectral thermal conductance and their corresponding cumulative conductance for various stacking sequences at the Gr/h-BN heterointerface (see Methods section), revealing a pronounced dependence of ITC on stacking sequences, as presented in Fig. 5b, c. The interlayer coupling strength and distance can act as straightforward metrics to interpret this phenomenon. The van der Waals potential energy and interlayer distance30 obtained after geometric optimization (Fig. 5d, e) reveal that stacking sequences with lower/higher ground-state energies correspond to smaller/larger interlayer distances, thereby allowing more/less out-of-plane phonon energy to transmit through the heterointerface. For different stacking sequences, the interlayer potential energy and distance are ranked as ABp > DW > AB > AA, which corresponds to the ITC ranking of ABp > DW > AB  > AA.

To further elucidate the stacking-sequence-dependent ITC behavior, it is instructive to examine the vibrational spectra along the Γ–A path36,55,56. Using the spectral energy density (SED) technique, we calculated the energy-weighted phonon dispersion and derived the phonon lifetimes for various stacking configurations (see Methods section). In Fig. 6a–c, a pronounced collapse (softening) of the transverse acoustic (TA) and longitudinal acoustic (LA) phonon modes, which dominate interfacial thermal transport36,55,56, is observed in the AA and AB stacking types compared to the ABp one. Additionally, the ABp stacking type demonstrates larger group velocities (Fig. S18) and longer phonon lifetimes (Fig. 6d) than the other two. Figure S19 presents the phonon lifetimes derived from the Lorentzian fitting of the SED, vividly showcasing the softening (redshift) of the TA and LA phonon modes when transitioning from ABp to AA stacking. These observations collectively indicate that energy-favorable stacking facilitates out-of-plane phonon energy transmission through reducing interlayer distances, manifesting in higher group velocities and prolonged phonon lifetimes, thereby driving the stacking-sequence-dependent ITC.

Fig. 6: Phonon SED and the corresponding lifetimes.
Fig. 6: Phonon SED and the corresponding lifetimes.The alternative text for this image may have been generated using AI.
Full size image

a–c Phonon SED (from Γ to A, ω/2π < 5 THz) for different stacking types of the Gr/h-BN heterostructure. The grey solid line represents lattice dynamics calculations at 0 K. d Phonon lifetimes for different stacking configurations, obtained through Lorentzian fitting of the SED peaks (see Fig. S19).

We averaged the spectral thermal conductance calculated over the entire heterointerface and across all distinct stacking configurations, achieving quantitative agreement (Fig. S20). In light of this, to obtain ITC values for 2D heterostructures from the NEMD+NEP method that are comparable to experimental measurements, it is insufficient to consider only a single stacking type in the computational model. Instead, at least one moiré period (see Fig. 3) containing multiple stacking types should be included. During the NEMD process at 300 K, the moiré pattern at the heterointerface remains clearly visible (see Fig. S21).

The stacking-sequence-dependent ITC can be further correlated with moiré patterns56. In 2D layered materials, the essence of moiré patterns stems from lattice mismatch caused by heterogeneous stacking, twisted angles, interlayer sliding, and heterostrain engineering. To minimize system energy, atomic reconstruction redistributes interlayer stacking sequences and induces diverse moiré patterns. It can be anticipated that different moiré patterns exhibit varying stacking proportions, leading to moiré-pattern-dependent ITC governed by distinct stacking sequences. This scenario invites a more comprehensive exploration of diverse layered materials.

Discussions

In summary, given the minimal lattice mismatch between Gr and h-BN, we use the Gr/h-BN system as a prototype to probe the ideal limit of ITC in 2D heterostructures. Building on the PBE functional combined with a state-of-the-art MBD correction, we develop a highly efficient NEP model for the Gr/h-BN heterostructure with near-quantum-mechanical accuracy. With this MLP model, we successfully reproduce the experimentally observed moiré patterns and periods, and perform NEMD simulations on a large Gr/h-BN heterostructure system with over 300,000 atoms. Through a proper quantum-statistical correction to the classical MD results, we achieve ITC values quantitatively consistent with experimental ones, thereby determining the ideal limit of ITC to be 58.4 ± 11.0 MW m−2 K−1 for 2D heterostructures at room temperature. This ideal limit is expected to serve as a standard reference for thermal design in 2D-based nanoelectronics.

Beyond this breakthrough, we unveil a novel stacking-sequence-dependent ITC hierarchy: ABp  > DW  > AB  > AA in the Gr/h-BN heterostructure. The mechanism behind this picture can be attributed to stacking sequences with larger interlayer distances being unfavorable for out-of-plane phonon energy transmission, as evidenced by lower group velocities (the collapse of TA and LA phonon modes from ABp to AA stacking) and reduced phonon lifetimes, which ultimately drive the stacking-sequence-dependent ITC. This phenomenon can be further associated with moiré patterns, potentially giving rise to moiré-pattern-dependent ITC. Inspired by these observations, we consider that this stacking-sequence-dependent thermal transport scenario is likely universal in van der Waals layered materials and may extend to other properties such as ferroelectricity, spintronics, and magnetism57,58,59.

Methods

Construction of training dataset

Constructing the training structures of layered materials like Gr and h-BN for the NEP model is relatively straightforward, with sampling primarily from structural perturbations to capture diverse stress conditions and MD simulations to account for various temperature regimes37. Starting from the optimized structure, we introduced random cell deformations (−3% to 3%) and atomic displacements (within 0.1 Å) for structural perturbations. For MD simulation, DFT-driven MD simulation (see Supplemental Note S1 for details) is the primary approach, while empirical potential-driven MD simulation serves as the auxiliary one.

The training dataset comprises four types of structures: monolayers, bilayers, trilayers, and bulk materials (see Fig. 1a for typical structures). For monolayers, we generated 50 perturbed structures each for Gr and h-BN and uniformly sampled 200 additional structures from each DFT-MD simulation, yielding a total of 500 monolayer structures. Using the same methods, we generated 500 bulk structures for graphite and h-BN. To accurately describe the binding energy curves, we varied the interlayer spacing to produce 59 additional bulk structures for graphite and h-BN, resulting in a total of 618 bulk configurations.

Due to the large number of atoms in trilayer structures (see Fig. 1a), we performed empirical potential-driven (Tersoff+ILP)23,27,28,29,30 MD simulations in the isothermal-isobaric (NpT) ensemble, generating 787 structures by linearly increasing the target temperature from 10 K to 1500 K over 200 ps.

For bilayer graphene, we considered AA and AB stackings (see Fig. S1c) and generated 200 structures through perturbations and varying interlayer distances. In the case of bilayer h-BN, we explored five stacking types-AA, AAp, AB, ABp, and ApB (see Fig. S1d)-producing 500 structures. Regarding the Gr/h-BN heterostructure, we considered AA, AB, and ABp stackings (see Fig. 2a), producing a total of 900 structures using perturbations, DFT-driven MD, and varied interlayer spacings.

Additionally, rigid sliding configurations of bilayer Gr (AB stacking), bilayer h-BN (AA stacking), and bilayer Gr/h-BN heterostructure (AA stacking) were incorporated into the training set. Starting from the optimized configurations, we constructed a set of laterally shifted structures by rigidly translating the top layer60. In each case, a 6 × 11 grid of shifted configurations was generated along the zigzag and armchair directions, with each frame including two perturbed structures, yielding 198 structures per case and 594 configurations overall.

In total, we generated a comprehensive dataset of 4099 structures (containing 229,760 atoms), all of which are utilized as the training set for the NEP model.

DFT calculations for reference data generation

After preparing the initial training configurations, the Vienna Ab initio Simulation Package (VASP) that implemented the projected augmented wave (PAW) method61 and a plane-wave basis set62 was used to perform quantum-mechanical calculations to obtain reference data, including the energy and virial of each structure, as well as the forces on each atom in the structure. The electronic exchange-correlation interaction was described by the PBE functional44 within the generalized gradient approximation. Single-point calculations were performed with an energy cutoff of 850 eV and a Γ-centered k-point mesh, using a spacing of 0.15/Å to sample the Brillouin zone. The electronic self-consistent loop was converged to a threshold of 10−8 eV. For accurate total energy calculations, the tetrahedron method with Blöchl corrections was employed. To prevent periodic interactions in the z-direction, a vacuum layer of 50 Å thickness was applied. The state-of-the-art MBD correction45,46 was employed to accurately capture the long-range dispersion interactions. Details of the INCAR file used for generating the reference data are provided in Supplemental Note S2.

The NEP model training

We trained the NEP model for Gr/h-BN heterostructure by employing the fourth-generation (NEP4) scheme34. The NEP approach32, based on a feedforward neural network (NN), uses a single hidden layer with Nneu neurons to represent the site energy Ui of atom i as a function of a descriptor vector with Ndes components,

$${U}_{i}=\mathop{\sum }\limits_{\mu =1}^{{N}_{{\rm{neu}}}}{w}_{\mu }^{(1)}\tanh \left(\mathop{\sum }\limits_{\nu =1}^{{N}_{{\rm{des}}}}{w}_{\mu \nu }^{(0)}{q}_{\nu }^{i}-{b}_{\mu }^{(0)}\right)-{b}^{(1)},$$
(1)

where \(\tanh (x)\) is the activation function, w(0), w(1), b(0), and b(1) are the trainable weight and bias parameters in the NN. The local atom-environment descriptor \({q}_{\nu }^{i}\) is constructed from a set of radial and angular components33.

Optimizing the NEP model parameters involves minimizing a loss function that combines a weighted sum of RMSEs for energy, force, and virial, along with regularization terms. After extensive benchmarking, we identified the optimal hyperparameters for the NEP model training, available in the nep.in file (see Supplemental Note S3). To more accurately describe the binding energy curves and sliding PESs, we increased the weight of the 816 configurations for the four atoms in the training set from 1 to 8, highlighting their relative importance in the total loss function.

NEMD simulations

We primarily employed the NEMD method to calculate the ITC of the Gr/h-BN heterostructure. Given the inherent in-plane lattice mismatch of approximately 1.8% between Gr and h-BN12,13,21,22, we constructed supercells by expanding the Gr unit cell 56-fold and the h-BN unit cell 55-fold along both the armchair and zigzag directions. This resulted in in-plane dimensions of 13.8 × 23.9 nm, significantly reducing the in-plane size effects on ITC. Subsequently, we constructed a multilayer Gr/h-BN heterostructure model featuring a single heterointerface, comprising 12 layers of graphene (12,544 atoms per layer) and 13 layers of h-BN (12,100 atoms per layer), totaling 307,828 atoms. Through testing, this layer configuration effectively minimizes the impact of out-of-plane dimensions on ITC while maintaining computational efficiency. The initial interlayer spacing was set to 0.35 nm, and subsequently optimized during equilibration to a more physically realistic value. Detailed information on the model can also be found in Fig. S7.

The initial structure was first equilibrated for 0.5 ns under an NpT ensemble (with a zero target pressure) at the target temperatures. This was followed by a 0.25 ns run in the NVT ensemble. We adopted a time step of 0.5 fs, which has been proven sufficiently small for accurate results. During the production stage lasting 2.5 ns, two Langevin local thermostats63, specifically a heat source and sink (with a temperature difference of 60 K), are employed to generate a nonequilibrium steady state with constant heat flux. At this stage, a noticeable temperature jump ΔT, occurs at the Gr/h-BN heterointerface (see Fig. S7c), representing the average interfacial temperature difference between Gr and the neighboring h-BN layer. One can calculate the ITC (G) according to Fourier’s law:

$$G=\frac{dE/dt}{S\Delta T},$$
(2)

where S is the cross-sectional area in the transverse directions and dE/dt is the average energy exchange rate between the thermostats and the thermostatic regions.

AEMD simulations

The AEMD technique64,65 is used to cross-validate the results from the NEMD method, and the ITC value reported in previous work17. After equilibrating the constructed Gr/h-BN bilayer model (see Fig. S8), we switched to an NVE ensemble and applied two Langevin thermostats: one set to 350 K for Gr and another set to 250 K for h-BN. This step was performed for 1 ns to establish a stable temperature difference of ΔT0 = 100 K. After removing the thermostats, we monitored the temperature difference ΔT(t) over 2 ns, which decayed exponentially over time t according to \(\Delta T(t)=\Delta {T}_{0}\exp (-t/\tau )\), where τ is the decay time to be determined (see Fig. S8c). Based on the lumped heat-capacity model64, the ITC can be extracted as

$$G=\frac{{C}_{V}}{S\tau },$$
(3)

where S is the cross-sectional area. Here, and CV is the effective constant-volume heat capacity of the Gr/h-BN bilayer system, calculated as \({C}_{V}=\frac{{C}_{{\rm{Gr}}}\cdot {C}_{{h}\text{-}BN}}{{C}_{{\rm{Gr}}}+{C}_{{h}\text{-}BN}}\), where both CGr and Ch-BN are given by 3NkB, with kB being the Boltzmann constant.

In the NEMD and AEMD simulations, periodic boundary conditions were applied in the in-plane directions, while non-periodic boundary conditions were used along the heat transport direction. The MD simulations with the NEP model were performed using the GPUMD package66 (version 3.9.2), while those using empirical potentials were carried out with LAMMPS67 package (version 2 Aug 2023).

Spectral thermal conductance and quantum-statistical corrections

Within the NEMD framework, one can first calculate the virial-velocity correlation function in the nonequilibrium steady state68,

$${\boldsymbol{K}}(t)=\sum _{i}\left\langle {{\bf{W}}}_{i}(0)\cdot {{\boldsymbol{v}}}_{i}(t)\right\rangle .$$
(4)

where vi and Wi are the velocity and virial tensor of atom i, respectively. Then, the classical thermal conductance can be spectrally decomposed:

$$G(\omega ,T)=\frac{2}{V\Delta T}\mathop{\int}\nolimits_{-\infty }^{\infty }{\rm{d}}t{{\rm{e}}}^{i\omega t}K(t),$$
(5)

where V is defined as the product of the cross-sectional area and the interlayer distance across the Gr/h-BN interface. With the classical spectral G(ω, T) available, a quantum-corrected spectral thermal conductance Gq(ω, T) can be obtained by multiplying it with a ratio between quantum and classical modal heat capacity38,39,40,

$${G}^{{\rm{q}}}(\omega ,T)=G(\omega ,T)\frac{{x}^{2}{e}^{x}}{{({e}^{x}-1)}^{2}},$$
(6)

where x = ℏω/kBT, ℏ is the reduced Planck constant, and kB is the Boltzmann constant. Then, the quantum-corrected thermal conductance is then obtained as an integral over the entire frequency range as

$${G}^{{\rm{q}}}(T)=\mathop{\int}\nolimits_{0}^{\infty }\frac{{\rm{d}}\omega }{2\pi }{G}^{{\rm{q}}}(\omega ,T).$$
(7)

For the quantum-corrected spectral thermal conductance of different stacking types in Gr/h-BN heterostructure, to prevent displacement of the stacking regions, a few atoms at the edges of the interfacial layers were fixed after optimizing the multilayer model. The regions corresponding to different stacking types in the two-layer heterointerface (see Fig. 3a and Fig. S7) were then selected for calculation. In this case, the volume V is calculated as the product of the area of the selected region and the average interlayer distance30 between the two atomic layers. The reported spectral thermal conductance results for different stacking types (Fig. 5) are the averages of ten independent simulations.

Phonon spectral energy density calculations

The phonon SED technique69 is employed to examine the phonon dispersion and corresponding phonon lifetimes for the AA, AB, and ABp stacking configurations in the Gr/h-BN heterostructure. It represents the intensity of the kinetic energy distribution of a system over different phonon modes and can be calculated in MD simulations using the following formula69:

$$\begin{array}{lll}\Phi ({\bf{q}},\omega )=\displaystyle\frac{1}{4\pi {\tau }_{0}{N}_{T}}\mathop{\sum }\limits_{\alpha }^{3}\mathop{\sum}\limits_{b}^{n}{m}_{b}\left|\mathop{\displaystyle\int}\nolimits_{0}^{{\tau}_{0}}\mathop{\sum}\limits_{l}^{{N}_{T}}{\dot{u}}_{\alpha }(l,b,t)\right.\\\qquad\quad\quad\,\,\,\,\,\times\exp{\left.(i{\bf{q}}\cdot{{\bf{r}}}_{0}(l)-i\omega t)dt\vphantom{\mathop{\displaystyle\int}\nolimits_{0}^{{\tau}_{0}}}\right|}^{2}.\end{array}$$
(8)

Here, τ0 is the total MD simulation time, while n and NT represent the total basis atoms and the total unit cells in the system, respectively. The mass of the b-th basis atom is denoted by mb, while \({\dot{u}}_{\alpha }(l,b,t)\) represents the α-th Cartesian direction of the velocity of the b-th basis atom in the l-th unit cell at time t, and r0(l) indicates the equilibrium position of the l-th unit cell.

The phonon lifetime can be obtained by fitting the SED curve with the Lorentzian function69,

$$\Phi ({\bf{q}},\omega )=\frac{I}{1+{\left[(\omega -{\omega }_{c})/\gamma \right]}^{2}},$$
(9)

where I is the peak magnitude, ωc is the frequency at the peak center, and γ is the half-width at half-maximum. Finally, one can define the lifetime τph at each wave vector q and frequency ω as:

$${\tau }_{ph}({\bf{q}},\omega )=\frac{1}{2\gamma }.$$
(10)

Starting from the equilibrium configurations of the primitive cells for different stacking types (with interlayer distances of 3.56 Å for AA, 3.49 Å for AB, and 3.37 Å for ABp), a 6 × 6 × 50 supercell was constructed for each configuration. All systems were first equilibrated at 300 K for 0.5 ns in the NVT ensemble, followed by a 3 ns production stage in the NVE ensemble. During this period, the velocities and positions of all atoms were collected every 50 fs, for use in the subsequent SED calculations. As the DW stacking represents a transitional state and cannot exist independently, it was excluded from the SED calculations. The SED calculation and Lorentzian fitting were performed using the PYSED package70, while the lattice dynamics calculation was carried out with the PHONOPY package71.