Microscopic Mechanism of Van der Waals Heteroepitaxy in the Formation of MoS2/hBN Vertical Heterostructures.

Recent studies have revealed that van der Waals (vdW) heteroepitaxial growth of 2D materials on crystalline substrates, such as hexagonal boron nitride (hBN), leads to the formation of self-aligned grains, which results in defect-free stitching between the grains. However, how the weak vdW interaction causes a strong limitation on the crystal orientation of grains is still not understood yet. In this work, we have focused on investigating the microscopic mechanism of the self-alignment of MoS2 grains in vdW epitaxial growth on hBN. Using the density functional theory and the Lennard-Jones potential, we found that the interlayer energy between MoS2 and hBN strongly depends on the size and crystal orientation of MoS2. We also found that, when the size of MoS2 is several tens of nanometers, the rotational energy barrier can exceed ∼1 eV, which should suppress rotation to align the crystal orientation of MoS2 even at the growth temperature.

Introduction

Two-dimensional (2D) materials have been attracting much attention due to their fascinating properties and possible applications in nanoelectronics and photonics.[1,2] 2D group-VI transition metal dichalcogenides (TMDs: MoS2, WS2, MoSe2, WSe2, etc.), in particular, are outstanding because they can have a sizable bandgap (1.5–2.0 eV) that is absent in graphene.[3,4] In addition, TMDs with the 2H phase, at the monolayer limit, possess a direct gap,[4,5] where the spin and valley degrees of freedom are coupled due to inversion symmetry breaking.[6] These properties make TMDs excellent platforms for exploring the fundamental physics at the 2D limit and for next-generation optoelectronic devices.[7]

One of the unique growth modes of 2D materials is van der Waals (vdW) heteroepitaxial growth, where the vdW force significantly limits the crystallographic orientation of 2D materials grown on crystalline substrates.[8,9] In the past few years, chemical vapor deposition (CVD) growth of monolayer TMDs on various crystalline substrates, including Au (111), Al2O3, graphene, graphite, and hexagonal boron nitride (hBN), has been reported.[10−20] Although there is a significant lattice mismatch between TMDs and the substrates (e.g., ∼26% in the case of MoS2 and hBN), the substrate–TMD interaction, vdW force, precisely aligns the TMD domains, limiting the stacking angles to 0 and 60°, that is, vdW heteroepitaxial growth.[11−13] The vdW heteroepitaxial growth is a universal phenomenon, which has widely been observed not only in the growth of TMDs but also in the growth of other layered materials, such as hBN and graphene.[21,22]

“How the vdW force can strongly restrict the crystal orientation in vdW heteroepitaxial growth?” is one of the mysteries in the growth of 2D materials. In the heteroepitaxial growth of three-dimensional (3D) systems such as III–V semiconductor compounds, chemical bonds between a grown material and a corresponding substrate determine the crystallographic orientation between them.[8] In this case, it is natural that strong chemical bonds perfectly define the crystal orientation to form commensurate heterostructures.[9] On the other hand, the weak nonbonding vdW force determines the stacking angle in the formation of incommensurate stackings in the vdW heteroepitaxial growth.[8,9] Because there is a significant lattice mismatch in most vdW heteroepitaxial growth, atoms in the growing layer cannot always locate at stable positions; the configuration of atoms between a substrate and a growing layer gradually and continuously changes due to the difference in lattice constants.[18] In this case, the stabilization energy would not simply increase depending on the size of a growing layer, and the optimum crystal orientation is not trivial. Moreover, growing layers are expected to move and rotate easily because there are no direct bonds between layers and substrates.[23]

In this work, we focus on the microscopic mechanism of the substrate-induced alignment in the CVD growth of monolayer MoS2 on hBN substrates. Detailed analyses on the relative orientation of MoS2 crystals and hBN substrates by electron diffraction have confirmed that relative orientations between CVD-grown MoS2 and the underlying hBN substrate are limited to 0 and 60°. To elucidate the mechanism of self-alignment, we have calculated cluster-size-dependent interlayer-energy landscapes. In addition to the density functional theory (DFT), we have also employed Lennard–Jones (LJ) pair potential to calculate the interlayer-energy landscape of the MoS2/hBN system because the LJ potential has allowed us to treat a system containing over 104 atoms, which is far beyond the scope of the DFT calculation. In both cases, we found that the interaction energy shows significant stacking angle and cluster size dependence. Although the stacking angle of 0/60° is not always the most stable configuration, we found, up to the size of Mo972S1944 (an edge length of 5.7 nm), that the most stable configurations with the stacking angle of 0/60° appear periodically, whose period is close to the moiré period between MoS2 and hBN. We also found that the rotational energy barrier (the difference between the interaction energy at 0° and the angle corresponding to the nearest local maximum to 0°) is comparable to thermal energy at growth temperatures when the cluster size reaches several tens of nanometers. These indicate that small clusters, at the early stage of growth, rotate to change the stacking angle, and then, the rotational degree of freedom is gradually lost as the cluster size increases. This finding provides a basis to get deep insights into the growth mechanism of 2D films onto these crystalline substrates.

