Precise structure and energy of group 6 transition metal dichalcogenide homo- and heterobilayers in high-symmetry configurations

Two-dimensional group 6 transition metal dichalcogenide (2D TMDC) bilayers show various high-symmetry stacking configurations, which have also been observed in extended domains formed in their twisted homo- and heterobilayers. The interlayer energy varies for these stacking configurations, and the energy differences determine the relative size of the stacking domains. Therefore, the precise prediction of the composition- and stacking-dependent interlayer energy is crucial to model the domain structure of 2D TMDCs in their twisted bilayer homo- and heterostructures. For the validation of approximate methods that are necessary to tackle these systems encompassing thousands of atoms precise reference data is still lacking. Here, we employ the random phase approximation (RPA) on previously validated SCAN-rVV10 geometries to obtain interaction energies of state-of-the-art accuracy on the six high-symmetry stacking configurations ( Hhh , HhM , HhX , Rhh , RhM , and RhX ) of MX2 (M = Mo, W; X = S, Se) bilayers and compare them with the dispersion-corrected density-functional theory (DFT) functionals Perdew–Burke–Ernzerhof (PBE)+D3(BJ), PBE-rVV10L, and SCAN-rVV10. We identify SCAN-rVV10 as most reliable DFT variant with an average deviation of 1.2 meV/atom in relative energies from the RPA reference, and a root mean squared error of less than 2 meV/atom for interlayer interaction energies. We find interlayer distances obtained by PBE+D3(BJ) as being too short, with an impact on the electronic structure, resulting in the incorrect prediction of the band gap character in some cases. A further result of this work is the significant lowering of the interlayer energy and increasing of the interlayer distance in the high-energy stacking configurations. These stackings can be accessible via shear strain and promote exfoliation.


Introduction
Despite the absence of covalent or ionic bonds, interlayer interactions in stacked two-dimensional (2D) materials can significantly affect their electronic and optoelectronic properties.This is often referred to as proximity effect.The proximity effect can be particularly pronounced in transition metal dichalcogenides (TMDCs).For example, a semiconductormetal transition is observed when the number of layers in noble metal dichalcogenides increases from one to three or more [1,2].Even if the impact on the electronic structure seems to be marginal at first glance, it can have a qualitative effect on the properties.A striking example is the strongly enhanced photoluminescence in MoS 2 monolayers [3,4] as result of a tiny shift of the valence band near the Γ and K points, which causes the transition from being an indirect (bilayer to bulk) to being a direct (monolayer) 2D semiconductor.The effect has also been observed in other group 6 TMDCs [5][6][7].Recent progress in processing TMDCs down to transfer of single-layer flakes enabled the fabrication of homo-and heterobilayers with well-defined interlayer twist.The eclipsed (rhombohedral) configuration (R stacking, twist angle θ = 0 • , see figure 1) is slightly higher in energy than the dominant hexagonal configuration (H stacking, twist angle θ = 60 • , figure 1), which also constitutes the layered bulk structures.Small shear shifts of the layers with respect to one another result in three high-symmetry configurations each for the R and H stackings (figure 1), thus yielding a total of six high-symmetry stacking configurations [8].
In twisted TMDC bilayers a moiré pattern is observed.For small twist angles in the vicinity of 0 • (R stacking) and 60 • (H stacking), however, structure relaxation results in the formation of superlattices composed of high-symmetry stacking domains [9][10][11][12].Large domains of low-energy stacking and small nodes of high-energy stacking emerge, while so-called stacking solitons connect the nodes and separate the low-energy domains.The relative size of the domains is given by the relative energy of the high-symmetry stackings, delicately balanced with the strain energy necessary to form the solitons [13,14].The quantitative theoretical assessment of the domain energy and interlayer distance therefore is crucial to correctly model the domain reconstruction of twisted homoand heterobilayers.
Unfortunately, density-functional theory (DFT) as most popular method in the field lacks an intrinsically accurate treatment of London dispersion, which governs the interlayer interactions.Moreover, the performance of a posteriori London dispersion corrections varies strongly for different TMDCs [15].Therefore, a more rigorous, albeit computationally demanding, theory is required, at least for benchmark calculations of the high-symmetry stacks, which can then serve for validation purposes.
Here, we employ the random phase approximation (RPA) as benchmark method for the six high-symmetry configurations of all combinations of {Mo,W}{S,Se} 2 bilayers (figure 1).Then, we validate the popular London-dispersion-corrected DFT functionals SCAN-rVV10, Perdew-Burke-Ernzerhof (PBE)-rVV10L, and PBE+D3(BJ) for these systems.As the prediction of interlayer distance varies among the DFT variants, we study the impact of their geometries on the resulting electronic structure by comparing their band structures, calculated using the hybrid HSE06 functional, which is known to be a reliable method for this purpose.We identify the SCAN-rVV10 functional as the best-performing DFT variant, providing excellent relative energies and exfoliation energies, whereas PBE+D3(BJ) significantly overestimates the interlayer interaction and thus predicts too small interlayer distances, with consequences on the prediction of the band gap character.

