Predicting the thermal conductivity of two-dimensional (2D) heterostructures is challenging and cannot be adequately resolved using conventional computational approaches. To address this challenge, we propose a new and efficient approach that combines first-principles density functional theory (DFT) calculations with a machine-learning interatomic potential (MLIP) methodology to determine the thermal conductivity of a novel 2D van der Waals TiS2/MoS2 heterostructure. We leverage the proposed approach to estimate the thermal conductivities of TiS2/MoS2 heterostructures as well as bilayer-TiS2 and bilayer-MoS2. A unique aspect of this approach is the combined implementation of the moment tensor potential for short-range (intralayer) interactions and the D3-dispersion correction scheme for long-range (interlayer) van der Waals interactions. This approach employs relatively inexpensive computational DFT-based datasets generated from ab initio molecular dynamics simulations to accurately describe the interatomic interactions in the bilayers. The thermal conductivities of the bilayers exhibit the following trend: bilayer-TiS2 > bilayer-MoS2 > the TiS2/MoS2 heterostructure. In addition, this work makes the case that the 2D bilayers exhibit considerably higher thermal conductivities than bulk graphite, a common battery anode material, indicating the potential to utilize 2D heterostructures in thermal management applications and energy storage devices. Furthermore, the MLIP-based methodology provides a reliable approach for estimating the thermal conductivity of bilayers and heterostructures.

The increase in power density in modern electronics has led to a growing demand for efficient heat removal in electronic devices, including communication, information, and energy storage technologies, where a higher cooling rate enables faster operation and improved device performance.1 The excessive heat produced during operating cycles can affect these electronic and energy storage devices, such as batteries, resulting in heat concentration in localized regions called hotspots.1 Devices employing materials with low thermal conductivity can experience deteriorated performance or failure due to insufficient mitigation of increased heat concentration in the hotspots. Thermal management encompasses the implementation of temperature monitoring and device cooling to regulate the operational temperature of devices. Thermal management techniques for eliminating hotspots involve using materials with high thermal conductivity to ensure the safe and efficient operation of electronic and energy storage devices. Two-dimensional (2D) materials offer promising prospects for thermal management.2,3

Two-dimensional (2D) materials exhibit distinct thermal properties compared to their bulk counterparts. These include enhanced thermal conductivity due to the ballistic4 and hydrodynamic5 transport of phonons.6 Hence, 2D materials, such as graphene, boron nitride, and transition-metal dichalcogenides (TMDs), have demonstrated excellent thermal performance.1,7 Similarly, several other 2D materials, including black-phosphorene,8 C3N,9 and MoS2,10 have exhibited superior thermal properties than graphite, a commonly employed anode material in batteries. Building heterostructures using 2D materials offers a practical solution to address the limitations of monolayers, while still retaining some of the advantages of pristine 2D materials. 2D heterostructures also demonstrate exceptional properties in comparison with their individual monolayers, including improved electrochemical activity, structural stability, superior mechanical strength, and enhanced ionic conductivity.11 In addition, introducing lattice mismatches at the interfaces in heterostructures leads to variations in electrical current and thermal transport that are distinct from monolayers.

Estimation and optimization of thermal conductivity in 2D materials and their heterostructures offer significant potential in thermal management applications. Improving thermal conductivity in 2D materials enables effective heat dissipation and reduces energy losses. The thermal conductivity of 2D materials and heterostructures can be predicted by solving the Boltzmann transport equation (BTE) using density functional theory (DFT),12–15 which is a highly accurate method but also highly computationally intensive, particularly for lattice structures with large unit cells. Classic molecular dynamics (MD) simulations are more cost-effective in predicting the thermal conductivity of nanostructured materials with complex and large lattice structures.16 However, the accuracy of MD simulations depends significantly on the choice of interatomic potentials, which in turn restricts its widespread use because accurate potentials are not readily available for novel 2D materials and heterostructures. For example, experimental measurements of thermal conductivities of graphene monolayers reported in the literature range from 1500 to 5300 W/mK, while thermal conductivities predicted from MD simulations using different interatomic potentials range from 350 to 3000 W/mK.17 These findings demonstrate significant discrepancies among experimental and MD-predicted thermal conductivities of up to an order of magnitude. Hence, it is essential to develop efficient interatomic potentials for 2D materials and their heterostructures. Machine learning (ML) has demonstrated its effectiveness in creating accurate interatomic potentials called machine learning interatomic potentials (MLIPs) for MD simulations.18 MLIPs employ a regression algorithm to compute the ab initio potential energy surface (PES). This approach utilizes atomic configurations as input features, enabling the determination of interatomic force constants with minimal computational cost. This allows for a more rapid solution of the phonon Boltzmann transport equation (BTE) to evaluate the thermal conductivity. Previous studies have extensively investigated thermal conductivity in 2D monolayers using machine-learning-derived interatomic potentials.17,19 These potentials have demonstrated remarkable accuracy compared to experimental data and other theoretical works. For the estimation of thermal conductivity in 2D bilayers and heterostructures, various factors, including different stacking possibilities, grain boundaries, moiré patterns, computational costs, and lattice mismatches, have hindered the development of an effective interatomic potential.