Results and Discussion

To grow MoS2 onto hBN, we used the multifurnace CVD method with elemental sulfur and molybdenum trioxide as precursors. The hBN substrate is obtained using the mechanical exfoliation method with a bulk hBN crystal grown using the high-pressure, high-temperature method.[34] Details of the growth procedure are shown in the Methods section and a previous report.[24]Figure a shows a typical optical image of MoS2 grown on hBN (MoS2/hBN), where the hexagonal-shaped contrasts correspond to single-layer MoS2.[24,35] The grown MoS2 grains with a typical size of ∼5 μm are placed randomly on the hBN substrate; a schematic image of MoS2/hBN is shown in Figure b. To confirm the layer number of MoS2, we measured the Raman spectra at room temperature with an excitation energy of 2.33 eV. Figure c shows a Raman spectrum, where two characteristics peaks (383.7 and 404.7 cm–1) originating from in-plane (E′ mode) and out-of-plane vibrational modes (A′1 mode) exist.[36] The frequency difference between these two peaks is 21.1 cm–1, which is slightly larger than the reported value of monolayer MoS2.[37]Figure d shows a typical photoluminescence (PL) spectrum of MoS2/hBN at room temperature. As seen in Figure d, a single-peak PL emission centered at 1.884 eV is observed, consistent with previous reports.[38,39] The full width at half maximum of the obtained PL peak is about 43 meV, which is significantly smaller than that of monolayer MoS2 grown onto a SiO2/Si substrate (55 meV),[39] indicating that inhomogeneous broadening arising from substrates is suppressed in MoS2/hBN.[24] Atomic force microscopy (AFM) observations shown in Figure S1 are also consistent with monolayer MoS2.[12]

Figure 1

(a) Optical image of CVD-grown MoS2 on hBN; (b) schematic of MoS2/hBN; (c) typical Raman spectrum, (d) PL spectrum of the MoS2 crystal shown in the upper left side of Figure a; and (e) a typical SAED pattern of MoS2/hBN. Green and blue arrows indicate diffraction spots from MoS2 and hBN in MoS2/hBN, respectively.

As shown in Figure a, we can see that all MoS2 domains have only two orientations, where the 60° rotation of one orientation matches the other orientation. The observed limitation in the orientation of MoS2 results from a strict relationship in the crystallographic orientation between MoS2 and hBN. In fact, a typical selective-area electron diffraction (SAED) pattern shown in Figure clearly shows two sets of six-fold-symmetric spots with almost the same orientation; the larger (smaller) six-fold-symmetric pattern originates from MoS2 (hBN).[10,11] This means that the relative crystallographic angle between MoS2 and hBN is limited to two orientations, 0 and 60°, and statistics on the frequency of the two stacking angles are summarized in Table (a corresponding optical image is shown in Figure S2). As clearly seen, the ratio between the two stacking angles is nearly 1:1, which is consistent with previous studies.[11] The slight deviation from 1:1 probably originates from statistical error or defect-controlled nucleation of MoS2 on hBN.[40]

Table 1

Distribution of MoS2 Grain Orientation

structure	aligned—0°	aligned—60°	others
number of grains	60	46	0

