Atomic relaxation and electronic structure in twisted bilayer MoS2 with rotation angle of 5.09 degrees

It is now well established theoretically and experimentally that a moir\'e pattern, due to a rotation of two atomic layers with respect to each other, creates low-energy flat bands. First discovered in twisted bilayer graphene, these new electronic states are at the origin of strong electronic correlations and even of unconventional superconductivity. Twisted bilayers (tb) of transition metal dichalcogenides (TMDs) also exhibit flat bands around their semiconductor gap at small rotation angles. In this paper, we present a DFT study to analyze the effect of the atomic relaxation on the low-energy bands of tb-MoS2 with a rotation angle of 5.09 degrees. We show that in-plane atomic relaxation is not essential here, while out-of-plane relaxation dominates the electronic structure. We propose a simple and efficient atomic model to predict this relaxation.


Introduction
The broad family of transition metal dichalcogenides (TMDs) [1][2][3] offers the possibility to stack two layers with a small angle of rotation θ to each other, thus forming a moiré pattern superstructure. These twisted bilayers have given rise to numerous experimental [4][5][6][7][8][9][10][11][12][13] and theoretical  studies to understand electronic states that are confined by a moiré pattern in semiconductor materials. Many of these studies analyze the interlayer distances, the possible atomic relaxation, the transition from a direct band gap in the monolayer system to an indirect band gap in bilayer systems, and more generally the effect of interlayer coupling in these twisted bilayer 2D systems at various rotation angles θ. For small values of θ, the emergence of flat bands has been established from first-principles density functional theory (DFT) calculations [21,31] and tightbinding (TB) calculations [28][29][30] in twisted bilayer MoS 2 (tb-MoS 2 ), and observed in a 3 • twisted bilayer WSe 2 sample by using scanning tunneling spectroscopy [13]. It has been also shown numerically [26] that Lithium intercalation in tb-MoS 2 increases interlayer coupling and thus promotes flat bands around the gap. There is also experimental evidence that moiré patterns may give rise to confined states due to the mismatch of the lattice parameters in MoS 2 -WSe 2 heterobilayers [12].
Most theoretical investigations of the electronic structure of bilayer MoS 2 are density-functional theory (DFT) studies [4-9, 14-21, 31, 36-39] with eventually a Wannier wave function analysis [15]. To provide systematic analysis as a function of the rotation angle θ, in particular for small angles, i.e., very large moiré pattern cells for which DFT calculations are not feasible, several TB mod-els, based on Slater-Koster (SK) parameters [40], have been proposed for monolayer MoS 2 [41][42][43][44][45] and multi-layer MoS 2 [14,15,29,30,41,43]. Following these efforts, we have proposed [28] a SK-TB set of parameters for non-relaxed structures, i.e., rigidly twisted bilayers, that match correctly the DFT bands around the gap of tb-MoS 2 with rotation angles θ > 7 • . This SK-TB model, with the same parameters, is then used for smaller angles in order to describe the states confined by the moiré pattern. For θ < θ C 5 • , the valence band with the highest energy is separated from the other valence states by a minigap of a few meV. In addition, the width of this band decreases as θ decreases so that the average velocity of these electronic states reaches 0 for θ 2 • such that almost flat bands emerge at these angles. This is reminiscent of the vanishing of the velocity at certain "magic" rotation angles in twisted bilayer graphene [46][47][48]. However, in bilayer MoS 2 it arises for an interval of angles and not a set of specific values. Other minigaps and flat bands are also found in the conduction band. The confined states that are closest to the gap are localized in the AA stacking regions of the moiré pattern, as in twisted bilayer graphene. However, for small angles, it has been shown [21,31,33] that atomic relaxation strongly modifies these low-energy bands, in particular their approximate degeneracy. A better understanding of atomic relaxation and its effects on the electronic structure is therefore essential not only for very small angles but also for larger ones (θ 5 • ).
In this paper, we present a DFT study of tb-MoS 2 with a rotation angle of θ = 5.09 • . This angle is approximately the angle below which the highest energy valence band (just below the gap) is isolated from the rest of the valence states by a minigap for non-relaxed structures. We show arXiv:2303.07958v2 [cond-mat.mes-hall] 12 Jul 2023 that, unlike the case of very small angles [21,31,33], atomic relaxation amounts to essentially out-of-plane atomic displacements along the z-direction that can be simulated simply as a function of the atomic positions in the moiré cell. This simple atomic model allows to understand the origin of the modifications of the electronic structure induced by the relaxation.

