Intermediates of forming transition metal dichalcogenide heterostructures revealed by machine learning simulations

Domain: physics
Published: 2026-02-23
Editorial outcome: major_revision

Abstract
Two-dimensional (2D) transition metal dichalcogenide (TMD) van der Waals heterostructures (vdWHs) hold promise for high-performance electronics, but their large-scale synthesis remains limited by size constraints and alloying contaminations. Recently, a two-step vapor deposition method was reported for growing wafer-size TMD vdWHs with minimal impurities. In this study, we develop a machine learning potential (MLP) that captures the atomic-scale dynamic growth process of bilayer MoS2/WS2 vdWHs under feasible growth conditions. Our simulations uncover a crucial metastable SMMS (M = Mo or W) intermediate structure that facilitates metal atom swap and alloying. Eliminating the alloying contamination requires preventing the embedding of bare metal atoms. The results also show that the SMMS structure exhibits favorable electronic properties and emerges as a low Schottky barrier contact electrode for MoS2 field-effect transistors (FETs).

Introduction
2D TMDs have been extensively studied due to their semiconducting band gaps and high carrier mobility1–4, strong nonlinear optical response5,6, and ease of layer stacking assembly. The integration of vdWHs by pristine TMD layers can further enhance their properties, enabling various potential applications in microelectronics, optoelectronics, and nonlinear optics7–14.

Despite their promise, the controlled growth of TMD vdWHs faces numerous challenges. While commonly used mechanical assembly methods can achieve high-quality TMD vdWHs with atomically sharp interfaces15,16, they often struggle to achieve wafer-sized dimensions and can be prohibitively expensive. In contrast, chemical vapor deposition (CVD) has achieved significant success in growing wafer-sized monolayer TMD12,17,18. However, growing TMD vdWHs still faces limitations in size constraints and a strong tendency toward alloy formation19,20. Among these scalable methods, metal-organic chemical vapor deposition (MOCVD) has emerged as a promising alternative for growing large-size TMD heterojunctionswith improved interfaces20–22.

Recently, it was reported that a two-step vapor deposition process with a high-to-low temperature strategy was used to synthesize wafer-sized TMD vdWHs23. In this method, a monolayer WS2 film was first achieved by depositing a W film on a sapphire substrate and sulfurizing it at 900 °C, the highest temperature among the four stacked materials. Next, a Mo film was deposited on the WS2 monolayer via magnetron sputtering and sulfurized at 800 °C to form MoS2. Subsequently, a Nb film was deposited on the WS2/MoS2 film via magnetron sputtering and selenized at 700 °C to form NbSe2. Finally, PtTe2 was grown on the WS2/MoS2/NbSe2 vdWHs at 350 °C, the lowest temperature in the sequence. The final structure was a wafer-sized vdWHs several centimetres in size, consisting of four layers: WS2/MoS2/NbSe2/PtTe2. The authors further proposed that, during the metal deposition process, the metal atom films coat the TMD surface23.

Compared with experimental trials that involve many parameters and various growth conditions, accurate atomic simulations can certainly provide an in-depth understanding of the growth of TMD vdWHs. The growth of TMD vdWHs is always accompanied by rapid formation and breaking, involving various intermediate motifs with complex chemical bonds, such as metallic bonds, covalent bonds, both covalent and ionic characters of Mo-S bonds, and layered vdW forces. Density functional theory (DFT) is well suitable for describing the complex chemical bonding, but is limited to afford large-scale simulations. In contrast, classical molecular dynamics (MD) can simulate large systems but lacks the capability to capture complex chemical bond recombination. Although various machine learning potentials (MLPs) have been developed and reported to simulate large systems with accuracy comparable to DFT24–29, only a few MLPs are capable of describing intricate growth processes of TMD vdWHs.

In this study, we developed an MLP trained on extensive DFT data and utilizing a revised equivariant graph neural network implemented in the NequIP package24. The MLP was then implemented into MD (MLP-MD), enabling both large-scale and accurate simulation of the MoS2/WS2 vdWHs growth process. Our results indicate that a bare metal atomic layer is unstable on TMD layers; it spontaneously sinks into the S layer, forming a crucial intermediate structure (SMMS) with high stability. This behavior facilitates the exchange between Mo and W atoms, thereby revealing the atomic-scale mechanism of TMD alloys formation. To grow pristine TMD vdWHs, it is essential to suppress the SMMS formation, which can be achieved by preventing bare metal atom adsorption on pre-existing TMD layers. On the other hand, the SMMS structure serves as an ideal metallic electrode with a low Schottky barrier height, enabling MoS2 integrated circuits through the planar deposition of metal atoms on uncovered MoS2.

Results
MLP Development and DFT Validation
To enhance the complexity and versatility of our MLP model, we designed a diverse dataset encompassing a wide range of Mo, W, S, and mixed structures. We also included Se atomic structures to broaden the scope of potential research applications. Through an iterative learning process, we carefully selected and balanced the weights of different structural configurations. The initial training set was derived from on-the-fly machine learning MD simulations30 and was then continuously enriched with new structures throughout successive iterations. The final dataset included various TMD layers, MoS2/WS2 vdWHs, MoS2 or WS2 with homogeneous and heterogeneous metal clusters, S clusters, diverse MoWS alloy configurations, and growth intermediate structures. This comprehensive approach resulted in a dataset of approximately 26,000 DFT-calculated entries. Detailed procedures for both the initial training and subsequent iterative refinement are provided in the Supplementary Information.

Deposition Dynamics of Mo Atoms on MoS2
Before investigating TMD vdWHs, the developed MLP-MD model successfully simulated the growth of bilayer MoS2 through a two-step vapor deposition method23. First, the adsorbed single Mo atom was unstable on the MoS2 layer and, at the typical growth temperature of 1100 K, quickly sank beneath the S atom layer within tens of picoseconds. DFT calculations confirm that this embedding process releases 1.45 eV of energy. Figure 2c shows a significant drop in energy at around 50 ps, indicating that the Mo atom has been embedded into the MoS2 layer. Previous studies have shown that single or paired metal atoms suspended on bare surfaces exhibited significant catalytic enhancement33–36. This indicates that the configuration of metal atoms needs to be carefully considered in single-atom catalysis.

Subsequently, individual Mo atoms were simulated to be continuously sputtered onto the existing MoS2 layer with a kinetic energy of 0.12 eV. To simulate the stability during the metal atom deposition and pre-heating stage, simulations were performed at 900 K. The deposition process of Mo atoms that resulted in 0.25 monolayers (MLs) and 1.0 ML demonstrated the behavior of Mo atoms on the MoS2 layer. Figure S3 depicts the evolution of Mo-Mo bonds during continuous deposition, indirectly reflecting the embedding process. We observed that the number of Mo-Mo bonds increased gradually with continued deposition and reached saturation once deposition stopped. As shown in Fig. S4, when Mo atoms were deposited at higher kinetic energies, such as 1 and 5 eV, the Mo atoms were still embedded in MoS2. This indicates that the embedding behavior of Mo atoms is not sensitive to kinetic energy. Throughout the MLP-MD simulation, no Mo atoms formed the MoSMoS structure. Instead, all deposited Mo atoms spontaneously sank into the MoS2 layer. These findings suggest that under practical deposition conditions, Mo atoms preferentially form an SMoMoS structure. This SMoMoS structure was then constructed and optimized using DFT, revealing that it exhibits an energy 0.080 eV/atom lower than that of the directly formed embedded structure from MD simulations.

Heterogeneous Deposition Dynamics of Mo Atoms on WS2
Following the two-step vapor-deposition process23, we continuously deposited Mo atoms onto the WS2 monolayer in MLP-MD simulations, leading to the formation of MoS2/WS2 vdWHs. Figure 3(a–c) display snapshots from the simulation at 0.06 ns, 0.15 ns, and 1.1 ns, showing that the Mo atoms do not remain on the WS2 surface but are embedded within the WS2 monolayer. Figure S8 shows the evolution of Mo/W-Mo/W bonds during continuous Mo deposition, indirectly reflecting the degree of Mo atom insertion. The number of these bonds gradually increases with deposition and saturates once deposition stops. Importantly, during the simulation, Mo and W atoms in the SMoWS intermediate structure layer are able to exchange, resulting in alloying. At 300 K, simulation results show that while Mo atoms continue to embed into the WS2 monolayer to form the SMoWS structure, no Mo-W atomic exchange occurs over a 4.5 ns trajectory, indicating that alloying is kinetically suppressed under low-temperature conditions.

Sulfurization of Intermediate Structures
We now pose the question of what happens when S is deposited on the SMoMoS and alloyed SMMS intermediate phases during the second step. For SMoMoS, a sufficient amount of S atoms was further deposited on the top surface of the SMoMoS intermediate phase. The simulation shows that the SMoMoS intermediate phase was initially very stable. After 531 ps, we observed that S atoms penetrated the SMoMoS structure and pulled some Mo atoms to the surface, gradually forming a bilayer of MoS2. As shown in Fig. S12a, there are very few initial Mo-S bonds because the S atoms deposited on the SMMS surface cannot bond with Mo. As the MD simulation progresses, some Mo atoms are pulled to the surface to bond with the deposited S atoms, leading to an increase in the number of Mo-S bonds. Between 700 and 1000 ps, a large number of Mo atoms are pulled to the surface, resulting in a substantial increase in Mo-S bonds. The number of Mo-S bonds tends to be saturated after about 1.2 ns, when a bilayer MoS2 is fully formed.

Similarly, for the alloyed SMMS intermediate phase, a sufficient amount of S atoms was further deposited on the top surface. As shown in Fig. 4a, the alloyed SMMS intermediate phase was also very stable initially. After 608 ps of simulation, S atoms penetrated the SMMS structure and pulled metal atoms from the top surface. However, in the alloyed intermediate phase, S atoms did not selectively extract either Mo or W atoms, instead, they pulled both Mo and W atoms from the upper layer to the surface. Therefore, the resulting structure was not a distinct MoS2/WS2 vdWHs, but rather an alloyed MoxW1-xS2/Mo1-xWxS2 vdWHs. As shown in Fig. S12b, the evolution of Mo/W-S bonds clearly illustrates this process. According to the experimental results of Zhou et al.23, the successful synthesis of non-alloyed MoS2/WS2 heterostructures suggests that the formation of the SMMS intermediate state may have been avoided through process control during the experiment.

Electronic Properties and Device Potential of Intermediate Structures
As previously discussed, the SMMS intermediate structure impedes the growth of non-alloyed MoS2/WS2 vdWHs. Figure 5a and Fig. 5b show the electronic band structures. Both SMoMoS and SMMS exhibit metallic properties.

We further investigated the contact characteristics of these metallic intermediate structures with the semiconductor MoS2. In both cases, the contact between the metallic intermediate structures and the semiconductor MoS2 forms a p-type Schottky contact. The Schottky-Mott limit predicts p-type SBHs of 0.55 eV for SMoMoS and 0.69 eV for alloyed SMMS. For conventional metal electrodes (e.g., Ti, Cr, Au, Pd) interfaced with MoS2, the predicted SBH ranges from 0.56 eV to 1.86 eV42. However, Fermi-level pinning induces strong n-type behavior at MoS2 interfaces. Experimentally observed p-type SBHs range from 1.56 to 1.75 eV. We constructed SMMS-MoS2 and SMoMoS-MoS2 interfaces. Compared to Schottky-Mott predictions, we observed significantly reduced Fermi-level pinning. The p-type SBH increased by only approximately 0.1 eV in the same theoretical framework. Moreover, continuous SBH tuning can be achieved by adjusting metal deposition parameters, enabling tailored contact properties for electronic or optoelectronic applications. These results suggest that SMMS and SMoMoS are promising candidates for MoS2 FET electrodes.

Methods
DFT
DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP)44,45 version 6.3.0. A plane wave basis set was employed, utilizing the projector augmented wave (PAW) method and standard pseudopotentials. The Perdew-Burke-Ernzerhof (PBE)46 exchange-correlation functional within the Generalized Gradient Approximation (GGA)47 was employed to compute the electronic structure and energy. The DFT-D3 (Becke-Johnson)48,49 van der Waals correction was chosen to account for dispersion interactions. A plane-wave cutoff energy of 500 eV and no symmetry constraints were applied. To ensure accuracy, the electronic self-consistent loop was converged to a tolerance of 10−5 eV. Gaussian smearing with a width of 0.05 eV was employed to facilitate the convergence of calculations. Spin-polarized calculations were performed. For all periodic structures, a center-symmetric k-point mesh with a density of 0.3 Å−1 was employed. The electronic band structure calculations were enhanced through the implementation of the HSE06 hybrid functional for both band structures and electrostatic potentials. A Fermi-Dirac smearing parameter of 0.2 eV was employed in the calculations. MD simulations were conducted using the on-the-fly machine learning potential from the VASP package, extracting frames from DFT calculations as data to accelerate the construction of the initial training set.