To investigate the mechanism of the vdW heteroepitaxial growth of MoS2/hBN, we calculated the cluster-size-dependent interaction-energy landscapes by DFT. We calculated the interaction-energy landscapes for four MoS2/hBN with different cluster sizes; the corresponding models are shown in Figure a. The detailed calculation method is provided in the Methods section. Figure b shows the total energy of MoS2 on hBN as a function of the stacking angle. As clearly seen, the energy is sensitive to the cluster size and the stacking angle. For the smallest MoS cluster, the cluster prefers a stacking angle of approximately 15°, reflecting the local atomic arrangement between MoS and hBN. The increase of the cluster size causes additional energy minima, owing to increased preferential interlayer atomic arrangements arising from the lattice mismatch. Indeed, for the largest cluster, we found two global minima at the angles of 10 and 60° in addition to the two local minima, still reflecting the local atomic arrangement between them. Thus, with the further increase in the flake size, the MoS cluster exhibits some particular orientations by the averaged vdW interaction. Note that the stacking structures with a stacking angle of 0 or 60° do not correspond to the most stable structures since we have not optimized the relative position of MoS and hBN, such as interlayer distances and locations of the center of gravity.

Figure 2

(a) Schematics of the structure model used in DFT calculations. The coloring of the elements is the same as that used in Figure b. (b) Stacking-angle-dependent total energy calculated with different cluster sizes. The total energy at the most stable stacking angle is set to zero.

To have a more in-depth insight into the mechanism of the vdW-mediated epitaxial growth of MoS2 onto hBN, we extended the calculation of interaction-energy landscapes to larger clusters with the LJ pair potential. A simple pair potential, such as the LJ potential, allows us to calculate the interaction energy between a cluster of MoS2 and hBN, which are composed of a large number of atoms. In this work, we calculated large MoS2/hBN-stacked structures up to Mo972S1944/B5400N5400, which is beyond the scope of DFT calculations. In the LJ pair potential, the interaction energy between a TMD cluster and an hBN substrate, Einter, is described as a summation of the LJ potential.NTMD and NhBN are the number of the atoms in TMDs and hBN, respectively. ε and σ are the LJ parameters corresponding to interactions between the ith and jth atom in the system, and r is the distance between the ith and jth atom in the system. We calculated the stacking angle-dependent energy landscapes with various sizes of MoS2. Note that MoS2 and hBN are treated as rigid bodies, and intralayer energy, such as bending chemical bonds, is ignored. The LJ parameters in eq were estimated with the Lorentz–Berthelot combining rules, where σ and ε between different kinds of atoms, σAB and εAB, were calculated with those between the same kind of atoms, σAA, σBB, εAA, and εBB, as follows.

The LJ parameters of Mo–Mo, S–S, B–B, and N–N (σMoMo, σSS, σBB, σNN, εMoMo, εSS, εBB, and εNN), which are needed to calculate σMoB, σMoN, σSB, σSN, εMoB, εMoN, εSB, and εSN, are taken from references;[41,42]Table shows the obtained LJ parameters.

Table 2

LJ Parameters Used to Calculate the Potential Curves and Images Shown in Figures –5, S3–S8

interaction	ε (meV)	σ (Å)
Mo–B	58.73	3.002
Mo–N	72.56	2.958
S–B	15.79	3.411
S–N	19.51	3.367

To confirm the validity of the LJ parameters obtained, we have checked the consistency between the interaction-energy landscapes calculated with the LJ pair potential and the ab initio DFT calculation. Figure S3 shows the stacking angle dependencies of interaction energy derived with the LJ pair potential; a large hexagonal bilayer hBN (B2700N2700 × 2, an edge length of 7.5 nm) and MoS clusters used in DFT calculations are employed in this calculation. As clearly seen, calculations based on the LJ potential qualitatively reproduce the DFT results (Figure b). Having reproduced the DFT results successfully, we conducted LJ potential-based calculations to obtain the interaction energy of MoS2/hBN up to Mo972S1944/B5400N5400.

Figure S4 shows the scheme to calculate the interaction-energy landscapes of MoS2/hBN. First, we put a hexagonal Mo3S6 (n = 1, 2, ..., 18; the relationships between n and the edge length of each cluster are shown in Table S1) cluster on a hexagonal bilayer hBN (B2700N2700 × 2) with a certain stacking angle (Figure S4 (i)), and then, we search for the ground minimum configuration at the stacking angle using the simulated annealing (SA) method (Figure S4 (ii, iii)). The parameters in the SA calculations are lateral positions and interlayer distances. The SA-based ground-minimum search has been repeated with various stacking angles to calculate an interaction-energy landscape for a certain cluster size of MoS2 (Figure S4 (iv)). This process has been repeated with various sizes of clusters and the obtained interaction-energy landscapes for 18 different clusters of MoS2 are shown in Figure a; the most stable configurations found for a stacking angle of 0° are shown in Figure S5 and Table S1.