Two-dimensional TiS2/MoS2 van der Waals heterostructures have shown promise as potential anode materials in alkali-ion batteries due to their excellent electrochemical characteristics, including strong alkali-ion adsorption and low diffusion energy barriers. Furthermore, these heterostructures also exhibit good thermal stability at high temperatures.20 However, estimating the thermal conductivity of TiS2/MoS2 heterostructures using DFT alone is computationally prohibited. In addition, the lack of an interatomic potential for novel materials makes it challenging to perform classical MD simulations. This work presents the development of an interatomic potential using machine learning to determine the thermal conductivity of a TiS2/MoS2 heterostructure. We also predicted the thermal conductivity of individual MoS2 and TiS2 bilayers. Many machine learning methodologies require a complex implementation involving numerous single-point DFT calculations during the training process to derive the interatomic potential.17 Here, we develop highly accurate interatomic potentials with machine learning techniques using cost-effective ab initio molecular dynamics (AIMD) trajectories over less than 1 ps, eliminating the need for additional DFT calculations. Our study aims to explore the potential of machine learning in accurately determining the thermal properties in 2D heterostructures and bilayers. Furthermore, this study can greatly contribute to investigating advanced thermally conductive materials for thermal management applications.

In summary, this work’s most relevant contribution is the implementation of the moment tensor potential in conjunction with the D3-dispersion correction scheme to approximate the particular interlayer interactions of 2D bilayers and heterostructures. This approach is more effective than the typical Lennard–Jones (L-J) potential for long-range (interlayer) van der Waals (vdW) interactions. This work develops highly accurate interatomic potentials with machine learning techniques using cost-effective ab initio molecular dynamics (AIMD) trajectories over less than 1 ps. With this machine-learning interatomic potential (MLIP) methodology, we investigate the in-plane thermal conductivity of a novel two-dimensional TiS2/MoS2 van der Waals heterostructure.

The computational simulations in this study are conducted in the following sequential steps: We initially constructed the TiS2/MoS2 heterostructure and the bilayers of TiS2 and MoS2 using the appropriate lattice configurations. We then prepared the training datasets by conducting ab initio molecular dynamics simulations (AIMD) for the structures at low (50 and 100 K) and high (300 and 600 K) temperatures. The moment tensor potential (MTP) parameters were then obtained using the machine-learning interatomic potentials (MLIPs) scheme. We then performed classical MD simulations, utilizing the generated MLIPs, to determine the thermal conductivity of bilayer-TiS2, bilayer-MoS2, and TiS2/MoS2 heterostructures.

We performed density functional theory (DFT) calculations based on first-principles as implemented in the QuantumATK software package.21 The exchange–energy correlation interactions were simulated utilizing the generalized gradient approximation (GGA)22 in conjunction with the Perdew–Burke–Ernzerhof (PBE) functional.22 Using the linear combination of atomic orbitals (LCAO) basis sets23 with a maximum energy cutoff of 150 rydbergs, the wave functions were expanded within a double-polarized basis set subject to periodic boundary conditions.