Structure and DFT calculations
The construction of the commensurable twisted bilayer we study is explained in detail in Refs. [28,49]. It corresponds to (n, m) = (6, 7) with 762 atoms per unit cell. Starting from an AA stacked bilayer (where Mo atoms of a layer lie above a Mo atom of the other layer, and S atoms of a layer lie above an S atom of the other layer), one layer (layer "+") is rotated with respect to the other layer (layer "−") by the angle θ = 5.09 • around an axis containing two Mo atoms. Thus, AA stacking regions are located at the corner of the moiré cell ( Fig. 1(bottom)). BA' stacking regions (where S− atoms lie above a Mo+, and Mo− (S+) do not lie above any atom of layer + (layer −)), and AB' stacking regions (regions where Mo− lie above an S+, and S− (Mo+) do not lie above any atom of layer + (layer −)) are located at 1/3 and 2/3 of the long diagonal of the moiré cell, respectively.
The DFT calculations were carried out with the ABINIT software [50][51][52]. We have checked previously [28,49] that LDA [53] and GGA [54] + Van der Waals exchange-correlation functionals yield very similar results, so all the results presented here are based on LDA calculations, which require less computation time for large systems. While the results presented here may not be quantitatively very accurate, they are sufficient to qualitatively show the importance of out-of-plane atomic displacements on the electronic structure of the low-energy flat bands. Likewise, the exact value of the band gap may not be very precise in LDA, but the evolution of this gap and the existence of minigaps are significant. The Brillouin zone was sampled by a k-point mesh in reciprocal space within the Monkhorst-Pack scheme [55]. One k-point is used for atomic relaxation and a 2×2 k-grid for the self-consistency procedure of the electronic structure calculation. The kinetic energy cutoff was chosen to be 408 eV. We checked that a 3×3 k-grid for the self-consistency procedure and an energy cutoff of 544.22 eV yield very similar bands for the relaxed structure. The structural optimization of atomic positions is done by using the Broyden-Fletcher-Goldfarb-Shanno minimization. A vacuum region of 10Å was inserted between the MoS 2 bilayers to avoid spurious interactions between periodic images. In our calculations, the spin-orbit coupling (SOC) is not taken into account in order to reduce the calculation time. SOC is important for TMDs, as it introduces some band splittings close to the gap [21,29,35,44,45], however it will not change the existence or not of low-energy flat bands in tb-MoS 2 , so it is not essential for the present study.

Atomic relaxation
The DFT-relaxed atomic structure of tb-MoS 2 with rotation angle θ = 5.09 • is shown in Fig. 1. Remarkably, the in-plane displacement τ i of each atom i with respect to the rigidly twisted structure is rather small. Indeed, the average τ i is 0.04 ± 0.02Å, 0.03 ± 0.02Å, and 0.05 ± 0.02Å, for atoms S ±ext , Mo ± , and S ±int , respectively. Such displacements are almost not visible in Fig. 1(bottom), and they have little effect on the electronic structure (see next section). This result shows that the strong in-plane displacements obtained for the smallest angles [21,31,33] are not too important for θ 5 • . However, it is interesting to note that these small displacements are precursors to the larger displacements and shear solitons obtained for the smallest angles [21], as shown in Fig. 2. These in-plane displacements tend to reduce the AA stacking regions with respect to the AB stacking regions to minimize the energy [21,31,33].

S -ext S +ext
Mo -Mo + S -int S +int On the other hand, the displacements along the zdirection are important. For each atomic layer, the mean values of the z-coordinate and the corresponding standard deviation are given in Table 1. As expected given the interplane distances in the simple stacking cases (see Ref. [38] and Table 1), the distance between layers is greater in the AA stacking regions than in AB' (BA') stacking regions. Inspired by the work of Koshino et al. [57] for twisted bilayer graphene, we propose the following atomic model ("z-mod" model) where xy in-plane atomic coordinates are those of the rigidly twisted bilayer and the modulation of the atomic z-coordinates are calculated by where r are the non-relaxed positions of the atoms in the rigidly twisted moiré cell, and G i are the six vectors of the reciprocal lattice that define the first Brillouin zone of the moiré pattern. For each atomic layer, i.e., for the six Mo and S atomic layers of tb-MoS 2 , z AA and z AB are the z values of the DFT-relaxed structure at r = 0 (AA stacking region) and r at AB' or BA' stacking regions (Table 1). Figure 1 shows that DFT-relaxed positions and z-mod positions fit well together. Table 1. Atomic positions: Averagez value per layer and the corresponding standard deviation σz, for layer + (z > 0) and layer − (z < 0) in non-relaxed, relaxed, and z-mod tb-MoS2 structures. The atoms of the z-mod structure have the same in-plane xy coordinates as the non-relaxed structure and their z coordinates are calculated from Eq. (1) with zAA and zAB (i.e., in AB' and BA' stacking regions) of the DFT-relaxed structure. For each layer δz = |zAA − zAB|. The four last lines correspond to simple AA and AB' (BA') bilayer stackings for whichz = zAA and zAB, respectively. All calculations were performed using the LDA pseudopotential, except z GGA AA and z GGA AB that are calculated with the GGA PBE pseudopotential [54] including Van der Waals corrections (D2) [52,56]. The latter calculations show that the displacements obtained by the LDA are qualitatively correct. Thez values for the two layers are symmetric with respect to z = 0. Distances are in A.