Computational details
For all systems, the Vienna Ab Initio Simulation Package was employed [16][17][18][19].The projector augmented wave method [20] was used.The plane wave basis cutoff to the kinetic energy was set to 388 eV.We used a Gaussian smearing of 0.05 eV in the self-consistent field calculations, and a total energy convergence criterion of 10 −8 eV.The out-of-plane vacuum space for both monolayers and bilayers was set at a fixed distance of 18 Å.For the RPA reference calculations, we used the geometries optimized using the SCAN functional [21] augmented with the revised Vydrov-van Voorhis nonlocal correlation functional (rVV10) [22] to account for London dispersion.We have confirmed SCAN-rVV10 as being a state-ofthe-art technique for predicting interlayer distances in layered materials in alignment with experimental values earlier [15], and use it as reference method throughout this work.It has been previously shown that the geometries obtained from different London dispersion corrected DFT variants produce similar RPA relative energies [15], and for systems suggesting strong differences we carried out explicit checks which confirmed the reliability of the SCAN-rVV10 structures.For geometry optimizations, a Γ-centered k-grid of 12 × 12 × 1 (c = 1 is out-of-plane) was chosen.Full structural relaxations, incorporating the consideration of space group symmetries, were carried out using the atomic simulation environment [23] to optimize both lattice vectors and atomic positions until the forces and stresses were reduced below a specified threshold of 0.01 eV Å −1 .
RPA calculations [24][25][26][27][28][29] have been carried out on Kohn-Sham bands obtained using the PBE functional [30].In order to ensure Hartree-Fock convergence, a more refined set of Γ-centered k-grids with dimensions of 21 × 21 × 2 was utilized for the RPA Hartree-Fock (E EXX ) energy calculations, while for the RPA correlation energies (E C ) the same k-grids as those used in the geometry relaxations were used.
For better comparability, we define the d X−X int as the vertical separation of the inner chalcogen atoms (X) of separate layers, measured projected on the direction perpendicular to the layers (z-direction).For straight-forward comparison with the literature, distances between the metal layers are shown in the supporting information (figures S1 and S2).Energy differences (∆E) are calculated with reference to the most stable configuration for each studied system and method and is given per unit cell.We distinguish between interlayer interaction (E int ) and exfoliation energies (E exf ), where E int is defined as the energy gained by separating the fixed bilayer-optimized monolayers, while E exf accounts for the interaction between fully optimized bi-and monolayers: Here, E opt BL represents the total energy of the fully optimized bilayer, while E BL ML1 and E BL ML2 denote the total energies of monolayer 1 and monolayer 2 within the optimized bilayer lattice constant.
denote the total energies of the fully optimized monolayer 1 and monolayer 2, respectively.