Monolayers TiS2 and MoS2 possess two distinct phases: trigonal (1T) and hexagonal (2H). These phases are distinguished by the stacking of the S–TM–S atomic layers (TM = Ti/Mo). In the 1T phase, the atomic layers are stacked in an ABC sequence, while in the 2H phase, the stacking follows an ABA sequence. Previous studies confirmed that the 1T phase is more energetically stable than the 2H phase in monolayer TiS2,24,25 whereas MoS2 is only stable in the 2H phase.20,26 Hence, for the present investigation, we constructed a bilayer of 1T-TiS2, a bilayer of 2H-MoS2, and a heterostructure of 1T-TiS2 and 2H-MoS2. Each bilayer structure contains 72 atoms. Due to the different lattice constants of MoS2 (3.169 Å) and TiS2 (3.412 Å) monolayers, a lattice mismatch of 5.12% occurs while stacking MoS2 on TiS2. Since the lattice mismatch between two constituent materials in heterostructures may result in stacking disorder or moiré superstructures, we calculated and compared the strain and adhesive energies for the TiS2/MoS2 heterostructure (Fig. S1 of the supplementary material), which are determined by stacking one monolayer on top of another. As the strain energies (energy difference between the freestanding state and the strained state) of MoS2 (0.15 eV/unit cell) and TiS2 (0.18 eV/unit cell) in the TiS2/MoS2 heterostructure are relatively small compared to the adhesive energy (energy difference between the heterostructure and corresponding strained monolayers of TiS2 and MoS2) at the interface (0.60 eV/unit cell), the lattice mismatch will not persist. Therefore, this work’s TiS2/MoS2 heterostructure will not exhibit moiré patterns as the adhesive energy will compensate for the raised strain energies. However, in specific instances, such as in the study conducted by Kou et al.,27 the strain energies of MoS2/MoTe2 heterostructures were determined to be greater than the adhesive energy. This causes the lattice mismatch between MoS2 and MoTe2 to persist, and strain will be released by forming moiré patterns with different stacking configurations. We constructed the training dataset by utilizing structurally relaxed supercells of bilayer-MoS2 [Fig. 1(a)], bilayer-TiS2 [Fig. 1(b)], and TiS2/MoS2 heterostructures [Fig. 1(c)]. We employed a Monkhorst–Pack grid28 with a 3 × 3 × 1 k-point mesh for geometry optimization. The van der Waals (vdW) interactions were approximated using the Grimme DFT-D3 correction.29 A large vacuum space of 30 Å was applied in the out-of-plane direction to prevent any interaction between the mirror structures. The bilayers were fully relaxed until an energy convergence of 10−5 eV and a force convergence of 0.01 eV Å−1 were attained.

FIG. 1.

Schematic illustration of top and side views of the atomic structures of (a) bilayer-MoS2, (b) bilayer-TiS2, and (c) TiS2/MoS2 heterostructure. (d) Bilayer configuration depicting the different atomic interactions within the intralayer and interlayer. (e) Interlayer energy vs interlayer distance for three distinct bilayers.

FIG. 1.

Schematic illustration of top and side views of the atomic structures of (a) bilayer-MoS2, (b) bilayer-TiS2, and (c) TiS2/MoS2 heterostructure. (d) Bilayer configuration depicting the different atomic interactions within the intralayer and interlayer. (e) Interlayer energy vs interlayer distance for three distinct bilayers.

Close modal

Interlayer atomic interactions arise from the coupling of atomic layers in two-dimensional (2D) bilayers and heterostructures, playing a crucial role in characterizing their diverse properties for technological applications. In 2D heterostructures, interactions can occur between atoms within the same layer (intralayer interactions) or between atoms in adjacent layers (interlayer interactions). These interactions can be classified as short-range covalent interactions or long-range van der Waals interactions, respectively. Figure 1(d) illustrates the various types of atomic interactions in a bilayer system. We determined the optimal interlayer distance between the layers of bilayer-TiS2, bilayer-MoS2, and TiS2/MoS2 heterostructures by plotting the binding energy (Eb) as a function of the interlayer distance (dInterlayer), as shown in Fig. 1(e). The optimal interlayer distances of bilayer-TiS2, bilayer-MoS2, and the TiS2/MoS2 heterostructure were found to be 2.7, 3.4, and 2.7 Å, respectively, which are consistent with prior results.20 