Figure 3

(a) Cluster-size and stacking-angle evolutions of interaction energy; (b) cluster-size dependencies of absolute values of interaction energy (i.e., absolute values of the difference between the maximum energy and the minimum energy) calculated with a stacking angle of 0 and 60° and the most stable stacking angle of each cluster; (c) map showing element-decomposed interaction energies of a MoS2 cluster (an edge length of 1.6 nm) with a stacking angle of 0°. Yellow circles correspond to S2 pairs; and (d) corresponding structure used to calculate the element-decomposed interaction energy in (c). Purple lines indicate the moiré superlattice period. Element coloring is the same as that of Figure b.

As shown in Figure a, the interaction-energy landscapes show that several local minimums locate at different stacking angles, and the different cluster sizes of MoS2 give different positions of local minimums. Notably, structures with stacking angles of 0 or 60° do not always correspond to the ground minimum in the corresponding interaction-energy landscapes. The cluster-size dependency of the interaction energy at stacking angles of 0 and 60° oscillates periodically (Figure b); clusters with an edge length of 0.3, 1.6, 2.8, 3.8, and 5.1 nm have considerable interaction energy compared with those of clusters with a similar size. It should be noted that the periodicity of oscillations in interaction energy seen in Figure b is roughly close to that of the moiré superlattice in MoS2/hBN with stacking angles of 0 and 60°; the periodicity of the moiré superlattice is 1.196 nm; the moiré lattice constant λ and the relative rotation angle θ of the moiré pattern with respect to the hBN lattice can be calculated with the following equations[43]where δ, a, and ϕ represent the lattice mismatch between MoS2 and hBN, the hBN lattice constant, and the relative rotation angle, respectively. This can be understood by visualizing the contribution of each atom in MoS2 clusters to Einter. Figure c shows the contribution of each S2 pairs to Einter in the cluster with an edge length of 1.6 nm; the contribution is evaluated with the sum of the LJ potential energy between a S2 pair in the cluster with the total atoms in hBN. As shown in Figure c,d, S2 pairs that contribute to enlarge Einter appear with the moiré superlattice periodicity because the stable configuration between S2 pair and hBN, which significantly contributes to Einter, appears with the moiré superlattice period.

In the calculation of the interlayer-energy landscapes, we searched for the most stable configuration at each stacking angle. If there is a substantial energy barrier between the most stable configuration and configurations corresponding to a local minimum, MoS2 clusters cannot find the most stable configuration. To see this, the interaction-energy landscape on the position of a cluster (an edge length of 5.1 nm) was calculated (Figure S7). As clearly seen in the figure, the potential barrier does not exist for the translational motion of the cluster, and this means that the MoS2 cluster can move to find the most stable configuration. This conclusion can be extended to larger MoS2 clusters, which contain a larger number of moiré units because interaction energy can be roughly proportional to the number of moiré units, and cluster size is not expected to alter the shape of the potential profile (related discussion is given in the next paragraph). The significantly small energy barrier for translational motion in large systems, called superlubricity, has been observed in incommensurate systems, such as twisted graphite flakes,[44,45] and large MoS2 clusters are also expected to move on an hBN substrate. These results strongly suggest that, at the early stage of the growth, small MoS2 clusters would rotate and move to find the most stable configuration. Further discussions about the energy landscape, which include results on typical triangular-shaped clusters, are provided in the Supporting Information (Figure S8 and discussion therein).