Structure
S±ext

Electronic band dispersion
The DFT band dispersions are shown for the DFT-relaxed bilayer and the rigidly twisted bilayer (non-relaxed bilayer) in Fig. 3(a) and 3(b), respectively. For the nonrelaxed bilayer the interlayer distance is that obtained for simple AA stacking, like in our previous calculations [28]. In the non-relaxed tb-MoS 2 , the minimum of the conduction band is at K like for monolayer MoS 2 , which is no longer the case after atomic relaxation. Fig. 3(c) shows a zoom of the highest-energy valence bands. The relaxation does not change the valence band maximum energy at K. However it leads to significant modifications of the bands close to the gap. For the non-relaxed structure [28], when θ < θ c , the highest-energy band is non-degenerate and isolated from the other valence bands by a minigap. For θ = 5.09 • , this minigap is equal to ∼ 0.5 meV (Fig. 3(c)). For the relaxed structure there is no isolated single valance band, but 2 bands which cross in K while remaining linear in k (around K). This result is similar to that of Naik et al. [21,31] and Vitale et al. [33], obtained by using a multiscale approach, a pair potential for relaxation, and DFT or TB calculations for the band dispersion. There are nevertheless small differences. For instance, in our calculation, the bandwidth for these first two valence bands is 84 meV with respect to ∼ 120 meV in Ref. [31]. Unlike that previous calculation, these two bands are isolated from the rest of the valence bands by a minigap of ∼ 20 meV. To test the validity of the structure with the z-modulation only (Eq. (1)), bands around the gap of DFTrelaxed and z-mod structures are compared in Fig. 4. The gaps of the two structures are slightly different (difference of 36 meV), which may be due to the difference between the average z values per atomic layer ( Table 1). The conduction bands are nevertheless almost the same and valence bands are very similar. In particular, the two valence bands closest to the gap that are characteristic of the atomic relaxation effect are well reproduced. Therefore, the simplified relaxation model given by formula (1) is sufficient to account for the low-energy flat bands due to a moiré pattern at not too small rotation angles.

Conclusion
We have performed a DFT study of the atomic relaxation and the electronic band dispersion of twisted bilayer MoS 2 with a rotation angle equal to 5.09 • . Contrary to what has been observed for very small angles (typically less than ∼ 2 • ) [21,31,33], the in-plane atomic displacements with respect to a rigidly twisted bilayer are very small, and they have almost no effect on band dispersion. However, the out-of-plane displacements are large and significantly alter the low-energy bands around the gap. These atomic displacements can be modeled by a simple formula that depends only on the interlayer distances in the AA and AB stacking regions.
The reduction of bandwidth and related emergence of flat bands identifies weakly doped MoS 2 bilayers as good candidates for the observation of strong correlation effects. For a complete theoretical study of electronic correlations in these complex systems, it is important to take into account the atomic relaxation. We offer here a simple outof-plane atomic displacements model for not too small rotation angles, typically a few degrees. Preliminary investigations (not shown here) indicate that for determining the electronic structure of low-energy bands, the Slater-Koster tight-binding models [28,29], that are efficient for rigidly twisted bilayers, would require a further adjustment of the parameters in order to be applicable to the z-modulated structures.

Authors contributions
S. Venkateswarlu, A. Misssaoui, and G. Trambly de Laissardière performed the DFT calculations and the numerical analysis. S. Venkateswarlu, A. Honecker, and G. Trambly de Laissardière wrote and revised the manuscript. All authors discussed the results and approved the final version of the manuscript.