We performed AIMD simulations on the fully relaxed bilayer supercells at temperatures of 50, 100, 300, and 600 K. Each simulation consisted of 1000 steps, with a time step of 1 fs, using the NVE velocity Verlet algorithm.30 We also evaluated the thermal stability of the structures. Figure S2 of the supplementary material shows the variation in the kinetic energy of bilayer-MoS2 as a function of time (fs), along with the final (t = 1000 fs) structures at the end of the simulation for 50, 100, 300, and 600 K. During the 1000 fs, the kinetic energies exhibited minimal fluctuations, indicating that the MoS2 bilayers achieved equilibrium. Furthermore, the overall structural integrity of the structures is preserved in all the cases (50, 100, 300, and 600 K), with negligible deviations in the mean lattice positions of individual atoms, without any bond breaking, demonstrating the thermal stability of bilayer-MoS2 at these temperatures. Similar observations were made regarding the thermal stability of bilayer-TiS2 and TiS2/MoS2 heterostructures at these temperatures. In addition, the vibrational frequency of phonons can provide valuable insights into the dynamic stability of the bilayer structures.31,32 Lower-temperature simulations are useful for capturing low-frequency atomic displacements, while higher-temperature simulations are advantageous for describing high-frequency displacements. A training dataset was generated by randomly selecting 80% of the 1000 steps, where each step represents an MD trajectory. The remaining 20% of the steps were utilized as a testing dataset to verify the correlation between the predicted and true values. The datasets included atomic positions (r), energy (E), forces (F), and stress (σ) as input for the fitting process implemented in the machine learning interatomic potentials (MLIP) package,33 as illustrated in Fig. 2(a). We generated three distinct interatomic potentials for bilayer-MoS2, bilayer-TiS2, and the TiS2/MoS2 heterostructure. Each interatomic potential was trained independently considering data from all four temperatures (50, 100, 300, and 600 K).

FIG. 2.

(a) Rectangular colored planes represent training data (r, E, F, and σ) generated from AIMD simulations performed at 50, 100, 300, and 600 K. Predicted (MTP) total energy values vs target (DFT) total energy values for (b) bilayer-MoS2, (c) bilayer-TiS2, and (d) the TiS2/MoS2 heterostructure, each containing 72 atoms, from the test set. (e) Schematic representation of a bilayer with neighbor atom selection for long-range interactions.

FIG. 2.

(a) Rectangular colored planes represent training data (r, E, F, and σ) generated from AIMD simulations performed at 50, 100, 300, and 600 K. Predicted (MTP) total energy values vs target (DFT) total energy values for (b) bilayer-MoS2, (c) bilayer-TiS2, and (d) the TiS2/MoS2 heterostructure, each containing 72 atoms, from the test set. (e) Schematic representation of a bilayer with neighbor atom selection for long-range interactions.

Close modal

In this work, we utilized the moment tensor potential (MTP) developed by Shapeev17,19,33,34 as a robust and computationally efficient model for describing interatomic interactions. MTPs utilize inertia tensors of different ranks multiplied by radial polynomial functions to represent atomic environments. This representation facilitates the incorporation of many-body interactions. Most machine-learning-derived interatomic potential models consider two-body (and often three-body) interactions.35 Takahashi et al. showed that incorporating explicit three-body terms can significantly enhance the accuracy of machine-learning-derived interatomic potentials for many elemental materials.36 Moment tensor potentials extend this concept by systematically incorporating two-body, three-body, and four-body interactions until the degree of the resultant polynomial becomes too large. This approach allows us to construct highly accurate potentials with sufficient training data. The functional form of the MTPs is determined by the positions and types of neighboring atoms within the cutoff radius. This functional form remains unchanged under orthogonal transformations, including rotations, reflections, and permutations of atoms of the same type.34 We employed a methodology similar to that described in previous studies17,33 to effectively train the MTPs for the generation of interatomic potentials.

The root mean squared error (RMSE) and R-squared (R2) are two relevant metrics for assessing the performance of an ML model. The RMSE metric quantifies the discrepancy between the predicted and true values in a dataset. A lower RMSE indicates a better fit between the model and the dataset. R2 determines the dispersion of the predicted data points relative to the true data. This value ranges from 0 to 1, where a higher R2 value indicates a stronger fit between the model and the dataset. We calculated the RMSE values for the energy (eV), force (eV/Å), and stress (eV/Å3) of the trained moment tensor potential (MTP) applied to bilayer-MoS2, bilayer-TiS2, and the TiS2/MoS2 heterostructure. The predicted MTPs exhibited negligible energy errors (<10−3 eV), with minimal force (<10−2 eV/Å) and stress (<10−2 eV/Å3) errors. Figures 2(b)2(d) depict the total energy plots obtained from the MTP and reference (DFT) methods for bilayer-MoS2, bilayer-TiS2, and TiS2/MoS2 heterostructures, respectively, demonstrating an excellent correlation between the predicted (MTP) and true (DFT) values. The target DFT total energy values are extracted from the AIMD simulation. In addition, we utilized 2D histograms (Fig. S3) to illustrate the error distributions of the total energies and forces for both the training and test sets of bilayer-MoS2, bilayer-TiS2, and TiS2/MoS2 heterostructures. As shown in Fig. S3, most training/test total energy and force errors for bilayer-MoS2, bilayer-TiS2, and TiS2/MoS2 heterostructures are predominantly clustered around the zero-error line. This indicates that our predicted interatomic potentials exhibit a high level of accuracy, consistent with other studies.37–40 