Results
The principal target results of this work are the energy differences (∆E) between high-symmetry configurations (H h h , H M h , H X h , R h h , R M h , and R X h ), interlayer interaction energies (E int ), exfoliation energies (E exf ), and interlayer distances between the chalcogen atoms (d X−X int ) of {Mo,W}{S,Se} 2 bilayers.The ∆E metric revealed four low-energy (H h h , H X h , R M h and R X h ) and two high-energy configurations (H M h and R h h ).This differentiation is consistently observed across all bilayers and for all methods employed here.The highenergy configurations typically have a ∼0.6 Å larger interlayer distance compared to their low-energy counterparts.The extent of this difference is consistent across all the methods and all {Mo,W}{S,Se} 2 homo-and heterobilayers.Table 1 lists the relative energies ∆E between low-and high-energy configurations, where the average ∆E is 12.08 meV/atom at the RPA level.
Among the dispersion-corrected DFT methods, SCAN-rVV10 correlates closest with the RPA values, with an average ∆E of 10.85 meV/atom.Conversely, PBE-rVV10L underestimates this value, yielding 8.01 meV/atom, while PBE+D3(BJ) significantly overestimates it (17.01meV/atom).Among lowenergy configurations, RPA energy differences vary by not more than 1.5 meV/atom, except for the H X h configuration, which is up to 4 meV/atom less stable than the other low-energy configurations.All methods employed here yield similar trends on ∆E between the stacking configurations (figures 2(a)-(d) for MoS 2 -MoS 2 and MoS 2 -MoSe 2 , and figures S1 and S2 for the eight remaining systems).To elaborate, PBE-rVV10L, which tends to underestimate ∆E, predicts smaller energy differences between configurations belonging to the same group (low energy/low energy, highenergy/high-energy), whereas PBE+D3(BJ) consistently overrates those differences compared to the RPA reference.Notably, only SCAN-rVV10 gives values close to the RPA (tables S1-S3).
To discuss the performance of the different methods for assessing interlayer distances, we focus on the generally most stable low-energy configuration (H h h ) and the highest-energy configuration (R h h ).Again, the performance of the different methods varies strongly.SCAN-rVV10 predicts average values of 3.21 Å for H h h and 3.80 Å for R h h .PBE-rVV10 estimates slightly larger distances of 3.24 Å for H h h and 3.84 Å for R h h .PBE+D3(BJ), on the other hand, produces considerably smaller values, finding 2.99 Å for H h h and 3.59 Å for R h h .Using SCAN-rVV10 geometries as a reference, we see that the results obtained with PBE-rVV10L are close to those of SCAN-rVV10, while PBE+D3(BJ) consistently yields smaller interlayer distances.
Table 1.Energy differences (∆E), interlayer interaction energies (E int ), exfoliation energies (E exf ), and interlayer chalcogen-chalcogen distances (d X−X int ) of six high-symmetry configurations of 10 bilayer systems.Interlayer distances are computed using SCAN-rVV10, while all energies are obtained using the RPA on SCAN-rVV10 geometries.All energies are given in meV/atom, all distances in Å.  (H h h and H X h ), and high-energy nodes (H M h ) emerge, while close twists with respect to the R M h stack result in low-energy R M h domains and R X h nodes.The significantly different interlayer distance of the low-and high-energy domains results in a corrugation of the systems, and ∆E relates directly to the relative domain sizes [9,12,40].Therefore, accurately predicting ∆E between the high-symmetry stacking configurations is vital, as domains and nodes govern the electronic properties of these supercells.
As already indicated by the shorter interlayer distance, low-energy configurations result in stronger interactions (E int has larger negative values, E exf has larger positive values).We generally observe about 1.5 times stronger E int and E exf in low-energy configurations compared to their high-energy counterparts, as depicted in figures 2(c)-(f) for MoS 2 -MoS 2 and MoS 2 -MoSe 2 bilayers, and in figures S1 and S2 for remaining systems.For the MoS 2 -MoS 2 bilayer, RPA-calculated E int and E exf are predicted to be ∼30 meV/atom for low-energy and ∼20 meV/atom for high-energy stackings.Methods like SCAN-rVV10 and PBE-rVV10L provide estimates close to these values, whereas PBE+D3(BJ) significantly overestimates both E int and E exf (figure 3).This general pattern is consistent across all other tested systems.Figure 3 details both signed errors as boxplot and the root mean squared error for E int and E exf across all systems and configurations, showing lower errors for SCAN-rVV10 and PBE-rVV10L methods compared to PBE+D3(BJ).
Throughout all systems, no significant differences in E int and E exf are observed as long as the stacking configurations are identical.E int is approximately −30 meV/atom for low-energy and −20 meV/atom for high-energy stackings.As also shown in figures S1 and S2, the interlayer interaction is slightly stronger for systems with stronger polarizable selenium atoms (E Se−Se int < E S−Se int < E S−S int ) with minor energy differences of ∼1 meV/atom between the systems.
Exfoliation energies, however, show notable weakening in comparison to interlayer interaction  S2).In the case of the latter system, we even observe negative E exf for the H M h and R h h configurations, implying a spontaneous exfoliation.This might be the reason why ultrasonication, producing shear stress, is an efficient method for exfoliating 2D materials [41,42].To prevent potential errors associated with structural dependence, we confirmed the negative exfoliation energies of these configurations by using the PBE+D3(BJ) geometries for RPA energy calculations.Thus, we can conclude that high-energy stackings show significantly lower, in one case even repulsive, exfoliation energies.
In slightly twisted {Mo,W}{S,Se} 2 bilayers, both low-and high-energy stacking domains are present for symmetry reasons.Thus, the samples contain areas of high-energy stackings that are unstable and cannot be realized in the untwisted systems, and the electronic structure has contributions of the lowenergy bilayer stackings with indirect band gaps, nodes with direct band gap characteristics as in the monolayers, and strain solitons creating superlattices imposing Dirac cones or flat bands [9,12].These variations can be exploited to design materials with specific electronic properties that show local variation.
We now turn to the electronic structure of the investigated bilayers.Given the giant spin-orbit coupling (SOC) in their monolayers [43] and SOC annihilation only by bilayers with inversion symmetry (H stacking) [44], we calculate all electronic structures at the HSE06+SOC level.A detailed examination of the band structures of the studied systems in both low-and high-energy configurations reveals a pattern match observed in ∆E metrics, as shown in figure 4. The energy difference between K and Γ positions in the valence band (E val K − E val Γ ) of the studied systems in low-and high-energy configurations follows a similar trend.However, a closer look reveals a larger E val K − E val Γ difference for geometries optimized with PBE+D3(BJ).These geometries also produce a larger deviation in band gap values, and for some systems different band gap types are predicted (MoSe 2 -MoSe 2 , MoSe 2 -WSe 2 , and WS 2 -WS 2 in low energy configurations, see figures 4 and figures S7, S9 and S10 as well as tables S4 and S5).
Taking the HSE06+SOC electronic structures of the SCAN-rVV10 geometries as reference, we find that heterobilayers with different chalcogen atoms always, that is, in all low-and high-energy stacking configurations, have a direct band gap.{Mo,W}S 2 bilayers exhibit an indirect band gap except for the high-energy configurations of MoS 2 -WS 2 and WS 2 -WS 2 , while for {Mo,W}Se 2 bilayers no general trend is observed.
The MoSe 2 -MoSe 2 bilayer displays an indirect band gap in low-energy configurations and a direct band gap in high-energy configurations.Therefore, in a MoSe 2 -MoSe 2 twisted bilayer superlattice, it is expected that domains would have an indirect, but nodes a direct band gap, endowing the superlattice with spatially varying electronic properties.These results are consistent with the literature where particular systems have been reported, including MoS 2 bilayers [12] and various TMDC heterobilayers in the H h h configuration [45][46][47].