MLP
The MLP was trained using the NequIP24 framework, which implements E(3)-equivariant graph neural networks. This network provides enhanced stability for MD simulations relative to other networks, enabling extended stable simulations. The model utilized 12 radial basis functions and a maximum angular momentum of 2 for the interatomic edges. The hidden layers included irreducible representations as follows: 128×0e, 64×0o, 128×1o, 64×1e, 32×2o, and 32×2e. Element information was embedded as 128-dimensional vectors. Five graph convolution layers were implemented, including self-connections and residual connections. The cutoff distance for constructing the graph was set at 6 Å. The loss function included both energy per atom and atomic force terms, each contributing equally to the total loss. Each loss term used mean squared error. During the initial training, we used a learning rate of 0.005 and a batch size of 5. For the iterative learning process, we employed a lower learning rate of 10−4 to fine-tune the model. To enhance the model’s robustness, we added a constant repulsive term to the output layer: E_repulsive = 0.5 * k * (r0/r)^12 - k * (r0/r)^6. This approach enhanced the model’s stability, with r0 = 1.8 Å in this study. The iterative learning process allowed us to continuously refine the model by incorporating new structures and configurations, ensuring comprehensive coverage of the complex structural landscape involved in TMD heterostructure growth.

MD
MD simulations were conducted using the Large-scale Atomic/Molecular Massive Parallel Simulator (LAMMPS)50, utilizing our developed MLP. We employed the bin algorithm for neighbor list construction with a cutoff radius of 6.0 Å and a skin distance of 2.0 Å. The neighbor list was rebuilt only when at least one atom had moved beyond half of the skin distance threshold. Simulations were conducted within the canonical ensemble (NVT) employing a Nosé-Hoover chain thermostat51,52, with a temperature damping time set to 0.1 ps. An integration timestep of 1.0 fs was employed to ensure performance while maintaining simulation stability. Initial velocities were sampled from a Gaussian distribution. The degrees of freedom contributing to the system temperature were dynamically updated with the deposition of new atoms. The simulations employed an ideal Morse substrate (V(r) = D0(1 - e^-α(r-r0))^2) with parameters D0 = 0.2 eV, α = 1.5, and r0 = 3.5 Å. This wall mimics the weak van der Waals interaction (<0.1 eV/atom) between Mo/W atoms and oxide substrates, such as sapphire, thus providing only mechanical support while avoiding spurious interfacial chemistry. For Mo atom deposition, we used a 4√3×7 MoS2 or WS2 supercell at 900 K. This temperature was chosen to represent the thermal state during the pre-heating/baking stage and the ramp-up process prior to full sulfurization. The SMoMoS and SMMS sulfurized process simulations were performed using a 4√3×7 supercell at 1100 K, aligning with the experimental sulfurization temperatures (approx. 1073 K). The diffusion behavior of Mo-S clusters on the MoS2 surface was also investigated using a 4√3×7 MoS2 supercell at 1100 K. For the simultaneous Mo and S deposition to grow the second MoS2 layer, we employed a larger 6√3×8 MoS2 or WS2 supercell at 1100 K. The atomic deposition was achieved through the fix deposit command in LAMMPS, which introduced a Mo atom with a downward velocity of 5 Å/ps at random positions within the top region of the simulation box every 6000 steps. The initial kinetic energy of deposited atoms was set to 0.12 eV. To ensure the reproducibility of the non-equilibrium processes, all key qualitative observations (including the sinking of Mo atoms and the sulfurization-induced extraction) were confirmed in three independent MD simulations initialized with different random velocity seeds.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In their manuscript, Zhao and coworkers conduct extensive atomistic simulations to explore the growth process of 2D transition metal dichalcogenide (TMD) van der Waals heterostructures (vdWHs). They employ a machine learning potential (MLP) based on the NequIP architecture, trained on DFT calculations, using an iterative training approach to generate the training and test sets. Non-equilibrium molecular dynamics simulations are utilized to examine the stability and growth of MoS2/WS2 heterostructures. The study identifies two distinct growth mechanisms: a sulfur-deficient growth via an SMeMeS (Me=W/Mo) intermediate, and a sulfur-rich growth where the MeS2 structure remains intact while MeSx clusters adsorb and agglomerate to form the next MeS2 layer. The former leads to alloy formation, while the latter prevents metal mixing in the vdWHs. These findings help rationalize a recently published \"two-step vapor-deposition process.\" Additionally, the authors propose that the SMeMeS intermediate could be used to create an atom-scale p-type Schottky contact, potentially significant for technological applications.\n\nThe manuscript provides valuable scientific insights into the growth process and should be of significant interest to the 2D-material community. However, there are several shortcomings that require major revision. Firstly, the explanation of the training procedure lacks critical information necessary for reproducibility. Secondly, the MLP validation does not ensure the required accuracy for the results obtained. Thirdly, the non-equilibrium simulations are poorly explained and seem to lack a systematic approach, which raises questions about their robustness. Additionally, the writing is sometimes vague and difficult to follow, this applies to both the language and the organization of the results. I recommend that the authors address these issues to provide a much needed confidence in the presented results (which I believe at the moment to be insufficient) and improve readability.\n\nMajor comments:\n1. Details of the iterative training procedure for machine learning potential are missing. I couldn't find any information on how they prepared the many different configurations in the training set – and what structures were specifically included. I believe that this also needs motivation since a repository alone leaves too many open questions. The current information on the procedures “on-the-fly MLPMD in VASP” and “Monte Caro method” are not clear and do not offer information on structural data. Also the data-selection procedure (each iteration with 50% of new data being included) is unclear since “new” data is not defined. Note, that 50% of 100,000 configurations per iteration clearly exceeds a final data set of 26,000 entries\n2. The MLP validation consists of fairly meaningless RMSE’s and some MD simulations validating(?) the MLP and AIMD validating the stability predictions. These calculations have limited certainty of the accuracy. What I believe is much more important is, e.g. comparing the predicted relative stability of optimized and non-optimized structures (e.g. SMeMeS vs. MeS2) via the MLP vs. DFT.\n3. The non-equilibrium MD simulations are also unclear. How exactly are new species introduced during the MD simulations? Why are they introduced with low velocities. Furthermore, why didn’t the authors collect statistics and clearly quantify e.g. mean times of incorporation. The efficient MLP would certainly allow for this. Currently, all results are qualitative observations without any quantitative indication.\n4. The discussion on of the growth mechanism should be more organized. Each paragraph appears as a loose description of what simulation was performed and what was seen lacking motivation that would help the reader gaining the very important scientific insight.\n5. In illustrate the latter point, data from various simulations appears sometimes erratically combined. An example is the discussion of the stability of the bilayer system. If I understood correctly, the structure for MoSMoS, SMoMoS, SMo3S, SMo4S (Fig2. (f-i)) are obtained by geometry optimization, not by MD simulation. Did the authors observe the SMoMoS (Fig.2 (g)) structure during the MD simulation starting from the MoSMoS structure (Fig.2 (f))? If not, the statement of “all deposited Mo atoms spontaneously sank into the MoS2 layer, tending to form the intriguing SMoMoS structure, as shown in Fig. 2(g).” and “In addition to the SMoMoS structure formed by depositing 1 ML of Mo atoms” might be misleading, because as seen in Fig. S5 (b) it is less likely that the formation of SMoMoS structure from MoSMoS structure happens.\nAlso, in this discussion, the authors showed a lot of data:\n○ energy landscape of MoSMoS, SMoMoS, SMo3S, SMo4S structures (Fig2. (f-i))\n○ phonon band structures (Fig2. (l))\n○ tendency that Mo atoms sink beneath the S layer during the MLP-MD simulation (Fig2. (e))\n○ tendency that Mo atoms sink beneath the S layer during the AIMD simulation (Fig S5(b))\n○ stability of SMoMoS during the MLP-MD simulation (Fig S3 (a)(d))\n○ stability of SMoMoS during the AIMD simulation (Fig S5 (c))\nAs such, the data from geometry optimization and MD simulations using DFT and MLP are given, and it is hard to understand how they complement each other in the current version of the manuscript. This should be more clarified for better discussions.\n\nMinor comments:\n● Please check the numbering for the figures. For example, Fig. S4 comes earlier than Fig. S3, which is not kind to readers.\n● Small comments on the figures:\n○ Fig.2 (j) it is stated that “with green dots” but I couldn't find the green point\n○ Fig.2 (c) what is the red line and blue line?\n○ For the caption for Fig.2, the labels (a) to (e) come after the caption but (f) to (l) come before the caption, which is, I think, very confusing.\n○ Fig.5 (b) SMOWS should be SMoWS\n● For the phonon band structures, the data from DFT and MLP are plotted together on the same figure, but I suggest to show only DFT results. The authors stated in page 3 “It is worth noting that the phonon dispersion simulated by the MLP accurately reproduced the DFT results, validating the accuracy of our MLP in terms of atomic forces.”, but some degeneracy at the gamma point are not well described by MLP, and it is hard to tell the accuracy of MLP from this figure. For Fig.3 (j-l), the authors don’t mention the differences between DFT and MLP, so I suggest removing the phonon band structures by MLP for better clarity.\n● Equations in Page 4\n○ Please define N and T.\n○ E_SMoWS is energy, then why is E formation energy?\n○ x should be “one Mo atom” not “one metal atom”?\n● Page 5: “Throughout the calculations, general precision was used”, it would be better to specify the precision. “general precision” is vague.\n● For the simulation given in video S4, the S layer looks floating above the SMoMoS layer. Is it constrained?\n● Page 5, what do you mean by “a specific rate of 516 ps per layer”?\n● For most of the systems, the system size is not clearly stated but it should be given.\n● Fig. S5(c) shows the AIMD simulation of SMoMoS structure until 10 ps, but what happens after 10 ps? Did the system retain the SMoMoS structure?\n● Again, are the structures of Fig.3 (g-i) obtained by geometry optimization and not appeared during the MD simulations?", "reviewer_2": "The manuscript by Zhao et al. presents a machine learning study investigating the atomistic details of van der Waals heterostructure (vdWH) formation during chemical synthesis. While the topic is inherently interesting and relevant, the paper does not meet the necessary criteria for publication due to several significant shortcomings.\n\nFirst, the language throughout the manuscript, particularly in the first two pages, is poorly constructed, making it difficult—at times even impossible—to discern the authors' intended meaning. This lack of clarity severely impacts the readability and comprehension of the work.\n\nSecond, the overall structure and logic of the manuscript are unclear. The progression of ideas is disorganized, making it challenging to follow the content from beginning to end. Additionally, several critical details are omitted, which undermines the rigor of the study. For instance:\n(i) The temperatures at which the calculations were conducted are not clearly specified.\n(ii) The substrates used for investigating the formation of vdWHs are not identified.\nThe authors motivate their theoretical study with experimental results published in Nature (621, 499, 2023). However, they fail to provide sufficient a meaningful comparison with their theoretical results with the experiment. Such a comparison is essential for evaluating the reliability of the study and its implications for the synthesis of pure vdWHs (without alloying).\n\nIn its current state, the manuscript lacks the clarity, structure, and critical details necessary for a proper review. Given these substantial deficiencies, I regret to conclude that the manuscript cannot be considered for publication and must be rejected.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "The authors have written an interesting paper on the formation of a stable intermediate state which they call SMMS that forms when depositing bare Mo atoms on MoS2 or WS2. While there isn’t much methodological innovation in the paper, at least as currently presented, it is a nice and fairly complex application of ML potentials to an important class of materials that sheds some light on the formation of MoS2/WS2 bilayers, namely that Mo must be deposited with sulfur to ensure that Mo doesn’t sink below the surface.\n\nI do have some concerns that I think should be addressed if the paper is to be published in Nature Communications.\n\nMajor concerns\n* The iterative learning approach presented in the paper is unclear. The results section states that the initial training set was generated with on-the-fly learning and then “continuously enriched with new structures throughout the iterations”. The methods and SI do not provide much additional information, beyond the number of iterations (17) and the number of “data points” per iteration (100,000). Are these data points entire structures or individual forces and energies? How exactly were these new structures generated? The SI alludes to a Monte Carlo method, but this is vague. In my opinion, the generation of the training set is central to the paper and will be of interest to others in the field. I suggest the authors describe their method in much more detail. As presented currently, it is not possible to reproduce it.\n* One of the central claims of the paper, appearing in the last sentence of the abstract, is that the SMMS phase that they identify is an “ideal electrode for MoS2 FETs”. This claim is weakly supported. The two pieces of evidence for the claim are 1) SMMS is a metal, as evidenced by the band structure in Fig. 5a-b, and 2) the p-type Schottky barriers seem to be in a reasonable range. To evaluate the claim that this structure is an “ideal electrode”, much more information would need to be presented, including comparison with existing electrodes, assessment of how easy SMMS is to synthesize relative to these other electrodes, etc. The authors also seem to suggest that the difference between 0.55 and 0.69 eV (the Schottky barriers for SMoMoS and SMMS respectively) is significant, but this seems hard to believe without additional context.\n\nMinor concerns\n* The authors claim in the SI that the accuracy of their potential is “outstanding”, but the RMSEs don’t strike me as particularly remarkable (151 meV/A on forces and 10.6 meV on energies). They also claim that their potential “outperforms several existing machine learning force fields in key aspects”, but do not elaborate. For these claims to stand the authors will need to provide direct comparisons with competing potentials. Otherwise, the claims should be removed or significantly weakened.\n* It seems the authors do not train on the stress tensor, which is unusual. Why were stress labels excluded from the training set?\n* In the results section, the authors claim that “the phonon dispersion of the SMoMoS structure showed no imaginary frequencies, further confirming that MoSMoS would spontaneously transform into SMoMoS.” Lack of imaginary frequencies only confirms that SMoMoS is stable, not that MoSMoS would transform into it.\n* In the results section, when comparing to AIMD simulations, the authors claim that the time for Mo to embed below the surface is significantly shorter than in MLMD simulations “due to the high chemical reactivity of the Mo atoms from vacuum”. I don't understand this explanation. Why isn’t this effect captured in the MLMD simulations, which are supposed to reproduce AIMD results?\n* The convex hull presented in Fig. 2j is unclear. The caption and main text highlight “green dots” that don’t appear in the figure. Four points are labelled with letters, but the letters aren’t referred to anywhere. And it’s hard to tell which point corresponds to SMoMoS.\n* The authors claim that “the Mo-S1 structure always floats on the surface without sinking”. But the simulation is quite short, only 1.4 ns according to Fig. 4. The authors should perhaps claim instead that “the Mo-S1 structure remains on the surface without sinking for the duration of the 1.4 ns simulation”.\n* The authors mention that “the 1T phase present in the early stages mostly disappears, indicating that MoS2 transitions from the 1T phase to the more stable 1H phase as growth proceeds.” It’s very difficult to see this in Fig. S8. The authors should consider pointing out where the 1T phase appears in the figure, and perhaps give schematic examples of the pristine 1H and 1T phases.\n* There are grammatical mistakes throughout the paper. For example, in the title, \"transition metal dichalcogenides heterostructures\" should be \"transition metal dichalcogenide heterostructures\". In the abstract, “a highly stable intermediate easily introduces” should be “a highly stable intermediate that easily introduces”, and “Eliminating the alloying contamination… is avoiding SMMS structure” should be “Eliminating the alloying contamination involves avoiding the SMMS structure”."}, "author_response": "The reviewer’s positive comments and valuable suggestions are highly appreciated.\nRegarding the main points reviewer made, we would like to respond as follows:\nComment 1:\nMajor comments:\n1. Details of the iterative training procedure for machine learning potential are missing. I\ncouldn't find any information on how they prepared the many different configurations in the\ntraining set – and what structures were specifically included. I believe that this also needs\nmotivation since a repository alone leaves too many open questions. The current information on\nthe procedures “on-the-fly MLPMD in VASP” and “Monte Caro method” are not clear and do\nnot offer information on structural data. Also the data-selection procedure (each iteration with\n50% of new data being included) is unclear since “new” data is not defined. Note, that 50% of\n100,000 configurations per iteration clearly exceeds a final data set of 26,000 entries.\n\nOur reply: Thank you for your valuable comments. The sampling procedure of the training\ndata and detailed explanations are shown below.\nQ1: Details of the iterative training procedure for machine learning potential are missing. I\ncouldn't find any information on how they prepared the many different configurations in\nthe training set?\nIt starts with generating an initial dataset from quantum mechanical calculations, which is then\nused to train a machine learning model. Next, we validate the reliability of the model and apply it\nto problems to be investigated, and if it becomes problematic, e.g., unstable molecular dynamics\nand poor phonon spectral correspondence, fine-tune it by adding relevant DFT constitutive data.\nThis data expansion and model retraining process are repeated iteratively until the model achieves\nsatisfactory accuracy. The final model can then be deployed for efficient and accurate large-scale\natomistic simulations.\n(1)what structures were specifically included?\nOur dataset includes:\n• Single-layer and multilayer TMD structures\n• MoS /WS heterostructures\n2 2\n• Configurations of TMD interactions with metal clusters (Mo, W), including embedded\nclusters\n• MoWS alloy configurations\n2\n• Intermediate structures may appear during the growth process, including the SMMS\ntype and MSMS.\n• Sulfur cluster and sulfur on TMD surface\nThe dataset has been uploaded to the Zenodo platform\n(https://doi.org/10.5281/zenodo.13777360).\n(2) The current information on the procedures “on-the-fly MLPMD in VASP” and “Monte Caro\nmethod” are not clear and do not offer information on structural data.\n“On-the-fly MLPMD in VASP” is an adaptive methodology developed by the VASP vendor\nthat dynamically trains machine-learned force field during MD simulations. At each simulation\nstep, the machine-learned force field predicts the energy and forces and determines whether the\nuncertainty exceeds a threshold through Bayesian error estimation. If the uncertainty exceeds a\npredetermined threshold, an initial DFT calculation is performed, and the resulting configuration\nis added to the training set to update machine-learned force field. If it is exceeded, an initial DFT\ncalculation is performed, and the resulting configuration is added to the training set to update the\nmachine-learned force field. This iterative process enables the machine-learned force field to learn\nthe potential energy surface in real time, concentrating computational effort on poorly predicted\nregions. Consequently, this approach accelerates the exploration of the potential energy surface\nduring the training set construction.\nIn our study, the “Monte Carlo (MC) method” is implemented through our self-developed\nprocedure, where we employ a random atom-swapping strategy to rapidly perturb the system, e.g.,\nby randomly swapping the metal atoms in the upper and lower layers of the SMMS, thus\naccelerating the exploration of the potential energy surface.\n(3) Also the data-selection procedure (each iteration with 50% of new data being included) is\nunclear since “new” data is not defined. Note, that 50% of 100,000 configurations per iteration\nclearly exceeds a final data set of 26,000 entries.\nIn each iteration cycle, we employed a specialized balanced sampling strategy to construct the\ntraining dataset. Specifically:\n1. “new” data definition: We randomly sampled 50,000 data points from newly generated\nconfigurations (some structural repetitions were inevitable during sampling), which is defined\nas new data.\n2. Simultaneously sampled 50,000 data points from the existing training data (some structural\nrepetitions were inevitable during sampling)\n3. Combined these to form a training batch of 100,000 data points\nThis 1:1 sampling ratio was designed with two key objectives:\n1. To ensure the model maintains accuracy on previously learned features while acquiring\nnew ones\n2. To enhance training stability by minimizing the impact of dataset size fluctuations\nSince our total training set size was smaller than 100,000 configurations, some structural\nrepetitions were inevitable during sampling. These repetitions effectively served as data weighting,\naiding the model in better learning key structural features. Through this process, we ultimately\nobtained approximately 26,000 unique structural configurations, encompassing all critical data\nfrom the initial training to the final iteration.\nDuring the iterative process, we employed a reduced learning rate (1e-4) for fine-tuning, which\nfurther ensured training stability. This combination of balanced sampling and cautious fine-tuning\nallowed us to effectively expand the training set's coverage while maintaining model stability and\naccuracy.\nAccordingly, we added these discussions in the revised Supplementary Information (Page 1):\n“In this study, we adopted a basic iterative training procedure that begins with generating an\ninitial dataset from quantum mechanical calculations, which is then used to train a machine\nlearning model. Next, we validate the reliability of the model and apply it to problems under\ninvestigation. If it becomes problematic, e.g., unstable molecular dynamics and poor phonon\nspectral correspondence, fine-tune it by adding relevant DFT constitutive data. This data\nexpansion and model retraining process are repeated iteratively until the model achieves\nsatisfactory accuracy. The final model can then be deployed for efficient and accurate large-scale\natomistic simulations. Then, we employed on-the-fly molecular dynamics (MD) simulations to\ngenerate the initial training set, effectively exploring unknown configuration spaces while\nminimizing sampling costs. Using VASP's on-the-fly MD framework, we implemented the\nfollowing parameters to optimize efficiency: a timestep of 3.5 fs ensured accurate atomic motion\nsampling, 20,000 maximized simulation steps, fixed unit-cell shapes, Nosé-Hoover thermostats,\nand a temperature range (500 K to 1300 K) covered typical growth conditions. The\nML_MB=2000 parameter constrained the number of local reference configurations in the on-the-\nfly MLP. Initial configurations comprised monolayer MoS /WS systems and their combinations\n2 2\nwith homogeneous/heterogeneous metal clusters.\nWe expanded the training set through ten iterative epochs with distinct sampling strategies:\n• Epoch 1: On-the-fly MD sampling of MSMS structure dynamics\n• Epoch 2: Monte Carlo (MC) sampling of MSMS-to-SMMS transitions via random metal-\nsulfur exchange\n• Epoch 3: MC sampling of alloyed SMMS intermediates through metal atom swapping\n• Epoch 4: MC sampling of a sulfur layer on SMMS via metal-S surface atom exchange\n• Epoch 5: Perturbed configurations near SMMS equilibrium states\n• Epoch 6: On-the-fly MD sampling of sulfur molecular dynamics (matching initial\nprotocol)\n• Epoch 7: S-S/Mo-S/Mo-W dimer configurations at varied distances\n• Epochs 8–17: MLPMD simulations of SMMS growth and Mo deposition (key structures\nextracted every 10,000 steps)\nIterative model refinement incorporated both historical and new data. We randomly sampled\n50,000 data points from newly generated configurations (with some structural repetitions being\ninevitable during sampling), which we defined as the new data. During each epoch, we constructed\na training set of 100,000 configurations through 1:1 random sampling of historical data\n(cumulative from prior epochs) and newly generated data. The MLP fine-tuning employed a\nlearning rate of 1×10−4 to balance convergence speed and accuracy. The final dataset contained\napproximately 26,000 configurations from initial training and all epochs, ensuring structural\ndiversity and representativeness critical for robust MLP development. The dataset has been\nuploaded to the Zenodo repository (https://doi.org/10.5281/zenodo.13777360).”\nComment 2: The MLP validation consists of fairly meaningless RMSE’s and some MD\nsimulations validating(?) the MLP and AIMD validating the stability predictions. These\ncalculations have limited certainty of the accuracy. What I believe is much more important is, e.g.\ncomparing the predicted relative stability of optimized and non-optimized structures (e.g. SMeMeS\nvs. MeS ) via the MLP vs. DFT.\nresponse to the comprehensiveness of MLP validation, we added the comparative analysis of the\nrelative stabilities of critical configurations. By systematically comparing the energy predictions\nof optimized and non-optimized structures obtained from MLP and DFT (see the details in Table\nS1 shown below), it can be seen that a high consistency between MLP and DFT in the relative\nstabilities of key configurations, such as SMMS and MS . Specifically, for the SMMS-like\n2\nconfiguration generated by direct molecular dynamics deposition and the optimized SMMS\nstructure, the energy differences are 0.080 eV and 0.068 eV, obtained from DFT calculations and\nMLP predictions, respectively. Similarly, for the SMoWS system, the corresponding energy\ndifferences are 0.090 eV from DFT simulation and 0.077 eV from MLP prediction. These results\ndemonstrate that our MLP can accurately compute the energy differences between different\nstructures and reliably describe the relative stabilities of materials.\nTo increase our manuscript's readability, we added these discussions in the revision (page 5)\nand a new table (S1) in the revised Supplementary Information (page 4):\n“By systematically comparing the energy predictions of key structures obtained from MLP and\nDFT shown in Table S1, it can be seen that a high consistency between MLP and DFT in the\nrelative stabilities of key configurations such as SMMS and MS . Specifically, for the SMMS-like\n2\nconfiguration generated by direct MD deposition and the constructed and optimized SMMS\nstructure, the energy difference is 0.080 eV and 0.068 eV, obtained from DFT calculations and\nMLP predictions, respectively. Similarly, for the SMoWS system, the corresponding energy\ndifference is 0.090 eV from DFT simulation and 0.077 eV from MLP prediction. These results\ndemonstrate that our MLP can accurately compute the energy differences between different\nstructures and reliably describe the relative stabilities of materials.”\n“Table S1. Comparison of energies (eV/atom) for critical structures calculated by DFT\nand MLP.\nStructure DFT (eV/atom) MLP (eV/atom) ΔE (eV/atom)\nS solid -4.323 -4.311 -0.013\n1H-MoS -7.582 -7.583 0.001\n2\n1T-MoS -7.298 -7.300 0.002\n2\n1T'-MoS -7.390 -7.380 -0.010\n2\n1H-WS -8.214 -8.144 -0.070\n2\n1T-WS -7.893 -7.899 0.005\n2\n1T'-WS -8.026 -8.017 -0.009\nMoSMoS -7.978 -8.007 0.029\nSMoMoS -8.409 -8.382 -0.027\nSMoWS -8.895 -8.892 -0.002\nSMMS -8.897 -8.895 -0.002\nSMoMoS (from MD) -8.329 -8.314 -0.015\nSMoWS (from MD) -8.804 -8.815 0.011\n”\nComment 3: The non-equilibrium MD simulations are also unclear. How exactly are new\nspecies introduced during the MD simulations? Why are they introduced with low velocities.\nFurthermore, why didn’t the authors collect statistics and clearly quantify e.g. mean times of\nincorporation. The efficient MLP would certainly allow for this. Currently, all results are\nqualitative observations without any quantitative indication.\nOur reply: Manythanks for the reviewer’s thorough review and insightful comments.\n(1) How exactly are new species introduced during the MD simulations?\nThe atomic deposition was achieved using the fix deposit command in LAMMPS, which\nintroduces a Mo atom with a downward velocity of 5 Å/ps at random positions within the\ntop region of the simulation box every 6000 steps. The initial kinetic energy of deposited\natoms was set to 0.12 eV. Regarding the implementation details of non-equilibrium MD\nsimulations, we have provided an expanded description in the Methods section (page 18):\n“The atomic deposition was achieved through the fix deposit command in LAMMPS, which\nintroduces a Mo atom with downward velocity of 5 Å/ps at random positions within the top region\nof the simulation box every 6000 steps. The initial kinetic energy of deposited atoms was set to\n0.12 eV.”\n(2) Why are they introduced with low velocities?\nIn non-equilibrium MD simulations, a lower kinetic energy was chosen to show the\nintrinsic nature of the embedding behavior, as a better representation of the embedding behavior\ndoes not necessitate a high kinetic energy. As shown in Fig. R1, we simulated cases with 1 eV\nand 5 eV of incident kinetic energy, exhibiting similar embedding behavior.\nFig. R1: (a) Simulated 1 eV kinetic energy deposition of a single Mo layer at 900 K on MoS .\n2\n(b) Simulated 5 eV kinetic energy deposition of a single Mo layer at 900 K on MoS .\nSupplementary Information:\n“As shown in Fig. S4, when Mo atoms were deposited at higher kinetic energies, such as 1 eV\nand 5 eV, they remain embedded in MoS . This indicates that the embedding behavior of Mo atoms\n2\nis not sensitive to kinetic energy.”\n(1) why didn’t the authors collect statistics and clearly quantify e.g. mean times of\nincorporation?\nWe thank the reviewer for the valuable comments. Quantitative data are essential for\nunderstanding the MD process. Therefore, in the revised manuscript, we systematically quantified\nthe evolution of metal-metal and metal-sulfur bond counts [Fig. R2, Fig. R3 and Fig. R4]. These\nquantitative data reveal the atomic incorporation process and sulfurization reaction kinetics.\nSpecifically, the progressive increase in metal-metal bonds indicates the continuous embedding of\ndeposited atoms. The sudden surge of the metal-sulfur bonds during sulfurization corresponds to\na critical stage of coordination reconstruction. These quantitative analyses provide crucial\nevidence for understanding the formation mechanism of intermediate structures.\nFig. R2: Number of Mo/W-Mo/W bonds (including Mo-Mo, Mo-W, and W-W bonds) during\nMo deposition on a WS substrate. A bond cutoff distance of 2.85 Å was applied for bond\n2\nidentification in both cases.\nFig. R3: Number of Mo/W-Mo/W bonds (including Mo-Mo, Mo-W, and W-W bonds) during\nMo deposition on a WS substrate. A bond cutoff distance of 2.85 Å was applied for bond\n2\nidentification in both cases.\nFig. R4: (a) Evolution of the number of Mo-S bonds during sulfur deposition on the SMoMoS\nintermediate structure. Excess sulfur atoms extract Mo atoms from the SMoMoS structure, leading\nto the formation of a bilayer MoS heterostructure. (c) Evolution of the number of Mo/W-S bonds\n2\nduring sulfur deposition on the alloyed S(Mo W )(W Mo )S structure. Sulfur atoms facilitate\n0.5 0.5 0.5 0.5\nsimultaneous extraction of both Mo and W atoms, ultimately forming an alloyed bilayer\nMo W S configuration. A cutoff distance of 2.6 Å was applied for both Mo-S and W-S bond\n0.5 0.5 2\nidentifications.\nThree new figures (Fig. S3, Fig. S8 and Fig. S12) were added in the Supplementary Information.\nRelevant discussions were added in the revised main text:\nOn Page 7,\n“Fig. S3 depicts the evolution of Mo-Mo bonds during continuous deposition, indirectly\nreflecting the embedding process. We observed that the number of Mo-Mo bonds gradually\nincreased with continued deposition and reached saturation once deposition ceased.”\nOn page 9,\n“Fig. S8 shows the evolution of Mo/W–Mo/W bonds during continuous Mo deposition,\nindirectly reflecting the degree of Mo atom insertion. The number of these bonds gradually\nincreases with deposition and saturates once deposition ceases.”\nOn page 12,\n“As shown in Fig. S12(a), there are very few initial Mo-S bonds because the S atoms deposited\non the SMMS surface cannot bond with Mo. As the MD simulation progresses, some Mo atoms\nare pulled to the surface to bond with the deposited S atoms, leading to an increase in the number\nof Mo-S bonds. Between 700 and 1000 ps, a large number of Mo atoms are pulled to the surface,\nresulting in a substantial increase in Mo-S bonds. The number of Mo-S bonds tends to be saturated\nafter about 1.2 ns, when a bilayer MoS is formed.”\n“Therefore, the resulting structure was not a distinct MoS /WS vdWHs but rather an alloyed\n2 2\nMo W S /Mo W S vdWHs. As shown in Fig. S12(b), the evolution of Mo/W-S bonds\n2 2\nclearly illustrates this process.”\n𝑥 1−𝑥 1−𝑥 𝑥\nComment 4: The discussion on of the growth mechanism should be more organized. Each\nparagraph appears as a loose description of what simulation was performed and what was seen\nlacking motivation that would help the reader gaining the very important scientific insight.\nOur reply: We thank the reviewer for the thorough review and insightful comments. In the\nrevised manuscript, we have systematically restructured the mechanistic analysis framework to\nenhance logical flow and scientific clarity. The discussion now progresses through a well-defined\nhierarchy of fundamental processes: First, we employ molecular dynamics simulations to reveal\nthe intrinsic tendency of metal atoms to spontaneously embed beneath the sulfur layer during initial\ndeposition, forming subsurface metal-bonded configurations. This is followed by comprehensive\nthermodynamic calculations and phonon dispersion analysis that establish the SMMS intermediate\nas a structurally stable transition state. Next, simulations of hetero-deposition processes elucidate\nthe thermodynamic preference and stabilization mechanisms for alloyed SMMS formation. Non-\nequilibrium MD simulations then map the dynamic evolution pathways of these intermediates\nduring sulfurization, while comparative deposition studies demonstrate how co-deposition of Mo-\nS clusters effectively suppresses SMMS formation. Finally, electronic structure calculations unveil\nthe unique advantages of these intermediates for device applications. To improve navigability, we\nhave introduced descriptive subheadings to guide readers through each key mechanistic stage.\nComment 5: In illustrate the latter point, data from various simulations appears sometimes\nerratically combined. An example is the discussion of the stability of the bilayer system. If I\nunderstood correctly, the structure for MoSMoS, SMoMoS, SMo S, SMo S (Fig2. (f-i)) are\n3 4\nobtained by geometry optimization, not by MD simulation. Did the authors observe the SMoMoS\n(Fig.2 (g)) structure during the MD simulation starting from the MoSMoS structure (Fig.2 (f))?\nIf not, the statement of “all deposited Mo atoms spontaneously sank into the MoS layer, tending to\n2\nform the intriguing SMoMoS structure, as shown in Fig. 2(g).” and “In addition to the SMoMoS\nstructure formed by depositing 1 ML of Mo atoms” might be misleading, because as seen in Fig.\nS5 (b) it is less likely that the formation of SMoMoS structure from MoSMoS structure happens.\nAlso, in this discussion, the authors showed a lot of data:\n○ energy landscape of MoSMoS, SMoMoS, SMo S, SMo S structures (Fig2. (f-i))\n3 4\n○ phonon band structures (Fig2. (l))\n○ tendency that Mo atoms sink beneath the S layer during the MLP-MD simulation (Fig2. (e))\n○ tendency that Mo atoms sink beneath the S layer during the AIMD simulation (Fig S5(b))\n○ stability of SMoMoS during the MLP-MD simulation (Fig S3 (a)(d))\n○ stability of SMoMoS during the AIMD simulation (Fig S5 (c))\nAs such, the data from geometry optimization and MD simulations using DFT and MLP are\ngiven, and it is hard to understand how they complement each other in the current version of the\nmanuscript. This should be more clarified for better discussions.\nOur reply: We sincerely appreciate the reviewer's thorough critique. The original description\nregarding the formation mechanism of the SMoMoS structure indeed requires clarification. The\nfour configurations shown in Fig. 2(f-i) were obtained through DFT geometry optimizations, rather\nthan directly from MD simulation trajectories. The key evidence for the SMoMoS formation\npathway originates from MLPMD simulations, which reveal that Mo atoms spontaneously migrate\ninto the MoS₂ layer when directly deposited on its surface. MLPMD simulations of sequential Mo\ndeposition suggest that SMoMoS serves as a naturally formed SMMS (S-M-M-S) intermediate\nduring metal atom deposition. Crucially, the MoSMoS configuration was never observed as a\ntransition state in simulations. Formation energy calculations further indicate its high-energy state,\nconfirming negligible likelihood of existing as a metastable phase. We acknowledge that the\nstatement of “all deposited Mo atoms spontaneously sank into the MoS layer, tending to form the\n2\nintriguing SMoMoS structure, as shown in Fig. 2(g).” and “In addition to the SMoMoS\nstructure formed by depositing 1 ML of Mo atoms” might be misleading. We have revised the\nrelevant discussion in the manuscript. Throughout the MLPMD simulation, no Mo atoms formed the\nMoSMoS structure [Fig. R5(f)]. Instead, all deposited Mo atoms spontaneously sank into the\nMoS layer. These findings suggest that under practical deposition conditions, Mo atoms\n2\npreferentially form an SMoMoS structure [Fig. R5(g)]. This SMoMoS structure was then\nconstructed and optimized by DFT. The results show this SMoMoS structure exhibits an energy\n0.080 eV/atom lower than that of the directly formed embedded structure from MD simulations\n[Fig. R5(e)].\nTo reveal the dynamic behavior during the actual growth process, the MLP-MD simulations\n(Fig. R5(e)) demonstrate the deposition dynamics of metal atoms and the pathways of structural\nevolution on a nanosecond timescale. Motivated by these simulation insights, we systematically\nconstructed and optimized the structures (Fig. R5(f-i)). Fig. R5(j) shows the relative\nthermodynamic stability of these structures. Additionally, the phonon spectrum analysis (Fig.\nR5(l)) further validates the dynamical stability of SMoMoS, ruling out the possibility of imaginary\nfrequencies. Long-term MLPMD simulations (5.5 ns) further validate the structural robustness of\nSMoMoS under thermal fluctuations. Furthermore, a shorter period of AIMD simulation was\nconducted to verify the reliability of the embedding behavior and the stability of the SMMS\nstructure, serving as a validation and complement to the MLP results.\nFig. R5: (a) The structure of Mo atoms deposited and (b) embedded into the MoS layer, (c)\n2\nalong with the corresponding energy plots (blue line represents raw energy, while red line shows\nthe low-pass-filtered result). (d-e) Snapshots of 0.25 monolayer (ML) (d) and 1.0 ML (e) Mo\natom deposition simulated at 900 K on the existing MoS layer. (f-i) Four possible structures:\n2\n(f) MoSMoS, (g) SMoMoS, (h) SMo3S, and (i) SMo4S with more Mo atoms embedded. (j)\nFormation energy convex hull of the considered MoS structures, with green dots corresponding\nto energies of structures found in the 2D Materials Database (MatHub-2d). (k) The unit cell\nstructure of SMoMoS and (l) its phonon dispersion relation.\nWe have modified these discussions in the revised main text (page 7):\n“\nThroughout the MLPMD simulation, no Mo atoms formed the MoSMoS structure [Fig. 2(f)].\nInstead, all deposited Mo atoms spontaneously sank into the MoS layer. These findings suggest\n2\nthat under practical deposition conditions, Mo atoms preferentially form an SMoMoS structure\n[Fig. 2(g)]. This SMoMoS structure was then constructed and optimized using DFT, revealing that\nit exhibits an energy 0.080 eV/atom lower than that of the directly formed embedded structure\nfrom MD simulations [Fig. 2(e)].\nWe conducted a detailed comparison of these two structures in terms of energy and dynamic\nstability. The energy significantly decreased by 1.73 eV per Mo atom when transitioning from\nMoSMoS to SMoMoS, indicating a strong driving force towards the formation of the SMoMoS\nstructure. In addition to the SMoMoS structure formed by depositing 1 ML of Mo atoms, we\nconstructed and optimized structures such as SMo S [Fig. 2(h)] and SMo S [Fig. 2(i)]. ], which\n3 4\nmay form with an increased number of Mo atoms. The convex hull of MoS compounds was\nplotted by varying the Mo ratio, referencing the bulk phases of elements Mo and S, as shown in\nFig. 2(j). Notably, MoS has the lowest energy, while SMoMoS is situated 196 meV/atom above\n2\nthe convex hull. We also found structures that match this elemental ratio from other 2D materials\ndatabases (MatHub-2d)37[Fig. S5], with their formation energies indicated by green dots in the\nfigure, showing that they are 253 meV/atom higher than the SMoMoS structure. SMo S and\n4\nwhen approximately 1 ML of Mo atoms is deposited on the MoS surface, SMoMoS is the most\n2\nlikely intermediate. The phonon dispersion of the SMoMoS structure showed no imaginary\nfrequencies [Fig. 2(k) and Fig. 2(l)], further confirming that MoSMoS would spontaneously\ntransform into SMoMoS. It is worth noting that the phonon dispersion simulated by the MLP\naccurately reproduced the DFT results, validating the accuracy of our MLP in describing structural\nstability. Long-term MLPMD simulations further validated the stability of SMoMoS, as shown in\nFig. S6.\nFurthermore, we used DFT to validate the key processes described above. Fig. S7(a) illustrates\nthe process of a single Mo atom into the MoS layer. Despite the relatively low initial kinetic\n2\nenergy, the Mo atom quickly embedded into the MoS layer within approximately 700 fs, as\n2\nindicated by the sudden drop in energy. Notably, compared to pre-adsorbed Mo atoms, vacuum-\ndeposited Mo atoms significantly shorten the embedding time to about 700 fs — much shorter than\nthe approximately 50 ps required for surface-adsorbed atoms in previous simulations. Fig. S7(b)\nexpands our investigation of the behavior of the Mo atomic layer on MoS . As the Mo atoms\n2\ngradually embed into the MoS structure, they descend accordingly. This energy decrease is\n2\nsignificantly lower than that of the initial configuration, further confirming the thermodynamic\nfavorability of Mo atom embedding. Fig. S7(c) focuses on the stability of the formed SMoMoS\nstructure over approximately 10 ps.\n”\nComment 6:\nMinor comments:\n● Please check the numbering for the figures. For example, Fig. S4 comes earlier than Fig.\nS3, which is not kind to readers.\nOur reply: Thank you for your careful review of our manuscript and the valuable comments.\nWe have repositioned all supplementary figures to ensure their numbering is consistent with the\nreferences in the main text.\nComment 7: ● Small comments on the figures:\n○ Fig.2 (j) it is stated that “with green dots” but I couldn't find the green point\n○ Fig.2 (c) what is the red line and blue line?\n○ For the caption for Fig.2, the labels (a) to (e) come after the caption but (f) to (l) come before\nthe caption, which is, I think, very confusing.\n○ Fig.5 (b) SMOWS should be SMoWS\nOur reply: Thank you for your detailed and constructive feedback, which has significantly\nenhanced the clarity and rigor of our manuscript. We have implemented systematic revisions to\naddress the issues related to the figures and ensure the accuracy of data presentation. We have\nredrawn the convex hull plot to better highlight the green dots [Fig. R8]. For Fig. 2(c), we added\nclarifications in the caption: the blue line represents the original energy data, while the red line\ndepicts the trajectory after low-pass filtering. This filtering process effectively suppresses high-\nfrequency thermal noise, providing a more intuitive representation of the energy transitions during\nthe atomic deposition process.\nFig. R6: Revised Fig. 2(j).\nComment 8: ● For the phonon band structures, the data from DFT and MLP are plotted\ntogether on the same figure, but I suggest to show only DFT results. The authors stated in page 3\n“It is worth noting that the phonon dispersion simulated by the MLP accurately reproduced the\nDFT results, validating the accuracy of our MLP in terms of atomic forces.”, but some degeneracy\nat the gamma point are not well described by MLP, and it is hard to tell the accuracy of MLP from\nthis figure. For Fig.3 (j-l), the authors don’t mention the differences between DFT and MLP, so I\nsuggest removing the phonon band structures by MLP for better clarity.\nOur reply: We sincerely thank the reviewer for the valuable suggestion. The MLP developed\nin this study is primarily aimed at dynamic simulations of growth processes. Although there are\nslight differences in the description of phonon degeneracy near the Γ point by the MLP, this does\nnot affect its reliability in simulating non-equilibrium processes such as atomic deposition and\nstructural reorganization. Furthermore, retaining the MLP results can demonstrate the reliability\nof the MLP in describing structural stability. Therefore, we believe that the MLP portion should\nbe retained.\nComment 9: ● Equations in Page 4\n○ Please define N and T.\n○ E is energy, then why is E formation energy?\nSMoWS\n○ x should be “one Mo atom” not “one metal atom”?\nOur reply: We sincerely thank the reviewer for the valuable suggestion. In this equation, N\nrepresents the total number of atoms in the system, T is the thermodynamic temperature, E\nSMoWS\ndenotes the energy of the unalloyed SMMS structure, E represents the energy of the alloyed\nconfiguration, and x is the fraction of Mo atoms in the upper metal layer.\nWe have added these descriptions in the revised main text (page 10):\n“To determine the thermodynamic driving force for alloying, we calculated the free energy\nchange for varying distributions of metal atoms between the upper and lower layers. Here,\nis defined as\n𝛥𝐹 𝛥𝐹\n1\nwhere denotes the ener𝛥g𝐹y o=f the[ 𝐸un−al𝐸lo S y M e o d W S S M−M𝑇S𝛥 𝑆st]r,ucture [Fig. 3(h)], is the energy\n𝑁\nof the alloyed configuration, is the temperature (300 K), and represents the mixing entropy\nfor Mo/W𝐸 S a M llo oW yi S ng. The mixing entropy is given by 𝐸\n𝑇 𝛥𝑆\n𝑆\nwhere is the Boltz\n𝑆\nm\n=\nan\n(\nn\n𝑁\ncMoon\n+\nsta\n𝑁\nnWt,\n)𝑘\ni𝑏s\n[\nt\n−\nh\n𝑥\ne\nl\nf\nra\n(𝑥\nct\n)\nio\n−\nn\n(\n1\no\n−\n𝑓 𝑥 M\n−\no\n𝑥 )\n𝑎\nl\n𝑡\n𝑛\no\n𝑚\n𝑠\n𝑖𝑛\n𝑡\nℎ\n]\n,\nthe upper metal layer, and\nand are the numbers of Mo and W atoms, respectively.\n𝑘𝑏 𝑥\n𝑁M”o 𝑁W\nComment 10: ● Page 5: “Throughout the calculations, general precision was used”, it would\nbe better to specify the precision. “general precision” is vague.\nOur reply: We appreciate the reviewer's valuable feedback. The phrase \"general precision\" in\nthe original text was indeed misleading and has now been removed. The details of the precision\ncan be found in the Methods part.\nComment 11: ● For the simulation given in video S4, the S layer looks floating above the\nSMoMoS layer. Is it constrained?\nOur reply: In the referenced simulation (Video S4), no positional constraints were applied to\nany atomic layers, including the sulfur layer. The molecular dynamics simulations were conducted\nunder fully unconstrained NVT ensemble conditions, allowing all atoms to move freely according\nto the interatomic forces described by MLP.\nComment 12: ● Page 5, what do you mean by “a specific rate of 516 ps per layer”?\nOur reply: We thank the reviewer for their valuable comment. The phrase has been revised to\nclarify that the deposition of the second layer requires a duration of 516 ps.\nWe have this description in the revised main text (page 14):\n“The Mo and S atoms required for the growth of one layer of MoS were deposited onto the\n2\nsurface within 516 ps, after which the deposition was stopped.”\nComment 13: ● For most of the systems, the system size is not clearly stated but it should be\ngiven.\nOur reply: We sincerely thank the reviewer for the constructive suggestion. In response, we\nhave added detailed system size parameters in the \"MD\" section of the Methods. Specifically, we\nhave included the dimensions of the supercell used for Mo atom deposition simulations,\nthe system sizes for SMoMoS and SMMS sulfurized process simulations, and the diffusion\nbehavior of Mo-S clusters on the MoS4√ s3ur×fa7ce. The supercell employed for studying the\n2\ngrowth of the second MoS₂ layer.\n6√3×8\nWe have this description in the revised main text (page 18):\n“For Mo atom deposition, we used a MoS or WS supercell at 900 K. The SMoMoS\n2 2\nand SMMS sulfurized process simulations were performed using a supercell at 1100 K.\n4√3×7\nThe SMoMoS and SMMS sulfurized process simulations were performed using a\nsupercell at 1100 K. The diffusion behavior of Mo-S clusters on4 √th3e ×M7oS surface was also\n2\ninvestigated by MoS supercell at 1100 K. For the simultaneous Mo and S depo4s√iti3on× t7o\n2\ngrow the second MoS layer, we employed a larger MoS or WS supercell at 1100 K.”\n2 2 2\n4√3×7\n6√3×8\nComment 14: ● Fig. S5(c) shows the AIMD simulation of SMoMoS structure until 10 ps, but\nwhat happens after 10 ps? Did the system retain the SMoMoS structure?\nOur reply: We thank the reviewer for raising this important question. Yes, the system retains\nthe SMoMoS structure after 10 ps. Due to computational resource limitations, our AIMD\nsimulation was indeed restricted to a short duration of 10 ps. However, during this simulation, no\nsignificant structural degradation or instability was observed in the SMoMoS configuration. To\ncomprehensively assess its thermodynamic stability, we conducted an extended simulation using\nMLPMD for a duration of 5.5 ns. The results demonstrate that the intermediate structure maintains\nits intact lattice configuration throughout the simulation, further confirming its thermodynamic\nstability as a growth intermediate.\nComment 15:● Again, are the structures of Fig.3 (g-i) obtained by geometry optimization and\nnot appeared during the MD simulations?\nOur reply: We appreciate the reviewer's careful observation. The structures shown in Fig. 3(g-\ni) were indeed constructed through geometry optimization as theoretical models and did not appear\nduring the MD simulations. To eliminate any confusion between these theoretical models and the\ndynamic evolution results observed in the MD simulations, we have added clarifying remarks in\nthe main text.\nWe added this description in the revised main text (page 10):\n“Due to the similar interaction between Mo atoms and WS , as previously observed with MoS ,\n2 2\nthe constructed SMoWS intermediate structure shown in Fig. 3(h) emerges as another crucial\nconfiguration to consider. Its energy is 0.090 eV/atom lower than that of the directly formed\nstructure from MD simulations [Fig. 3(c)].”"}, {"round": 2, "reviewer_comments": {"reviewer_1": "In this revision, Zhao and coworkers largely addressed inquiries and comments, as well as attempted to make appropriate changes to the manuscript. While many new details were conveyed in the revision, raised concerns about reproducibility, quality of the MLP, quality of the simulations were not satisfactorily rectified. Furthermore, with the new understanding, new questions emerged about the validity of the simulations in relation to the experimental insights (Nature 621, 499 15 (2023)), they are investigating. These concerns are paired with a remaining lack of readability/clarity of the manuscript, that despite being commented by other reviewers, has not been properly addressed. It would be advisable to put more care into the written text and scientific reporting which still appears to convey numerous mistakes. Compared to the first version, where the benefit of doubt about unknowns prevailed, their work appears now less rigor, raising doubt whether it is publishable.\n\nComments:\n- Details on MLP training procedure which the authors added are still vague. For the final data set, the authors should at least explain the model size (e.g. number of atoms & unit cells) and the number of configurations in the training and testing set (to assess the actual number of data-points = energies & forces). The manuscript should be self-contained and it can’t be expected of the reader to download the Zenodo repo to answer these questions. An exact understanding of the training set is critical, since by now it becomes clear that the training process was manually and not systematically performed.\n- The various methods to create training configurations need much more detailed explanation. For example, using VASP’s MLP-MD feature to get data points raises the questions, whether the data that went to training the Nequip potential were DFT calculations or VASP’s MLIP prediction (which would introduce considerable noise). Also, the explanation of the Monte Carlo method implemented by the authors remains unclear and should be described in more detail. Quite confusing is the also the wording and usage of words like “epoch” that appears out-of-place and mixes with terminology of NN training (e.g. as used in Nequip). The process of drawing 100k samples from a much larger structure pool, leading to excessive redundancy of 74% also raises questions of good-practice.\n- The MLP validation still falls short. From the revision it becomes apparent that no uncertainty metric was used during training (e.g. to track convergence), and instead qualitative criteria “unstable molecular dynamics” or “poor phonon spectra” called for “refining”. This should be clearly labelled as an ad-hoc approach and validation of the MLP needs to be especially meticulous. The comparison of “critical structures” is a good first step, but by far not enough to determine the MLPs accuracy over the investigated composition and structure space. A sign of significant discrepancy between MLIP and DFT is Fig. 3(j), and certainly raises doubt over the statement on page 8 “accurately reproduced DFT results”.\n- The authors did not implement the previous inquiry about systematic simulations and statistics. While a more thorough data analysis is appreciated, also more test simulations (starting from different initial conditions) should be performed – which is easily achievable with the highly performant Nequip potential. But also, in the newly provided analysis, many questions arise: Why do the number of bonds start at zero in Fig. S3 and S8? What about existing bonds? Those seem to be accounted for in Fig. S12. Note also, S3 and S8 have the same caption, which is probably another mistake. Further, it might make sense to distinguish the bonds between different elements, which could help understand the growth process.\n- The revised description of the simulations still requires careful refinement. It is not, yet, fully clear what is done and in the current description, terminology of ensembles is mixed up.\n- The discussion of the manuscript results still requires better structure, a lot of back-and-forth is contained around Fig. 3 making it hard to follow the authors argumentation. After revision, new questions about the simulations occurred. For example, a key result is the fact, that “metal films do not float” (compared to the original experimental data in Nature 621, 499 15 (2023)), what is the significance of this finding?\n- The results also contain quite a bit of inconsistent simulations: The MD simulations are conducted around deposition temperatures of 500K-1300K (?), but some analysis (e.g. mixing entropy), and other simulations (e.g. depicted in Fig. R7) are conducted at 300K. This seems inappropriate if the synthesis mechanism is investigated. Some arguments are made about alloying where WS2 and MoS2 slabs appear not systematically exchanged and kinetic stabilization is argued on long MD timescales (“stability during 1.4 ns”) but without a significant number of samples. These simulations thus only weakly support the authors’ conclusions.", "reviewer_2": "Report on the Revised Manuscript by Zhao et al.\n1. The language and overall presentation of the manuscript have been significantly improved, making the paper much more readable and understandable. However, the abstract and the first two paragraphs of the introduction still require further polishing for clarity and conciseness.\n2. The authors now clarify that they used an ideal Morse potential for the substrate simulation. Given that the study is motivated by experimental results on sapphire, it is essential that the manuscript discusses how well this idealized simulation approach approximates the actual experimental conditions. A brief justification or discussion in the main text would strengthen the manuscript.\n3. The authors have added two supplementary figures to clarify the heterostructure formation. Figure S10 includes the simulation temperature, but this detail is missing in Figure S11. For clarity, the temperature should also be included in Figure S11 to clearly indicate the conditions under which the sulfurization process was studied.\nOnce these remaining issues are addressed, I believe the manuscript will be suitable for publication in Nature Communications.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "The authors have adequately responded to my comments and I am happy to recommend the manuscript for publication.\nAs a small comment, I noticed that in revised Fig. S14, the figure and caption refer to a \"2H\" phase, but in the rest of the manuscript, this phase is referred to as \"1H\". I think the authors probably meant to write \"1H\" here.\n\n(Remarks on code availability)\nThe code contained in \"ocp_active.zip\" is unlikely to be a usable resource for the community, as it appears to be available only in the zenodo repository as a set of python files (as opposed to a maintained repository on github), and it is unclear what files from the original Open Catalyst repository were modified or how they were modified. I don't think this should necessarily preclude publication, however, as the authors describe their methods in the manuscript and do make the code available, and more importantly they also provide their training and test structures in extxyz format, which may very well be useful for the community."}, "author_response": "Author reply: We thank the reviewer for careful reading of our revised manuscript and response to reviewers’ comments. We also thank the reviewer for raising more questions about our revised manuscript, which helped us to improve the research greatly. After further study and careful revision of the manuscript, we hope the reviewer finds the new version of the manuscript acceptable. In the following part, we answered the questions and comments by the reviewer point-by-point.\n\nComment 1: Details on MLP training procedure which the authors added are still vague. For the final data set, the authors should at least explain the model size (e.g. number of atoms & unit cells) and the number of configurations in the training and testing set (to assess the actual number of data-points = energies & forces). The manuscript should be self-contained and it can’t be expected of the reader to download the Zenodo repo to answer these questions. An exact understanding of the training set is critical, since by now it becomes clear that the training process was manually and not systematically performed.\n\nAuthor reply: We thank the reviewer for the suggestion on improving the clarity of our machine-learning potential training dataset description. We have carefully retrained our machine-learning potentials by including more systematic data. The corresponding SI-data.extxyz files of our dataset are provided. The dataset comprises 1.18×107 data points (26,434 energy points and 1.18×107 force points (3.9×106 atoms)), covering all key structures encountered during transition-metal dichalcogenide growth. Detailed statistics of our dataset are shown in Supplementary Table S2 and SI-data.extxyz files.\n\nIn the revised manuscript, we have added the following table, on Supplementary Information page 6:\nTable S2: Dataset statistics (including cell volumes)\nStructure Category Configurations Avg. Atom Count Typical Cell Size (Å)\nTMD surface with W-Mo-S 18396 155 22 × 22 × 30\ncluster interaction\nSingle-layer Mo on WS 668 196 22 × 22 × 21\n2\nsurface\nDisordered alloyed TMD 1,231 144 22 × 22 × 23\nstructures\nStructure Category Configurations Avg. Atom Count Typical Cell Size (Å)\n1T-phase TMD 1804 144 22 × 22 × 23\nFurther Sulfidized 360 280 21 × 21 × 21\nStructures of SMMS\nTMD nanotubes 403 62 32 × 32 × 5\nDisordered SMMS 900 196 22 × 22 × 21\nstructures\nOrdered SMMS structures 748 144 18 × 16 × 30\nMoSWS structures 96 143 19 × 16 × 22\nPristine mono-/bilayer 138 58 13 × 10 × 26\nMoS /WS and\n2 2\nheterostructures\nElemental sulfur 1101 35 14 × 14 × 19\nDimers 589 2 25 × 25 × 25\nTotal 26434\n\nComment 2: The various methods to create training configurations need much more detailed explanation. For example, using VASP’s MLP-MD feature to get data points raises the questions, whether the data that went to training the Nequip potential were DFT calculations or VASP’s MLP prediction (which would introduce considerable noise). Also, the explanation of the Monte Carlo method implemented by the authors remains unclear and should be described in more detail. Quite confusing is the also the wording and usage of words like “epoch” that appears out-of-place and mixes with terminology of NN training (e.g. as used in Nequip). The process of drawing 100k samples from a much larger structure pool, leading to excessive redundancy of 74% also raises questions of good-practice.\n\nAuthor reply: We thank the reviewer for the insightful comments on our training-set construction. We have revised the manuscript to provide a more detailed descriptions of how the training configurations were created. We address the technical details point by point as follows:\n1. Use of the VASP on-the-fly MD feature\nWe confirm that all data used to train the NequIP potential originate from full self-consistent DFT (SCF) calculations; no interpolated data predicted by any ML potential were included. In the VASP on-the-fly MD framework, if the predicted uncertainty exceeds the threshold value (ML_CTIFOR), a full DFT SCF calculation is automatically triggered. The ionic steps containing electronic self-consistency information are then extracted from vasprun.xml and added to the training set.\n2. Monte Carlo sampling method\nWe agree that the previous explanation of MC sampling method was insufficient and have expanded in the revised manuscript. Our method involves random atomic exchanges between Mo and W atoms across different layers, starting from a relaxed S-Mo-W-S configuration. Each exchange generates a new alloyed structure, and all generated configurations are accepted. While this approach does not follow the Metropolis algorithm, it is commonly referred to as a generalized Monte Carlo method in materials modelling, where the goal is to efficiently explore configuration space and enhance structural diversity. To avoid confusion, we now describe it as a “randomized atom-exchange sampling scheme” in the revised manuscript.\n3. Terminology normalization\nWe appreciate the reviewer’s comment regarding the confusion caused by the term “epoch”. To reduce ambiguity, all occurrences of “epoch” have been replaced by “training round” throughout the manuscript.\n4. Clarification of the 100k sampling strategy\nThe alleged “74% redundancy” does not occur in our workflow. In each training round, we randomly draw historical and newly generated data with a 1:1 ratio to fine-tune the MLP model. This balance prevents catastrophic forgetting while continuously introducing new chemical environments. The use of “100k” is simply a practical implementation detail of this sampling strategy. To avoid misunderstanding, we have revised the manuscript to clarify this point.\nWe also modified this in Supplementary Information on page 1:\n“In this study, we adopted a basic iterative training procedure of MLP that begins with generating an initial dataset from quantum mechanical simulations. Next, we validate the reliability of the model by comparing the data with DFT and try to apply it to simulate the 2-step growth of MoS . If the MLP results in unstable molecular dynamics and poor phonon spectral correspondence, we then fine-tune the MLP by adding more relevant DFT constitutive data. These data expansions and models retraining processes are repeated iteratively until the model achieves satisfactory accuracy. The final model can then be deployed for efficient and accurate large-scale atomistic simulations.\nThen, we employed on-the-fly molecular dynamics (MD) simulations to generate the initial training set, effectively exploring unknown configuration spaces while minimizing sampling costs. Using VASP’s on-the-fly MD framework, we implemented the following parameters to optimize efficiency: a timestep of 3.5 fs ensured accurate atomic motion sampling, 20,000 maximized simulation steps, fixed unit-cell shapes, Nosé-Hoover thermostats, and a temperature range (500 K to 1300 K) that covered typical growth conditions. The ML_MB=2000 parameter constrained the number of local reference configurations in the on-the-fly MLP. Initial configurations comprised monolayer MoS and WS and their combinations with homogeneous/heterogeneous junctions and adsorbed metal clusters. We expanded the training set through iterative training rounds with distinct sampling strategies:\n• Training round 1: On-the-fly MD sampling of MSMS structure dynamics\n• Training round 2: Randomized metal-sulfur exchange sampling — Starting from the relaxed MSMS structure, we performed random substitutions between metal (Mo/W) and sulfur atoms at specific interfacial sites to generate intermediate SMMS-like configurations. These exchanges emulate early-stage metal–chalcogen mixing observed during the deposition process.\n• Training round 3: Randomized metal-metal swapping — To construct a diverse set of alloyed SMMS structures, we performed layer-constrained Mo–W swaps across top and bottom layers. This preserves the overall composition but introduces configurational disorder relevant to alloying behaviour.\n• Training round 4: Surface metal–sulfur exchange sampling — Starting from SMMS structures, we randomly exchanged surface S atoms with underlying metal atoms (Mo or W) to simulate possible sulfurization pathways and local rearrangements during the growth process. atom exchange\n• Training round 5: Perturbed configurations near SMMS equilibrium states\n• Training round 6: On-the-fly MD sampling of sulfur molecular dynamics (matching initial protocol)\n• Training round 7: S-S/Mo-S/Mo-W dimer configurations at varied distances\n• Training rounds 8-17: MLP-MD simulations of SMMS growth and Mo deposition (key structures extracted every 10,000 steps).\nDuring each Training round, we fine-tune the MLP model through 1:1 random sampling of historical data (cumulative from prior training rounds) and newly generated data. The MLP fine-tuning employed a learning rate of 1×10−4 to balance convergence speed and accuracy. The final dataset contained approximately 26,000 configurations from initial training and all training rounds, ensuring structural diversity and representativeness critical for robust MLP development. The dataset has been uploaded to the Zenodo repository (https://doi.org/10.5281/zenodo.13777360) and included in SI-data.extxyz”\n\nComment 3: The MLP validation still falls short. From the revision it becomes apparent that no uncertainty metric was used during training (e.g. to track convergence), and instead qualitative criteria “unstable molecular dynamics” or “poor phonon spectra” called for “refining”. This should be clearly labelled as an ad-hoc approach and validation of the MLP needs to be especially meticulous. The comparison of “critical structures” is a good first step, but by far not enough to determine the MLPs accuracy over the investigated composition and structure space. A sign of significant discrepancy between MLIP and DFT is Fig. 3(j), and certainly raises doubt over the statement on page 8 “accurately reproduced DFT results”.\n\nAuthor reply: We thank the reviewer for raising this critical point regarding the MLP validation. The configuration shown in Fig. 3(j) is thermodynamically unstable due to its relatively high energy, leading to insufficient accuracy of the MLP force constants for this structure. We have carried out targeted re-sampling of this configuration and incorporated perturbed variants into the training set. Specifically, we added 668 configurations obtained from VASP on-the-fly AIMD sampling of a single-layer Mo on WS surface, as well as 96 perturbed configurations of the MoSWS structure, to enrich the relevant local configuration space. After retraining the MLP with the expanded dataset, the agreement between the MLP and DFT phonon spectra has improved substantially. Although minor frequency deviations remain, the overall dispersion now reproduces the DFT result reasonably well (see new Fig. R1).\nWe also replaced the Fig. 3(j) accordingly.\nFigure R1: Phonon dispersion of the MoSWS.\n\nComment 4: The authors did not implement the previous inquiry about systematic simulations and statistics. While a more thorough data analysis is appreciated, also more test simulations (starting from different initial conditions) should be performed – which is easily achievable with the highly performant Nequip potential. But also, in the newly provided analysis, many questions arise: Why do the number of bonds start at zero in Fig. S3 and S8? What about existing bonds? Those seem to be accounted for in Fig. S12. Note also, S3 and S8 have the same caption, which is probably another mistake. Further, it might make sense to distinguish the bonds between different elements, which could help understand the growth process.\n\nAuthor reply: We thank the reviewer for further raising the issue of testing simulations. We have performed more simulations to verify the methodology and the conclusions presented in the manuscript.\n1. More test simulations from multiple initial conditions: We thank the reviewer for the suggestion regarding statistical repetition. After careful consideration, we believe that, for the embedding process investigated in this work, performing multiple simulations from randomized initial velocity distributions is not scientifically necessary. First, from a theoretical standpoint, DFT energy calculations indicate that the embedding of an Mo atom in MoS is strongly thermodynamically favoured, with an energy release of ~1.45 eV. Second, in the present study, we have repeatedly observed consistent embedding behaviour across a variety of independent simulations, including different substrates (MoS and WS ), temperatures (300–900 K), and incident kinetic energies. While the initial velocities in these simulations were not explicitly constructed for statistical sampling, they effectively serve as physically distinct perturbations, under which the same qualitative behaviour was observed. This indicates that the embedding process is robust and not sensitive to the specific initial conditions. Based on these considerations, we believe that additional randomized-velocity simulations, while potentially helpful, would not change the physical conclusions we draw regarding the embedding mechanism. We nonetheless acknowledge the general value of such statistical tests and plan to explore them further in future studies involving more complex dynamical pathways.\n2. Clarification on Figs. S3, S8, and S12: Regarding the issue of zero points in bond counts, we provide a complete physical explanation: the adopted cutoff radius of 2.85 Å is larger than typical M-M bonds in M bulk (2.74 Å and 2.75 Å for Mo and W, respectively) but smaller than the metal-metal equilibrium spacing in MoS /WS (3.15 Å). The essential difference between Fig. S12 and S3/S8 lies in the different observation objects: S12 records the change in Mo/W-S bond counts during the sulfurization process, and its initial value >0 reflects the existence of inherent metal-sulfur bonds in the substrate. While S3/S8 records the number of newly formed metal bonds during the deposition process, which does not include the initial bonding state of the substrate.\n3. Distinguishing Bonds Between Different Elements: We thank the reviewer for suggesting the differentiation of metal–metal bonds (Mo–Mo, Mo–W, and W–W). We have added Fig. R2, which presents the time evolution of different types of metal–metal bonds during the deposition of Mo atoms onto a WS substrate. A cutoff distance of 2.85 Å was used for bond identification. As shown in Fig. R2, the number of Mo–Mo and Mo–W bonds increase rapidly during the early stages of deposition, reflecting the aggregation of Mo atoms and their bonding interactions with the underlying W atoms. In contrast, the number of W–W bonds increase more slowly. This increase does not represent the formation of new chemical bonds but rather results from local lattice rearrangements during the deposition process. Specifically, some W–W pairs that originally had interatomic distances greater than the cutoff (the W–W distance in pristine WS is approximately 3.15 Å) experience slight displacements or strain-induced contractions, causing their separation to fall below the 2.85 Å threshold and thus be counted as bonds in the analysis. Therefore, the gradual rise in W–W bond count should be interpreted as a reflection of subtle structural reconstructions in the static lattice induced by Mo atom deposition, rather than the creation of new W–W bonds.\nFigure R2: Number of Mo-Mo, Mo-W, and W-W bonds during Mo deposition on a WS2 substrate. A bond cut-off distance of 2.85 Å was applied for bond identification in both cases.\nIn the revised manuscript, we have modified these captions in Supplementary Information, on page 8:\n“FIG. S3: Number of Mo-Mo bonds during Mo deposition on a MoS substrate. Bond identification is based on a cut-off distance of 2.85 Å, which is larger than the Mo-Mo bond length in Mo bulk (~2.74 Å) but smaller than the equilibrium Mo–Mo spacing in pristine MoS (~3.15 Å).”\nOn page 12:\n“FIG. S8: (a) Evolution of the total number of Mo/W–Mo/W chemical bonds (including Mo–Mo, Mo–W, and W–W bonds) during Mo deposition on a WS substrate. (b) Individual counts of Mo–Mo, Mo–W, and W–W bonds. Bond identification is based on a cut-off distance of 2.85 Å, , which is larger than the Mo-Mo bond length in Mo bulk (~2.74 Å) and W bulk (~2.75 Å) but smaller than the equilibrium Mo–Mo spacing in pristine MoS (~3.15 Å)”\nOn page 15:\n“FIG. S12: (a) Evolution of the number of Mo-S bonds during sulfur deposition on the SMoMoS intermediate structure. Excess sulfur atoms extract Mo atoms from the SMoMoS structure, leading to the formation of a bilayer MoS heterostructure. (b) Evolution of the number of Mo/W-S bonds during sulfur deposition on the alloyed S(Mo W )(W Mo )S structure. Sulfur atoms facilitate simultaneous extraction of both Mo and W atoms, ultimately forming an alloyed bilayer Mo W S configuration. A cut-off distance of 2.6 Å was applied for both Mo-S and W-S bond identifications. The bond counts include both pre-existing and newly formed bonds in the system. The initial number of bonds reflects the presence of metal–sulfur coordination in the intermediate structures, and as the bilayer TMD forms, the number of Mo–S bonds nearly double—indicating the completion of coordination environments in the emerging bilayer structure.”\n\nComment 5: The revised description of the simulations still requires careful refinement. It is not, yet, fully clear what is done and in the current description, terminology of ensembles is mixed up.\n\nAuthor reply: We thank the reviewer for pointing out the ambiguity in our description of the simulation ensembles. After careful review, we confirm that, except for the AIMD validation of Mo-atom embedding shown in Supplementary Figure S7(a), which employs the NVE ensemble, all other MLP-MD and AIMD simulations were carried out in the canonical (Nosé–Hoover NVT) ensemble. To avoid any further confusion, we revised our manuscript to state clearly the ensemble and temperature for each simulation.\nIn the revised manuscript, we have added the following Table, on Supplementary Information page 7:\n“Table S3: Summary of simulation parameters used in each figure\nFigure Reference Temperature Ensemble\nFig. 2(d-e) 900 K NVT\nFig. 3(a-f) 900 K NVT\nFig. 4 1100 K NVT\nFig. S4 900 K NVT\nFig. S6 1100 K NVT\nFig. S7(a) 1100 K NVE\nFig. S7(b-c) 1100 K NVT\nFig. S9 1100 K NVT\nFig. S10 300 K NVT\nFig. S11 1100 K NVT\nFig. S13 1100 K NVT\nFig. S14 1100 K NVT\nFig. S15 900 K NVT\nFig. S16 900 K NVT\nFig. S17 1100 K NVT\nFig. S18 1100 K NVT\n”\n\nComment 6: The discussion of the manuscript results still requires better structure, a lot of back-and-forth is contained around Fig. 3 making it hard to follow the authors argumentation. After revision, new questions about the simulations occurred. For example, a key result is the fact, that “metal films do not float” (compared to the original experimental data in Nature 621, 499 15 (2023)), what is the significance of this finding?\n\nAuthor reply: Thank the reviewer for concern about our discussion. We have reorganized the description of Figure 3 in the “Results and Discussion” section into four steps—empirical observation, mechanistic discussion, kinetic validation and statistical analysis. Specifically, at 900 K, our MLP-MD simulations show that Mo atoms deposited on the WS surface quickly diffuse through the top S-layer, which contrast the sharply with the “coating” hypothesis of Nature 621, 499 (2023). This indicates that the growth mechanism presented in the previous study is questionable, and thus the route towards the synthesis of wafer-scale vdWHs requires more experimental validation.\nIn the revision, we expanded discussion on the “metal films do not float” finding to address its significance and possible reasons for divergence from experimental observations. On manuscript page 11:\n“A recent study reported23 that a two-step vapor deposition process was used to synthesize wafer-scale TMD vdWHs, where the key intermediate structure during growth is a Mo atomic monolayer on the WS surface. However, our MLP-MD simulations reveal that the Mo atomic monolayer on the WS surface is highly unstable, it can be easily transformed into the SMoWS structure by thermal annealing. Then the growth kinetics is actually different from that presented in the reference 23.”\n\nComment 7: The results also contain quite a bit of inconsistent simulations: The MD simulations are conducted around deposition temperatures of 500K-1300K (?), but some analysis (e.g. mixing entropy), and other simulations (e.g. depicted in Fig. R7) are conducted at 300K. This seems inappropriate if the synthesis mechanism is investigated. Some arguments are made about alloying where WS and MoS slabs appear not systematically exchanged and kinetic stabilization is argued on long MD timescales (“stability during 1.4 ns”) but without a significant number of samples. These simulations thus only weakly support the authors’ conclusions.\n\nAuthor reply: We thank the reviewer for this important and detailed comment. We understand the reviewer’s concerns regarding the use of different temperature settings in our simulations. This issue may stem from the presence of multiple simulation modules in the manuscript, each serving distinct purposes. We clarify the physical background and motivation for the different temperature choices as follows:\nThe 500–1300 K temperature range was used exclusively in the “on-the-fly” sampling phase of AIMD, aimed at constructing the training dataset for the MLP (as detailed in SI Section 1). These simulations were designed to sample a wide range of atomic configurations that may occur during growth processes and were not intended to reflect the actual physical conditions of deposition.\nSimulations conducted at 300 K serve two specific purposes. First, this temperature was used as the reference in our mixing free energy calculations, allowing us to assess the thermodynamic stability of alloyed versus non-alloyed structures under ambient conditions. Second, one MLP-MD simulation at 300 K (Fig. S10) was performed to examine whether Mo–W atomic exchange could occur at low temperature. The result showed that such exchange is kinetically suppressed at 300 K, which provides a plausible explanation for the formation of non-alloyed interfaces observed in experiments.\nRegarding the reviewer's comment that we lack systematic simulations of MoS and WS structure substitutions, we appreciate this insightful suggestion. We have supplemented the revision with simulations of the behaviour of Mo-based clusters on WS , and Fig. R3 presents the results under the same simulation conditions (1100 K) as those used for MoS .\nFurthermore, concerning the reviewer's concerns about the MD simulation time (1.4 ns) and the limited number of samples, we would like to further clarify the objectives and simulation strategy of this work. Our research is not aimed at accurately calculating diffusion coefficients or constructing complete statistical distributions but rather focuses on the differences in the dynamic behaviour of different precursor clusters on the TMD surface, thereby revealing their structure-dependent behaviour and potential reaction pathways.\nDespite the limited total simulation time, we still observed clear cluster-dependent dynamic behaviour. As shown in Fig. R3 (e, j, o, t), different types of Mo-based clusters exhibit noticeably different in-plane migration behaviours on the WS surface:\nMo single atom undergoes localized motion and are trapped in the initial region within a short time (Fig. R3e);\nMo–S and Mo–S clusters exhibit limited diffusion on the surface, with lower mobility (Fig. 2 R3j, Fig. R3o);\nIn contrast, Mo–S₃ clusters exhibit more significant surface mobility and irregular migration paths (Fig. R3t).\nTo further verify the influence of temperature on dynamic behaviour, we also conducted similar simulations at 900 K (see Fig. R4). The results show that the diffusion behaviour of Mo–S and Mo–S clusters is significantly reduced at the lower temperature, which contrasts with the migration trends under high-temperature conditions.\nFigure R3: MLP-MD simulations of various Mo-based clusters on WS (1100 K): (a-d) Single Mo atom deposition; (e) Mo atom trajectory on the xy plane; (f-j) Mo-S cluster; (k-o) Mo-S cluster; (p-t) Mo-S cluster. 3\nFigure R4: MLP-MD simulations of various Mo-based clusters on WS (900 K): (a-d) Single Mo atom deposition; (e) Mo atom trajectory on the xy plane; (f-j) Mo-S cluster; (k-o) Mo-S cluster; (p-t) Mo-S cluster. 3\nIn summary, although extending the simulation time or increasing the sampling of initial conditions could further improve quantitative accuracy, in the current setup, we have been able to identify the migration trends and relative activity of different cluster structures on the TMD surface. These results provide meaningful physical insights into the structural evolution and surface reconstruction of metal precursors in the early stages of deposition.\nIn the revised manuscript, we have added the following sentences, on manuscript page 14:\n“To deeply understand the dynamic behaviour of Mo atoms and their clusters with S atoms on various TMD substrates, we conducted further MLP-MD simulations to investigate the migration mechanism and stability of Mo clusters on WS surfaces. We found that their behaviour was similar to that observed on MoS substrate (Fig. S14). Fig. S15 and S16 illustrate the behaviour of single Mo atoms and Mo-S clusters on MoS and WS surfaces at a lower temperature (900 K). The results indicate that, even at reduced temperatures, single Mo atoms still tend to rapidly embed into the substrates. However, the surface mobility of Mo-S and Mo-S clusters is significantly reduced compared to their behaviour at 1100 K. In contrast, Mo-S clusters maintain relatively high surface stability and mobility, further validating the critical role of sulfur-containing clusters in suppressing Mo atom sinking (Fig. S15, Fig. S16). Fig. S17 further confirms the reliability and accuracy of these MLP-MD results through AIMD simulations at 1100 K for 10 ps. Although the simulation duration was shorter due to computational limitations of AIMD, the clusters displayed similar characteristics as observed in the MLP-MD simulations.”"}, {"round": 3, "reviewer_comments": {"reviewer_1": "In this revision, the authors have satisfactorily addressed all remaining concerns. The validation of the MLP and the strengthened statistics of the non-equilibrium MD raise confidence in the simulations and enhance reproducibility. The contextualization with respect to reference [23] improves the overall understanding of the predictions relative to experiment. Finally, I commend the improved writing in this revision. I am pleased to recommend this work for publication now.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Author reply: We thank the reviewer for the constructive feedback and for acknowledging the significant improvements in the manuscript’s transparency and traceability. We have carefully addressed the concerns regarding reproducibility, contextualization with experimental findings, and presentation. Below is our point-by-point response.\n\nComment 1: MLIP training and validation: The inclusion of the MLIP training and validation details is appreciated and now allows the reader to follow the methodology. The reported RMSEs appear acceptable considering the high-energy structures accessible at elevated temperatures. Nevertheless, to align with state-of-the-art reporting, I suggest the authors should include a measure of convergence — either RMSE versus “training round” or RMSE as a function of structure energy per atom — to assess model quality comprehensively. I also recommend revising the newly added paragraphs for clarity, as the current language is rough (especially the first one). Units for energy RMSE (meV vs. meV/atom) must be reported consistently.\n\nAuthor reply: We appreciate the reviewer’s valuable suggestions regarding the validation of our machine learning potential (MLP). To comprehensively assess the model quality, we adopted the reviewer’s recommendation to analyse the RMSE of energy as a function of the formation energy per atom (calculated relative to the elemental references), as shown in Fig. R1. We believe this metric is particularly informative as it demonstrates the model’s accuracy across a diverse configuration space, ranging from stable low-energy structures to high-energy non-equilibrium states encountered during growth. We have included this analysis in the Supplementary Information [Figure S1(c)]. The results reveal that the model achieves high accuracy for stable structures and maintains prediction errors within a reliable range (with energy RMSE generally remaining below 20 meV/atom) even for high-energy configurations. Furthermore, we have revised the paragraphs describing the MLP training in the Supplementary Information to improve flow and clarity. We have also carefully checked the entire manuscript and standardized all energy RMSE units to meV/atom and force RMSE units to meV/Å to ensure consistency.\n\nFigure R1: Distribution of energy RMSE for test set structures across varying formation energy intervals. The histograms display the prediction errors for energy in meV/atom. The formation energy per atom determines the binning interval.\n\nOn Page 5 of manuscript, we added the related discussions to assess model quality comprehensively.\n“To comprehensively assess the model's reliability across the energy landscape, we analysed the prediction errors as a function of structure formation energy [see Fig. S1(c)]. This analysis confirms that the MLP maintains high accuracy not only for stable low-energy configurations but also for high-energy states relevant to growth simulations.”\n\nOn Supplementary Information page 4:\n“FIG. S1(c) provides a more detailed error analysis by categorizing the test structures based on their formation energy per atom. The histogram displays the energy RMSE distribution across different formation energy intervals. This analysis reveals that the model maintains exceptional accuracy for structures in the lower energy regimes (typically corresponding to equilibrium or near-equilibrium configurations). Crucially, even for higher-energy states, which represent distorted configurations or transition states frequently encountered during high-temperature growth simulations, the energy prediction errors remain low and bounded. This comprehensive validation confirms the MLP's robustness and consistent accuracy across the entire relevant energy landscape.”\n\nOn Page 1 of Supplementary Information, we revised the related paragraphs as the reviewer suggested:\n“In this study, we employed an iterative training approach to develop the MLP. The process began with the generation of an initial dataset using ab initio molecular dynamics (AIMD) simulations. Subsequently, we assessed the model’s reliability by benchmarking it against DFT calculations and applying it to simulate the two-step growth of MoS /WS . If the MLP exhibited instabilities during molecular dynamics or failed to reproduce accurate phonon spectra, we refined the model by incorporating additional relevant DFT configurations into the training set. This cycle of data expansion and model retraining was repeated iteratively until the MLP achieved satisfactory accuracy and stability. The final optimized model was then deployed for efficient and accurate large-scale atomistic simulations.\nTo generate the initial training set, we employed on-the-fly molecular dynamics (MD) simulations, a method that effectively explores unknown configuration spaces while minimizing sampling costs. These simulations were conducted using the VASP on-the-fly MD framework with parameters optimized for efficiency and accuracy. Specifically, we utilized a timestep of 3.5 fs to capture atomic motions precisely, with a maximum of 20,000 steps per trajectory. The simulations were performed in the canonical (NVT) ensemble using Nosé-Hoover thermostats, covering a temperature range from 500 K to 1300 K to encompass typical growth conditions. The unit-cell shapes were kept fixed, and the number of local reference configurations in the on-the-fly MLP was constrained by setting ML_MB=2000. The initial configurations included monolayer MoS and WS , as well as their combinations featuring homogeneous/heterogeneous junctions and adsorbed metal clusters. Subsequently, the training set was expanded through iterative training rounds employing distinct sampling strategies:”\n\nComment 2: Reproducibility: I strongly disagree with the authors’ assertion that performing multiple simulations is unnecessary. Repetition under varying initial conditions is essential, particularly for non-equilibrium simulations, to validate qualitative insights and ensure reproducibility. This is not about obtaining quantitative observables but demonstrating that observed mechanisms are robust. Performing multiple simulations could result in a brief statement such as, “All qualitative observations were confirmed in x independent runs” and would meaningfully strengthen the manuscript. Testing one’s own hypotheses in this way is fundamental to good scientific practice.\n\nAuthor reply: We acknowledge the reviewer’s concern regarding reproducibility and have performed additional simulations to demonstrate the robustness of our qualitative insights. To address this, we performed three independent MD simulations with different random velocity seeds, focusing on the primary non-equilibrium dynamic processes that underpin our main conclusions. Regarding the deposition process, as shown in Fig. R2, Mo atoms invariably sank into the chalcogen layer within similar timescales across all three independent runs, confirming that the formation of the embedded intermediate is a robust energetic preference. Similarly, for the sulfurization process shown in Fig. R3, all independent runs confirmed that the deposition of excess S atoms successfully pulled metal atoms to the surface, forming the bilayer structure. These results confirm that our observed mechanisms are highly robust and reproducible. We have added a statement in the \"Results and Discussion\" and \"Methods\" sections to attest to this verification.\nFurthermore, a very recent study published in Science [Gao et al., Science 390, 813–818 (2025)] reported an \"atomic layer bonding (ALB)\" contact (Fig. R4). This study realized a direct metallic bonding interface between metal electrodes and the transition-metal layer of TMDs by removing the top sulfur layer. The physical nature of this ALB structure (direct metal-metal bonding) is highly consistent with our predicted SMMS intermediate. Crucially, the experimental ALB contact was shown to be metallic and exhibited superb thermomechanical stability and ultralow contact resistance. This experimental evidence strongly supports our theoretical prediction that the SMMS structure serves as a promising electrode candidate. We have incorporated this discussion into the revised manuscript to further substantiate the physical validity and application potential of our identified intermediate.\n\nFigure R2: Snapshots of Mo atom deposition on MoS from 900 K MLP-MD simulations using three different random seeds: (a-d) Run 1, (e-h) Run 2, and (i-l) Run 3.\n\nFigure R3: Snapshots of sulfurization on the alloyed S(Mo W )(W Mo )S structure from 1100 K MLP-MD simulations using three different random seeds: (a-d) Run 1, (e-h) Run 2, and (i-l) Run 3. The S atoms consistently extract Mo/W atoms to form a bilayer alloyed Mo W S structure.\n\nFigure R4: Comparison of our predicted structure with the experimental ALB contact. (a) The SMoMoS structure identified in our simulations. (b) The ALB structure reported in Gao et al., Science (2025). (c) Cross-sectional STEM image of the experimental ALB structure, showing direct metal-metal bonding similar to SMoMoS.\n\nIn the revision, we added the following sentence on manuscript page 20:\n“To ensure the reproducibility of the non-equilibrium processes described above, all key qualitative observations (including the sinking of Mo atoms and the sulfurization-induced extraction) were confirmed in three independent MD simulations initialized with different random velocity seeds.”\n\nOn manuscript page 17:\n“Notably, a very recent experimental study43 realized an “atomic layer bonding (ALB)” contact by establishing a metallic coherent bonding interface between the transition-metal layer of TMDs and metal electrodes. This ALB structure shares the key feature of direct metal-metal bonding observed in our predicted SMMS intermediate. Consistent with our findings that SMMS exhibits metallic character and improved contact properties, the experimental ALB contact demonstrated ultralow contact resistance and superb thermomechanical stability. This strongly corroborates our theoretical prediction that establishing coherent metal-metal interactions (as in SMMS) is an effective strategy to overcome the limitations of van der Waals contacts and achieve high-performance electronic devices.”\n\nWe have added the following new reference, Ref. [43] in the revised manuscript:\n“43. Gao, L. et al. Atomic layer bonding contacts in two-dimensional semiconductors. Science 390, 813–818 (2025).”\n\nComment 3: Temperatures in MD simulations: The selection and reporting of MD temperatures remain inconsistent. Temperatures listed in Table S1 should also appear clearly in the main text and relevant figure captions. Moreover, the rationale behind each temperature choice must be stated — currently, these choices appear arbitrary. This information is essential for readers to assess the physical relevance of the simulations.\n\nAuthor reply: We apologize for the confusion regarding the temperature settings. We have now explicitly stated the simulation temperatures in the main text and all relevant figure captions. Furthermore, we have clarified the rationale for these choices in the Methods section, ensuring they systematically cover the two-step vapor deposition method reported in Ref. [23]. Specifically, 1100 K was selected to align with the sulfurization temperatures (Step 2) of the experiment, where MoS are synthesized at 800 °C (1073 K). 900 K was chosen to correspond to the pre-heating or baking stage (experimentally >600 °C) and the heating ramp-up process prior to full sulfurization, allowing us to assess the stability of the metal layer as the system undergoes thermal treatment before reaching the peak reaction temperature.\n\nIn the revision, we revised the following sentences, on manuscript page 7:\n“Subsequently, individual Mo atoms were simulated to be continuously sputtered onto the existing MoS layer with a kinetic energy of 0.12 eV. To simulate the stability during the metal atom deposition and pre-heating stage, simulations were performed at 900 K.”\n\nOn manuscript page 19:\n“For Mo atom deposition, we used a 4√3×7 MoS or WS supercell at 900 K. This temperature was chosen to represent the thermal state during the pre-heating/baking stage and the ramp-up process prior to full sulfurization. The SMoMoS and SMMS sulfurized process simulations were performed using a 4√3×7 supercell at 1100 K, aligning with the experimental sulfurization temperatures (approx. 1073 K).”\nand added temperatures to all relevant figure captions (e.g., Fig. 2, Fig. 3, Fig. 4, etc.).\n\nComment 4: Reporting of reaction energies: Certain data remain ambiguously described or mislabeled. For instance, in Figure 2, the “sinking” of Mo is attributed to DFT calculations, although the associated energy change appears to originate from MD trajectories and was only later corroborated by DFT single-point calculations according to the SI(?). The authors must clarify which quantities stem from MD and which from DFT to prevent misinterpretation, alternatively, a description of the DFT calculations is required.\n\nAuthor reply: We thank the reviewer for pointing out this ambiguity regarding the reporting of reaction energies. In the revised manuscript, we clarified that while the dynamic sinking process is illustrated by MLP-MD, the cited energy release of 1.45 eV was derived from DFT relaxations of the initial and final structures to ensure quantitative accuracy. Furthermore, all static reaction energy values reported throughout the manuscript are based on DFT calculations to ensure reliability (e.g., the energy differences for MoSMoS to SMoMoS, and MoSWS to SMoWS or SMMS).\n\nIn the revision, we revised the following sentences, on manuscript page 7:\n“First, the adsorbed single Mo atom was unstable on the MoS layer [Fig. 2(a)] and, at the typical growth temperature of 1100 K, quickly sank beneath the S atom layer within tens of picoseconds [Fig. 2(b)]. DFT calculations confirm that this embedding process releases 1.45 eV of energy.”\nand\n“DFT calculations reveal that the energy decreases significantly by 1.73 eV per Mo atom when transitioning from MoSMoS to SMoMoS, indicating a strong thermodynamic driving force towards the formation of the SMoMoS structure.”\n\nOn manuscript page 10:\n“Conversely, the Mo atoms completely sink and embed under the top S layer of the WS 2 monolayer, forming the SMoWS intermediate structure [Fig. 3(h)], releasing 2.08 eV per Mo atom (based on DFT calculations).”\n\nComment 5: Interpretation of results: The authors’ additional discussion connecting this work to Ref. [23] is appreciated, but some confusion persists. It would be clearer to reference this work in the new paragraph as “In the above-mentioned study [23]” rather than “A recent study reported…” to guide readers properly, after all [23] is mentioned multiple times as a motivation for this study. More importantly, the authors’ findings still seem to contradict the experimental conclusions of Ref. [23], which reported successful wafer-scale synthesis of TMDs after direct dosing of Mo on WS . The authors could strengthen their argument by discussing potential reconciling factors. I would think a possible interpretation are sulfur-rich conditions that might have been present during experiments in [23] and thus stabilized against the intermediate SMeMeS configuration, in accordance with the authors discussion in another context and findings around Figures S13–S16.\n\nAuthor reply: We appreciate the reviewer’s insightful suggestion on reconciling our simulation results with the experimental findings in Ref. [23]. We have adopted the suggested phrasing “In the above-mentioned study [23]” in the revised manuscript to improve flow and clarity.\nRegarding the experimental observations in Ref. [23], we note that the Mo film is deposited after the growth of the first WS layer. Given the sulfur-rich environment of the preceding step, it is highly likely that residual sulfur remained in the growth chamber or on the substrate surface. This residual sulfur ensures sufficient sulfur availability, allowing Mo atoms to be partially sulfurized or pinned (forming Mo-S clusters) during deposition or the subsequent heating ramp-up, before they have a chance to sink. Our simulations (Fig. 4 and Fig. S13) explicitly confirm that once metal atoms bond with sulfur, their embedding behaviour is suppressed. Therefore, our work does not contradict the experimental results; rather, it provides an atomic-level mechanism: residual sulfur from the previous step is critical for suppressing the SMMS intermediate and preventing alloying. We have expanded the discussion in the revised manuscript to clarify this point.\n\nIn the revision, we revised the following sentence on manuscript page 11:\n“In the above-mentioned study23, a two-step vapor deposition process was reported to synthesize wafer-scale TMD vdWHs, where the key intermediate structure during growth is a Mo atomic monolayer on the WS surface. However, our MLP-MD simulations reveal that the Mo atomic monolayer on the WS surface is highly unstable, it can be easily transformed into the SMoWS structure by thermal annealing. This suggests that the experimentally observed clean interfaces must be governed by kinetics or environmental factors that prevent the formation of this bare metal intermediate.”\n\nOn manuscript page 14,\n“Experimentally, the sulfur source and molybdenum source ratios used during the MOCVD growth of MoS and WS are significantly higher than 2:1, such as 70:138, 660:139, 6400:140, and 11111:141. Similarly, regarding the successful synthesis in the above-mentioned study23, we propose that sulfur-rich conditions likely played a decisive role. Since the Mo film is deposited after the synthesis of the first WS layer, residual sulfur remaining in the growth chamber or on the surface from the preceding step is inevitable. This residual sulfur can react with deposited Mo atoms to form Mo-S clusters during the deposition or the initial heating phase, effectively suppressing the embedding of bare metal atoms and the formation of the SMMS intermediate, as predicted by our simulations. This mechanism explains how alloying is prevented in the experimental two-step process.”"}]

Decision letter