The MLIPs trained on computationally inexpensive DFT datasets accurately define the short-range intralayer interactions. To fully characterize the properties of bilayers and heterostructures, it is also necessary to consider the long-range interlayer interactions. In a prior study conducted by Javvaji et al.,33 the authors considered the interlayer van der Waals (vdW) interactions by employing a Lennard-Jones (L-J) potential for the bilayers. The long-range van der Waals (vdW) interlayer interaction parameters were obtained by comparing the interlayer energy measured using DFT with the energy recorded from separate single-step MD simulations. The MD simulations involved varying the interlayer distance and recording the energy for a range of guessed L-J parameters for each atom type. Their study proved the utility of this approach in investigating the bending flexoelectricity in 2D van der Waals bilayers. However, in a separate study, Ding et al. confirmed that the interatomic interactions in the L-J potential between atoms of two distinct layers have a minimal impact on the overall thermal conductivity of a bilayer.41 

In our study, we incorporated a D3-dispersion correction42 scheme into the MTP to approximate the interlayer interactions. Dispersion interactions are the important component of a van der Waals (vdW)-type interaction potential between atoms or molecules that are not directly bonded to each other. This proposed approach is more accurate than the LJ potential-based approach as it incorporates the missing dispersion contribution into the already accurate MTP. Furthermore, the D3-dispersion correction scheme has been validated by several benchmarks with accurate reference values, demonstrating its good performance in noncovalent interactions in materials.43 Another important benefit of the D3-dispersion correction scheme is its computational efficiency and robustness across multiple materials, making it particularly valuable for analyzing large systems.42 We utilized the MTP with D3-dispersion correction scheme to determine the optimal interlayer distance between the layers of bilayer-TiS2, bilayer-MoS2, and TiS2/MoS2 heterostructures, as illustrated in Fig. S4 of the supplementary material. The interlayer distances calculated using the MTP with D3-dispersion correction for bilayer-TiS2, bilayer-MoS2, and the TiS2/MoS2 heterostructure closely matched the values obtained using DFT, as shown in Table S1 of the supplementary material.

We utilized non-equilibrium molecular dynamics (NEMD) simulations with the MLIPs to estimate the thermal conductivity of the bilayer-MoS2, bilayer-TiS2, and TiS2/MoS2 heterostructures. Figure 3(a) illustrates the schematic model setup for the NEMD simulation. This method involves selecting two regions of atoms as the heat source and heat sink, respectively, where the momentum exchange occurs between the hottest atom in the sink and the coldest atom in the source. This exchange ensures the conservation of both total energy and momentum.41 This external heat flux is balanced by a heat current flowing from the source to the sink within the configuration, resulting in a measurable temperature gradient. Prior to applying heat flux, the initial configuration is equilibrated at a temperature of 300 K for 2 × 105 steps, with a time step of 1.0 fs using a constant pressure and constant temperature (NPT) molecular dynamics (MD) simulation. The configuration is then transitioned to the constant volume and constant energy (NVE) ensemble in order to maintain energy conservation.41 A heat flux is applied to the configuration by adding a small amount of heat dq to the hot region and removing it from the cold region continuously, as depicted in Fig. 3(a).

FIG. 3.

(a) NEMD setup for determining the thermal conductivity of the bilayer system having opposite hot and cold regions. (b) Temperature profile of the TiS2/MoS2 heterostructure. (c) Inverse of thermal conductivity for the bilayers as a function of the inverse of the distance between the hot and cold regions (l). The thermal conductivity at infinite sample length can be obtained by linear extrapolating to 1/l = 0.

FIG. 3.

(a) NEMD setup for determining the thermal conductivity of the bilayer system having opposite hot and cold regions. (b) Temperature profile of the TiS2/MoS2 heterostructure. (c) Inverse of thermal conductivity for the bilayers as a function of the inverse of the distance between the hot and cold regions (l). The thermal conductivity at infinite sample length can be obtained by linear extrapolating to 1/l = 0.