Conclusion
We present reference energies and interlayer distances for group 6 TMDC {Mo,W}{S,Se} 2 bilayers for their six high-symmetry stacking configurations.We endorse the SCAN-rVV10 method as a DFT variant giving very close results compared to our reference, the RPA, while the popular PBE+D3(BJ) method overestimates interaction energies, predicts too short interlayer distances, and as consequence is not reliable for predicting electronic properties based on the produced structures, in particular failing to assign the correct band gap character for some of the target systems.
Consistently, low-energy configurations show by ∼0.6 Å shorter interlayer distance compared to the high-energy counterparts, which is the origin of corrugation in the twisted bilayer systems with small twist angles with respect to the eclipsed (R) and hexagonal (H) stackings.The energy difference between low-and high-energy stackings is ∼10 meV/atom, which indicates that application of shear stress is a viable method to facilitate exfoliation.For heterobilayers containing different chalcogen atoms the exfoliation energy is significantly lowered due to the release of strain energy resulting from the lattice mismatch between the layers.Additionally, these systems generally exhibit direct band gaps, whereas systems containing only sulfur as chalcogen atoms generally have indirect band gaps.
We expect this work to be particularly useful for the validation and (re)parameterization of approximate methods, to allow for accurate calculations of complex group 6 TMDC structures including multiple twisted and untwisted stacks, which we expect to have intriguing optoelectronic properties due to domain formation and locally different electronic characteristics.

Figure 1 .
Figure 1.(a) Six high-symmetry stacking configurations of the 10 {Mo,W}{S,Se}2 homo-and heterobilayers.Different colors indicate separate layers.For homobilayers, the configurations R M h and R X h are equivalent.M and X atoms are displayed as large and small spheres, respectively.(b) Matrix showing studied bilayer systems.

Figure 3 .
Figure 3. Box plots summarizing the signed errors (SE) of E int and E exf with respect to the RPA reference for all six high-symmetry stacking configurations of the 10 {Mo,W}{S,Se}2 bilayer systems for the interaction energy (left) and the exfoliation energy (right).

Figure 4 .
Figure 4. Difference in energy between the Γ and K points in the valence band resulting from HSE06+SOC calculations performed on geometries optimized with SCAN-rVV10, PBE-rVV10L and PBE+D3(BJ).Homobilayers are shown in blue, heterobilayers with one type of chalcogen atoms in red, and with two types of chalcogen atoms in green; filled markers indicate indirect, empty markers direct band gaps.