As discussed above, MoS2 clusters on hBN move and rotate to find the most stable configuration. The next question is that “when does a MoS2 cluster stop rotating and align?” To address this question, we focused on cluster-size-dependent evolution of the rotational energy barrier at fixed stacking angles; rotational energies were calculated as the difference between the interaction energy at 0o and the angle corresponding to the nearest local maximum to 0°. As shown in Figure b, the stable and unstable atomic configurations appear periodically with the moiré period, and hence, rotational energy barriers can roughly be divided into the moiré contribution (Emoiré × number of moiré unit cell) and the edge contribution. While the edge contribution cannot be neglected compared to the moiré contribution when the cluster size is small, the moiré contribution becomes more dominant when the cluster size becomes large; the moiré/edge contribution should depend quadratically/linearly on the cluster size. Therefore, the rotational energy barriers of large MoS2 clusters and hBN can be estimated by the moiré contribution alone. Figure a shows the stacking angle dependence on the total energy of a MoS2 cluster with an edge length of 2.8 nm. As shown in the figure, there is a global minimum and a local maximum at 0 and 4.5°, respectively. Corresponding mappings of the interaction energy are shown in Figure b,c, and we extracted 10 moiré unit cells for each cluster; moiré unit cells are shown as dotted lines in Figure b,c. In this case, we can express the interaction energy as 10Emoiré + (the edge contribution), and Emoiré has successfully been extracted as ∼0.1 meV/(MoS2 unit cell). Based on the obtained Emoiré of ∼0.1 meV, we have estimated the rotational energy barrier of a large cluster containing ∼1 × 104 MoS2 unit cells (a cluster size of ∼40 nm or an edge length of ∼20 nm) as ∼1 eV. Compared with the thermal energy of a typical CVD growth temperature (∼100 meV), the rotational energy barrier is large enough, thereby suppressing the rotation of clusters.

Figure 4

(a) Stacking angle-dependent interaction energy of a Mo243S486 cluster (the interaction energy at the most stable stacking angle is set to zero) and (b, c) maps showing element-decomposed interaction energies of S2 pairs in a MoS2 cluster with an edge length of 2.8 nm with stacking angles of 0° and 4.5°. Yellow circles correspond to positions of S2 pairs and magenta dotted lines correspond to the moiré superlattice period of each structure.

A possible discrepancy between the calculations based on LJ potential and the experimental results is local/global energy minimums around 10–20 and 40–50° in the potential landscapes (as shown in Figure b). However, these stacking angles have not been experimentally observed as shown in Figure a and Table . This discrepancy is probably explained by the variability of the stable stacking angle and the translational position. As shown in Figure a, the stacking angles at the local/global minimums around 10–20 and 40–50o significantly change depending on the cluster size. Furthermore, the cluster size evolution also alters the stable translational position of MoS2 on hBN. To see this, we have calculated the stacking angle-dependent interaction energy evolution of the MoS2 cluster with fixed center positions; the center of the MoS2 clusters locates on B or N atoms (Figure b). As clearly demonstrated, both the stable stacking angle and translational position vary depending on the cluster size, while the stacking angle of 0/60° always corresponds to the local or global minimum (Figure a). Even when a cluster is trapped at the stacking angles of 10–20 and 40–50°, the energy minimum disappears when the cluster size slightly increases. For example, a cluster with an edge length of 2.5 nm can be trapped at a stacking angle of 16°, but the stacking angle does not correspond to one of the energy minimums anymore when the cluster size increased to have an edge length of 2.8 nm (Figure b). The cluster, therefore, is forced to move and rotate to find a stable configuration. Because the cluster keeps moving and rotating during a growth process unless the stacking angle locates at 10–20° and 40–50°, the cluster probably escapes the energy minima at 10–20° and 40–50°, leading to falling into the experimentally observed stacking angle of 0/60°. Note that the lateral positions of stable structures with a stacking angle of 0/60° also change as the cluster size increases (Table S1). Clusters are, however, expected to find the stable lateral position easily because of the absence of a translational energy barrier in this case (Figure S7). Furthermore, the stacking angle of 0/60° always corresponds to, at least, one of the local minimums, and there is no significant local energy minimum nearby the stacking angle of 0/60° (Figure a). In conjunction with the sizeable rotational energy barrier (∼1 eV in a cluster size of ∼40 nm), these characteristics of cluster-size- and stacking-angle-dependent energy landscape probably make the experimentally observed stacking angle of 0/60° stable.

Figure 5

(a) Cluster-size-dependent stable stacking angle evolutions around 40–50° (upper panel) and 10–20° (lower panel); (b) stacking angle and structure-dependent interaction-energy evolution of clusters with an edge length of 2.5, 2.8, and 3.2 nm. The left and right panels show the results calculated with stacking angles of 12–21° and 39–48°, respectively. Curves labeled as “SA result” are the same as the curves shown in Figure a. Black dotted lines correspond to an energy minimum position of a cluster with an edge length of 2.5 nm.

Conclusions