Close modal
The heat flux (Q) is determined by
(1)
where A is the cross-sectional area and dt indicates the time step. We performed 2 × 105 time steps to achieve a steady state condition and an additional 106 time steps to obtain the time-averaged temperature profile. This profile is then utilized to calculate the temperature gradient (Tl) with respect to the distance (l) along the direction of heat transfer, as depicted in Fig. 3(b). The thermal conductivity is calculated as
(2)
The thermal conductivity predicted by non-equilibrium molecular dynamics (NEMD) often exhibits a significant size effect caused by phonon scattering at the ends of the sample, particularly when the sample size is smaller than the mean-free path.41,44,45 The thermal conductivity at infinite sample length is determined using an inverse fitting approach.41 We selected the bilayer-MoS2, bilayer-TiS2, and TiS2/MoS2 heterostructure samples with a constant width of 51 Å and varied the length within the range of 118–189 Å. As a result, the distance between the hot and cold regions increased from 57 to 98 Å. The thermal conductivity at infinite length (λ) is derived through linear fitting the calculated thermal conductivities at various sample sizes using the following equation,41 as illustrated in Fig. 3(c):
(3)
where λ is the thermal conductivity at infinite length, L refers to the mean free path, and l represents the distance between the hot and cold regions. We use a similar approach to determine the thermal conductivity for both bilayer-MoS2 and bilayer-TiS2 samples.

Upon intercepting the fitting line along the vertical axis, the λ of bilayer-MoS2 was approximated to be ∼11.41 W/mK. Using a similar approach, the λ of bilayer-TiS2 was calculated to be 15.40 W/mK. Interestingly, the TiS2/MoS2 heterostructure exhibited a λ of 9.12 W/mK. This minor decrease in λ in comparison with bilayer-MoS2 and bilayer-TiS2 could be attributed to the lattice-mismatch between the two materials. To the best of our knowledge, there have been no previous efforts in the literature to predict the thermal conductivity of bilayer-TiS2 and TiS2/MoS2 heterostructures using either theoretical or experimental approaches. Therefore, to partially assess the validity of our machine learning approach, we employed the MTP technique to predict the size-independent thermal conductivity (λ) of bilayer-MoS2, for which theoretical and experimental values are available in the literature. Our calculated λ value for bilayer-MoS2 is in excellent agreement with an experimentally measured value, as shown in Table S2 of the supplementary material, and is, to some extent, in agreement with other theoretical predictions using other interatomic potentials. The thermal conductivities of bilayer-MoS2, bilayer-TiS2, and TiS2/MoS2 heterostructures are considerably higher than that of graphite (∼1.0 W/mK at ambient temperature46), an anode material typically employed in batteries. This indicates that these materials have the potential for use in thermal management and energy storage applications due to enhanced thermal properties.

This paper investigated the in-plane thermal conductivity of bilayer-TiS2, bilayer-MoS2, and TiS2/MoS2 heterostructures using an efficient machine learning interatomic potentials (MLIP) scheme. We utilized moment tensor potentials (MTPs) to accurately define short-range (intralayer) interactions. These MTPs were trained over relatively inexpensive computational DFT datasets. In addition, long-range (interlayer) van der Waals interaction parameters were calibrated using the D3-dispersion correction technique. We conducted non-equilibrium molecular dynamics simulations to determine the in-plane thermal conductivity of the heterostructures. This work illustrated how the thermal conductivities of bilayer-MoS2, bilayer-TiS2, and TiS2/MoS2 heterostructures are significantly greater than those of other materials used in battery energy storage applications, such as graphite. This paper’s methodology for deriving interatomic potentials using MTP can be extended to handle various conditions, such as grain boundaries and lattice-mismatched bilayer supercells. Therefore, the proposed methodology is a reliable tool for estimating the thermal conductivity of other van der Waals heterostructures.

The supplementary material contains the following plots: adhesive and strain energies of the TiS2/MoS2 heterostructure; kinetic energy variation of bilayer-MoS2 as a function of time (fs); 2D histogram plots showing the error distributions of total energies and forces for both the training and test sets of bilayer-MoS2, bilayer-TiS2, and the TiS2/MoS2 heterostructure; interlayer energy vs interlayer distance plot obtained using the MTP + D3-dispersion correction technique for three distinct bilayers; interlayer distances (Å) of the structures obtained utilizing the MLIP and DFT approaches; and estimated thermal conductivity values (λ∞) for bilayer-MoS2 using different approaches.

The authors acknowledge the Natural Sciences and Engineering Research Council of Canada (NSERC) and the New Frontiers in Research Fund-Exploration program for their financial support. The DFT simulations were performed using the Niagara Supercomputer from SciNet, which is supported by the Canada Foundation for Innovation, overseen by Digital Research Alliance of Canada, as well as the Ontario Research Fund—Research Excellence, the Government of Ontario, and the University of Toronto, Canada. The authors also acknowledge CMC Microsystems for the provision of products and services that facilitated this research, including the Synopsys QuantumATK Academic Software.

The authors have no conflicts to disclose.

A. K. Nair: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal). C. M. Da Silva: Formal analysis (equal); Project administration (equal); Writing – review & editing (lead). C. H. Amon: Funding acquisition (equal); Resources (equal); Supervision (equal); Writing – original draft (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

1.
H.
Song
,
J.
Liu
,
B.
Liu
,
J.
Wu
,
H.-M.
Cheng
, and
F.
Kang
,
Joule
2
,
442
(
2018
).
2.
Y.
Fu
,
J.
Hansson
,
Y.
Liu
,
S.
Chen
,
A.
Zehri
,
M. K.
Samani
,
N.
Wang
,
Y.
Ni
,
Y.
Zhang
,
Z.-B.
Zhang
et al,
2D Materials
7
,
012001
(
2020
).
3.
A. L.
Moore
and
L.
Shi
,
Mater. Today
17
,
163
(
2014
).
4.
Y.
Hu
,
L.
Zeng
,
A. J.
Minnich
,
M. S.
Dresselhaus
, and
G.
Chen
,
Nat. Nanotechnol.
10
,
701
(
2015
).
5.
S.
Lee
,
D.
Broido
,
K.
Esfarjani
, and
G.
Chen
,
Nat. Commun.
6
,
6290
(
2015
).
6.
C.
Da Silva
,
F.
Saiz
,
D. A.
Romero
, and
C. H.
Amon
,
Phys. Rev. B
93
,
125427
(
2016
).
7.
S.
Ghosh
,
I.
Calizo
,
D.
Teweldebrhan
,
E. P.
Pokatilov
,
D. L.
Nika
,
A. A.
Balandin
,
W.
Bao
,
F.
Miao
, and
C. N.
Lau
,
Appl. Phys. Lett.
92
,
151911
(
2008
).
8.
A.
Jain
and
A. J.
McGaughey
,
Sci. Rep.
5
,
8501
(
2015
).
10.
X.
Gu
,
B.
Li
, and
R.
Yang
,
J. Appl. Phys.
119
,
085106
(
2016
).
11.
Y.
Liu
,
S.
Zhang
,
J.
He
,
Z. M.
Wang
, and
Z.
Liu
,
Nano-Micro Lett.
11
,
13
(
2019
).
12.
Y.
Luo
,
M.
Li
,
H.
Yuan
,
H.
Liu
, and
Y.
Fang
,
npj Comput. Mater.
9
,
4
(
2023
).
13.
W.
Li
,
J.
Carrete
,
N. A.
Katcho
, and
N.
Mingo
,
Comput. Phys. Commun.
185
,
1747
(
2014
).
14.
A.
Taheri
,
C.
Da Silva
, and
C. H.
Amon
,
J. Appl. Phys.
123
,
215105
(
2018
).
15.
A.
Taheri
,
C.
Da Silva
, and
C. H.
Amon
,
Phys. Rev. B
99
,
235425
(
2019
).
16.
Z.
Fan
,
L. F. C.
Pereira
,
H.-Q.
Wang
,
J.-C.
Zheng
,
D.
Donadio
, and
A.
Harju
,
Phys. Rev. B
92
,
094301
(
2015
).
17.
B.
Mortazavi
,
E. V.
Podryabinkin
,
I. S.
Novikov
,
S.
Roche
,
T.
Rabczuk
,
X.
Zhuang
, and
A. V.
Shapeev
,
J. Phys.: Mater.
3
,
02LT02
(
2020
).
18.
S.
Arabha
,
Z. S.
Aghbolagh
,
K.
Ghorbani
,
S. M.
Hatam-Lee
, and
A.
Rajabpour
,
J. Appl. Phys.
130
,
210903
(
2021
).
19.
B.
Mortazavi
,
I. S.
Novikov
,
E. V.
Podryabinkin
,
S.
Roche
,
T.
Rabczuk
,
A. V.
Shapeev
, and
X.
Zhuang
,
Appl. Mater. Today
20
,
100685
(
2020
).
20.
A. K.
Nair
,
C.
Da Silva
, and
C.
Amon
,
J. Chem. Phys. C
127
,
9541
(
2023
).
21.
S.
Smidstrup
,
T.
Markussen
,
P.
Vancraeyveld
,
J.
Wellendorff
,
J.
Schneider
,
T.
Gunst
,
B.
Verstichel
,
D.
Stradi
,
P. A.
Khomyakov
,
U. G.
Vej-Hansen
et al,
J. Phys.: Condens. Matter
32
,
015901
(
2019
).
22.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
23.
J. R.
Chelikowsky
and
S. G.
Louie
,
Phys. Rev. B
29
,
3470
(
1984
).
24.
A.
Nair
,
C.
Da Silva
, and
C.
Amon
,
Appl. Surf. Sci.
600
,
154164
(
2022
).
25.
A.
Nair
,
C.
Da Silva
, and
C.
Amon
,
J. Appl. Phys.
133
,
064302
(
2023
).
26.
T.
Lorenz
,
D.
Teich
,
J.-O.
Joswig
, and
G.
Seifert
,
J. Chem. Phys. C
116
,
11714
(
2012
).
27.
L.
Kou
,
T.
Frauenheim
, and
C.
Chen
,
J. Chem. Phys.
4
,
1730
(
2013
).
28.
H. J.
Monkhorst
and
J. D.
Pack
,
Phys. Rev. B
13
,
5188
(
1976
).
29.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
,
J. Chem. Phys.
132
,
154104
(
2010
).
30.
J.
Jung
,
C.
Kobayashi
, and
Y.
Sugita
,
J. Chem. Phys.
148
,
164109
(
2018
).
31.
Y.
Cai
,
J.
Lan
,
G.
Zhang
, and
Y.-W.
Zhang
,
Phys. Rev. B
89
,
035438
(
2014
).
32.
J.
Yuan
,
Z.
Lv
,
Q.
Lu
,
Y.
Cheng
,
X.
Chen
, and
L.
Cai
,
Solid State Sci.
40
,
1
(
2015
).
33.
B.
Javvaji
,
X.
Zhuang
,
T.
Rabczuk
, and
B.
Mortazavi
,
Adv. Energy Mater.
12
,
2201370
(
2022
).
34.
I. S.
Novikov
,
K.
Gubaev
,
E. V.
Podryabinkin
, and
A. V.
Shapeev
,
Mach. Learn.: Sci. Technol.
2
,
025002
(
2021
).
35.
T.
Mueller
,
A.
Hernandez
, and
C.
Wang
,
J. Chem. Phys.
152
,
050902
(
2020
).
36.
A.
Takahashi
,
A.
Seko
, and
I.
Tanaka
,
J. Chem. Phys.
148
,
234106
(
2018
).
37.
N.
Wang
and
S.
Huang
,
Phys. Rev. B
102
,
094111
(
2020
).
38.
A. M.
Goryaeva
,
J.
Dérès
,
C.
Lapointe
,
P.
Grigorev
,
T. D.
Swinburne
,
J. R.
Kermode
,
L.
Ventelon
,
J.
Baima
, and
M.-C.
Marinica
,
Phys. Rev. Mater.
5
,
103803
(
2021
).
39.
G.
Dhaliwal
,
P. B.
Nair
, and
C. V.
Singh
,
npj Comput. Mater.
8
,
7
(
2022
).
40.
H.
Niu
,
J.
Zhao
,
H.
Li
,
Y.
Sun
,
J. H.
Park
,
Y.
Jing
,
W.
Li
,
J.
Yang
, and
X.
Li
,
Comput. Mater. Sci.
218
,
111970
(
2023
).
41.
Z.
Ding
,
J.-W.
Jiang
,
Q.-X.
Pei
, and
Y.-W.
Zhang
,
Nanotechnology
26
,
065703
(
2015
).
42.
L.
Goerigk
,
Non-Covalent Interactions in Quantum Chemistry and Physics
,
195
(
Elsevier
,
2017
).
43.
J.
Moellmann
and
S.
Grimme
,
J. Phys. Chem. C
118
,
7615
(
2014
).
44.
P. K.
Schelling
,
S. R.
Phillpot
, and
P.
Keblinski
,
Phys. Rev. B
65
,
144306
(
2002
).
45.
Y.
Wang
,
A. K.
Vallabhaneni
,
B.
Qiu
, and
X.
Ruan
,
Nanoscale Microscale Thermophys. Eng.
18
,
155
(
2014
).
46.
B.
Koo
,
P.
Goli
,
A. V.
Sumant
,
P. C.
dos Santos Claro
,
T.
Rajh
,
C. S.
Johnson
,
A. A.
Balandin
, and
E. V.
Shevchenko
,
ACS Nano
8
,
7202
(
2014
).