In summary, we investigate the mechanism of the vdW heteroepitaxy of MoS2/hBN through the calculation of interaction-energy landscapes. We found that the interaction-energy landscapes strongly depend on the size and stacking angle of clusters. We also found that the stacking angle of 0/60° always corresponds to a local or the ground minimum configuration. The energy barrier for the 0/60° configuration evolves as the cluster size grows, reaching ∼1 eV when the cluster size is several tens of nanometers. These findings suggest that (1) at the very early stage of the growth, a small MoS2 cluster can rotate to find stable configurations and (2) the rotational degree of freedom is gradually suppressed as the cluster grows and finally stops when the cluster size reaches several tens of nanometers. A similar phenomenon can be observed in other vdW heteroepitaxy when there is a lattice mismatch between a growing layer and a substrate. In lattice-mismatched systems, the interaction-energy landscape should show strong cluster-size dependence with a periodicity of the corresponding moiré lattice constant. Also, the interaction-energy landscape is expected to show significant stacking-angle dependence. These cluster-size and stacking-angle dependent landscapes should lead to the alignment of growing layers as already seen in the growth of TMDs onto graphite or the sapphire (0001) surface. Our results can be a basis of controlling and understanding the vdW heteroepitaxy of 2D films onto crystalline substrates.

Methods

CVD Growth of MoS2 on hBN

We have grown monolayer MoS2 onto exfoliated hBN flakes using the CVD growth method.[24] hBN flakes were prepared using the mechanical exfoliation method on a quartz substrate. As precursors, we used molybdenum trioxide (MoO3, Sigma-Aldrich, 99.5%) and elemental sulfur (Sigma-Aldrich, 99.98%). Sulfur powder and a quartz substrate with hBN flakes were loaded into a quartz tube with an inner diameter of 26 mm. MoO3 was placed in a quartz tube with an inner diameter of 10 mm, and the quartz tube was placed inside the larger diameter quartz tube to avoid unwanted reactions between S and MoO3, and then, under an Ar flow of 200 sccm, we heated the quartz tubes with a three-zone electric furnace at 200, 750, and 1100 °C for S, MoO3, and the substrate, respectively. The typical growth time is 20 min.

MoS2 Characterization

Optical images were taken with a typical optical microscope (Leica DM 2500 M and Nikon Eclipse ME600). Raman and PL measurements were performed using a confocal Raman microscope (Renishaw inVia) with an excitation energy of 2.33 eV. AFM measurement was performed using a standard AFM (Park Systems NX10) operating in a tapping mode. We used a transmission electron microscope (JEOL 2100) operating at 200 keV to obtain a SAED pattern. MoS2/hBN was transferred onto a copper grid using the standard polymer-based transfer method.

DFT Calculation

Theoretical calculations were performed using DFT[25,26] as implemented in the program package Simulation Tool for Atom TEchnology (STATE).[27] We used the generalized gradient approximation with the Perdew–Burke–Ernzerhof functional[28,29] to describe the exchange–correlation potential energy among interacting electrons. The weak dispersive interaction between TMD flakes and hBN was treated using the vdW-DF2 with the C09 exchange–correlation functional.[30,31]

Ultrasoft pseudopotentials generated using the Vanderbilt scheme were adopted as the interaction between electrons and nuclei.[32] Valence wavefunctions and the deficit charge density were expanded in terms of plane-wave basis sets with cutoff energies of 25 and 225 Ry, respectively. The atomic structures of TMDs were fully optimized until the force acting on each atom was less than 1.33 × 10–3 HR/au. We consider the triangular MoS flakes, Mo6S14, Mo10S24, Mo15S36, and Mo21S50, with S edges as the structural model of MoS flakes, which are adsorbed on monolayer hBN with a lateral supercell with the sizes of 5 × 5, 7 ×7, 8 × 8, and 9 × 9, respectively, possessing the lattice parameters of a0 = 1.256, 1.758, 2.009, and 2.226 nm (Figure a). To exclude the unphysical dipole interaction due to the electrostatic potential difference between MoS and hBN, we imposed an open boundary condition normal to the hBN sheet using the effective screening medium method,[33] with sufficiently large normal lattice parameters in which TMD/hBN is separated from adjacents by a 1.3 nm vacuum spacing. Integration over the Brillouin zone was carried out using an equidistance mesh of 2 × 2 × 1 k points. The calculation was performed as follows: (i) calculating stable MoS cluster structures; (ii) putting a MoS cluster on hBN at a certain position with a stacking angle of 0°; (iii) optimizing the interlayer distance between MoS and hBN; and (iv) calculating the stacking-angle dependence of the total energy by fixing the center of gravity of the MoS cluster and the interlayer distance.