Moire flat bands in twisted 2D hexagonal vdW material

Moire superlattices in twisted bilayer graphene (TBG) and its derived structures can host exotic correlated quantum phenomena because the narrow moire flat minibands in those systems effectively enhance the electron-electron interaction. Correlated phenomena are also observed in 2H-transitional metal dichalcogenides moire superlattices. However, the number of moire systems that have been explored in experiments are still very limited. Here we theoretically investigate a series of two-dimensional (2D) twisted bilayer hexagonal materials (TBHMs) beyond TBG at fixed angles of 7.34 and 67.34 degree with 22 2D van der Waals (vdW) layered materials that are commonly studied in experiments. First-principles calculations are employed to systemically study the moire minibands in these systems. We find that flat bands with narrow bandwidth generally exist in these systems. Some of the systems such as twisted bilayer In2Se3, InSe, GaSe, GaS and PtS2 even host ultra-flat bands with bandwidth less than 20 meV even for such large angles, which make them especially appealing for further experimental investigations. We further analysis the characters of moire flat bands and provides guidance for further exploration of 2D moire superlattices that could host strong electron correlations.


I. INTRODUCTION
Moiré superlattice (MSL) is a special type of 2D layered material, generated by stacking 2D vdW materials with a small lattice mismatch or with a twist angle, including graphene, hexagonal boron nitride (hBN), transition metal dichalcogenides (TMDs), various 2D magnets and superconductors [1]. Different from their parent 2D materials, MSLs with emerging global symmetry and periodicity exhibit fascinating quantum phenomena due to periodic moiré modulation of onsite potentials, interlayer coupling and intralayer atomic strain, such as the formation of second-generation Dirac cones [2], Hofstadter butterfly states [3] and shear solitons and topological point defects [4][5][6].
Except for the typical materials of graphene-based and TMDs, the twisted bilayer black phosphorus [73] and grey antimonene [74] are also studied in structural and electronic properties respectively by Guo et al. More recently, Liu et al. investigated antiferro-and ferroelectric bilayer In 2 Se 3 with large twist-angles [75] based on the first-principles calculations, where low-energy extremely flat band is found. Moreover, the studies of MSLs are further extended to the complex magnetic materials [76][77][78][79]. Balents et al. theoretically investigated the twisted bilayers of vdW magnets in the structures and phases [77], while Tong et al. considered the twisted bilayer 2D magnets CrX 3 (X=Br, I) from the magnetization textures aspect [76]. Furthermore, other moiré dimensionalities are also explored from quasi-one dimension [80] up to three dimension [81][82][83], which greatly extends the use of twistronic in multi-dimensional systems.
The highly tunable correlation and superconductivity properties of MSLs also make them appealing for future technology applications. A couple of the early attempts were made by Jarillo-Herrero et al. and Rickhaus et al., who made use of the gate-tunable correlated and superconductivity phase of magic angel TBG to fabricate Josephson junctions in singlecrystal nanostructures [84,85]. Moreover, Jarillo-Herrero et al. showed signatures of unconventional ferroelectricity in the bilayer graphene/boron nitride moiré system, which may lead to ultrafast, programmable and atomically thin carbon-based memory device applications [86].
The studies of the MSLs have provided powerful venue to explore the correlated physics and unconventional superconductivity [1,[87][88][89] as well as their applications in future technology. In the meanwhile, the emergence of thousands of new vdW layered materials [90] gives access to tremendous opportunities for the research of different types of MSLs. However, it remains unclear that whether moiré flat bands are generally exist in MSLs when the twist angle is small enough or the moiré periodic is large enough. Moreover, it is unclear what types of 2D materials are more susceptible to the formation of moire flat bands and how one could find moiré flat band systems with a relatively large twist angle or small system size, which may potentially give rise to stronger electron-electron correlation and/or higher transition temperature for unconventional superconductivity. To address these questions, we perform first-principle calculations to study systematically 22 twisted homo-bilayer superlattices constructed with 2D materials that are accessible in experiments. We constraint our calculations to systems with twist angles at 7.34 • and 67.34 • as we want to look for moiré flat-band systems with small system size. By analyzing the band structures, we find bands with significantly reduced bandwidths generally exist in these MSLs. Interestingly, we do find a few MSLs that can host ultra-flat bands with bandwidth less than 20 meV even for such a large twist angle and small system size. Together with the relatively strong Hubbard interaction, these systems are expected to host strong electron-electron correlations, which is extremely appealing for further experimental investigation. We further discuss the band characters of the parent 2D materials that lead to the formation of the ultra-flat bands and provide guidance for future exploration of numerical exotic strongly-correlated MSLs.

II. MODEL AND COMPUTATIONAL APPROACHES
The present calculations are done within density functional theory (DFT) using the Vienna ab initio software package (VASP) [91]. For these configurations of 2D twisted superlattice, the exchange correlation functionals of Perdew Burke and Ernzerhof (PBE) [92] are used, in conjunction with Tkatchenko-Scheffler (TS) [93,94] vdW corrections, which has been shown to give results well consistent with the experimental observations in our previous work on TMD-based MSLs [62]. An energy cutoff of 400 eV for the plane wave basis sets and the Γ-centered k-meshes of 1×1×1 are used for geometry optimization and electronic structure calculations (a 18 × 18 × 1 k-grid is used for 1 × 1 unit cell). A vacuum thickness larger than 15Å is used to avoid artificial interactions between periodic slab images. All atoms are fully relaxed with residual force per atom less than 0.01 eV/Å. Considering the computational cost of superlattice calculations, while the internal atomic positions are fully optimized, the lattice constant for moiré supercells is fixed to a value such that it corresponds to the experimental lattice constant for a 1x1 unit cell (See Table I for supercell lattice parameter for each MSL). For all calculations, due to the relativistic effect in heavy elements existing in most systems of TBHMs, the spin-orbit coupling (SOC) effect is considered while results without SOC are also calculated for comparison to estimate the effect of the SOC on the moiré flat bands.

A. Flat bands and electronic correlation in MSLs
We calculate the band structures of the 22 MSLs at the PBE+TS level with and without SOC. The results are summarized in Table I and the detailed band structures can be found in the SI. Except for the two twisted metallic systems (i.e., twisted bilayer PdTe 2 and PtTe 2 ), we find bands with narrow bandwidths commonly appear at the band edges in these MSL systems regardless of the variety of band structures in the original untwisted form. This suggests moiré engineering via creating MSL is quite effective in general in creating flat or narrow bands in 2D semiconductor systems. For the metal-  I. The summary for the calculations of 2D twisted materials including: materials index (i), structural formula, moiré supercell lattice constant, average interlayer distance (D), bandwidths at the VBM (W-VBM) without (left) and with SOC (right), bandwidths at the CBM (W-CBM) without (left) and with SOC (right), relative permittivity (ε), on-site coulomb repulsion energy (U) and the projected orbitals at the VBM and CBM (only the top contributions listed), and the cohesive energy defined as the energy difference between the total energy of a bilayer with that of two well-separated monolayer in the primitive 1 × 1 cell.
Bi2Se3 lic system, the bands near the Fermi level are highly entangled such that it is difficult to identify a well-defined flat band (see Fig. S23). We list the bandwidth of the flat bands appearing at the band edges in these systems in Table I. Not all systems have a value there either because for some systems the bands at the band edges are still quite dispersive (bandwidth W > 200 meV) or there is no well-defined isolated flat bands at the twist angles we study in this work. Furthermore, we calculate the unfolded band of PtSe 2 to further elaborate the moiré flat band. As shown in Fig. S25 of SI, the unfolded band structure (Fig. S25 (a)) in the twisted bilayer is drastically different from its counterpart (Fig. S25 (b)) in the untwisted structure. In particular, one can visualize a section of isolated flat band appearing at the VBM near the Γ point in the Brillouin zone of the primitive cell. This feature corresponds exactly to the flat bands that we calculated in the supercell Brillouin zone. The unfolded band structure reveals that the flat bands are originated from the linear recombination of the VBM states near the Γ point. Similar features of the flat bands have been observed in twisted bilayer graphene [127,128] and WS 2 /WSe 2 moiré superlattices [129]. From Table I, it is clear that it is easier to find flat bands at the VBM than at the CBM. For the flat bands at the VBM (CBM), the bandwidth ranges from 0.9 (42) meV to larger than 100 meV in our calculations with SOC. Nevertheless, it is surprising to find that the isolated band at the band edges in some twisted bilayer systems (such as twisted In 2 Se 3 , GaSe, GaS, InSe, PtS 2 ) can be so flat that its bandwidth is less than 20 meV, even the twist angles we study here are relatively large.
The existence of flat bands in these MSLs indicates the kinetic energy scale of the electron states in these bands is significantly quenched and the electron-electron interaction and correlation may become important. To further evaluate the correlation effects, we estimate the on-site coulomb repulsion energy U in these system as e 2 /(4πεε 0 a), where e is electron charge, ε 0 is the vacuum permittivity, ε is relative permittivity and a is the effective linear dimension of each site (here we take the length scale of the moiré pattern). Though a combination of relative permittivity ε from the data of other literature and the lattice parameter of the moiré unite cell a, we estimate U and the values are depicted in Table I. Then, we compare the energy scale of the bandwidths W in VBM/CBM and estimated Hubbard U of different MSLs by plotting them in Fig. 1b (for MSLs with the same materials, only the smallest value of bandwidth is shown). In the region on the top left of Fig. 1b (colored in pink), the Hubbard U is larger than the bandwidth W, which indicates correlation effects will be important; while in bottom right region, the bandwidth W is larger and the electron correlation effect will be weak. Fig. 1b shows that the twisted compounds of GaS, GaSe, InSe, In 2 Se 3 , PtS 2 all have narrow bandwidths and relative strong Hubbard interactions and they are expected to host strong electron correlations even at such a large twist angle of 7.34 • . It is noteworthy that twisted arsenene should have a U value near that of antimonene, comparable or even more larger to its relatively small bandwidth of 24 meV, although we couldn't find the value of its ε in literature and didn't list it in the figure. The other systems near the diagonal line, such as VBM of PtSe 2 , Bi 2 Se 3 , Sb 2 Te 3 and MoSe 2 and CBM of GaS, ZrS 2 , HfS 2 , SnS 2 and SnSe 2 , show comparable U and W values, indicating electron correlations are also important in these systems. For those systems locate at the lower right region in the diagram, such as MoTe 2 and ZrS 2 , the electron correlation may be less important. We need point out that we only compare the relevant energy scale for twisted systems at a relative large twist angle here, the locations of data points in this diagram (i.e., Fig. 1b) will change as the twist angle decreases. The twisted systems appear to have weaker electron correlations here could become strongly correlated systems at smaller angles.

B. Characters of the flat bands in MSLs
To better understand the formation of flat bands in 2D MSLs, we further analysis the characters of the flat bands. To this end, we conduct the calculation of the projected band structures and the results are shown in Figs. S1-S22 in the SI and the major orbital components for the states in the flat bands are summarized in Table I. The results reveal that for the flat band systems with smaller dispersion in VBM including the twisted materials of Bi 2 Se 3 , In 2 Se 3 , Sb 2 Te 3 , Bi 2 Te 3 , GaS, GaSe, InSe, PtS 2 and PtSe 2 , the valence band edges are associated with the p z orbital of the chalcogen atom predominantly. Apart from these materials, the characters are also predominantly p z orbitals for both of the twisted arsenene and antimonene in VBM. This is not surprising as we only study MSLs with a relatively large twist angle in this work, the atomic reconstruction has a minor effect. The moiré modulation of the band structure is dominated by the modulation of interlayer coupling, instead of the atomic reconstruction as reported in the work by Crommie and others [72] for MSLs with a much larger system size. Under such circumstance, as the electronic states with p z orbital are very sensitive to the interlayer coupling, they can be significantly modified by the moiré potentials created by the modulated interlayer coupling and turned into flat-band states.
We take the twisted Bi 2 Se 3 at 7.34 • (Fig. 2a) as a typical example to further elaborate the role of different atomic orbitals in the formation of flat bands. The calculated band structure and the corresponding projected band structure of twisted Bi 2 Se 3 are given in Fig. 2b and 2c. From these figures, it is clear that the states at the valence band edges are associated with Se p z and a fraction of Bi s states. Whereas, the states at the conduction band edges are predominantly Bi p z states. As the Se p z orbitals are located at the outermost Se atomic layer, with charge density extended towards the stacking interface (See lower panels Fig. 2d), they are sensitive to the modulation of the interlayer coupling. Thus those states at the VBM, with relatively large contribution from Se p z orbitals, are significantly altered by moiré interlayer potentials, forming a flat band with a small bandwidth of 21 meV. The charge density distribution of the states in this flat bands is very localized in the real space as shown in Fig. 2d. On the other hand, the Bi p z orbitals locate at the inner atomic layers and their wavefunction barely extend towards to interface region (see lower panel of Fig. 2e). Therefore, those states at the conduction band edge mainly contributed by the Bi p z orbitals are much less sensitive to the modulation of interlayer couplings and the bands at the CBM are still very dispersive. The corresponding charge density distribution of those states is delocalized over the whole moiré cell as shown in Fig. 2e.
Next, we look at those systems that host flat bands with extremely small bandwidths even at a relatively large twist angle. Here, we investigate twisted bilayer In 2 Se 3 and GaSe at 67.34 • as two typical examples. As shown in Figs. 3b and 3f, the flat bands appear at the VBM with a extremely narrow bandwidth W of 0.8 and 5 meV for twisted bilayer In 2 Se 3 and GaSe at 67.34 • , respectively, in the calculations without SOC. When including SOC, the band structures of both systems have significant modifications: for the In 2 Se 3 system, the band gap between the top and the lower flat bands increases and additional flat bands appear at higher energies; for the GaSe system, a Rashba type of splitting is introduced in the top flat band. Nevertheless, the bandwidths in both systems remain small after including SOC (0.9 meV for the In 2 Se 3 system and 6 meV for the GaSe system). Similar to what we discussed above for the case of twisted Bi 2 Se 3 , these ultraflat bands states are mainly derived from the p z orbital of the atoms in the outermost atomic layer (i.e., Se atoms in these cases), as shown in Fig. S24. Moreover, although the chemical composition and atomic structures are different, both systems share similar features in the band structure of the untwisted 1x1 form: that is, a relatively flat plateau at the VBM. This is also the case for the other ultra-flat band systems such as GaS and InSe. Such band plateau, when folded in the moiré supercell, naturally appears as flat bands. The moiré potential due to the modulated interlayer coupling further introduces band gaps between these flat bands at the supercell Bouillon Zone (BZ) boundary and confines those electronic states, leading to isolated ultra-flat bands. This is actually similar to the trilayer graphene/hBN moiré superlattice system [36,37], where the states at the band plateau at the VBM are confined by the moiré potential formed by the graphene/hBN superlattice. As the band plateau region is more extended in the BZ in these systems, a much smaller moiré length scale is sufficient for the formation of ultra-flat bands. It is expected that the moiré heterostructures of these 2D layers could also host flat bands at the VBM.
Twisted bilayer PtS 2 is also an interest system. Different from the systems discussed above, although untwisted prinstine bilayer PtS 2 does not has a flat band plateau at the band edges, its twisted form at 7.34 • also host a ultra-flat band with a bandwidth of 15 meV in the calculation without SOC as shown Fig. 4b. The band structure doesn't change much when including SOC (see Fig. S16 in SI). Although the bandwidth is larger than the cases discussed above, it is still considerable small and comparable to the bandwidth of twist bilayer graphene at the magic angle of 1.05 • , even at such a large twist angle. A noticeable feature in this twisted system is that the separation between flat bands is relatively large, indicating the strength of the interlayer moiré potential is relatively large. The relatively large moiré potential is likely to be related to the relatively strong interlayer hybridization of the S p z states. As shown in Fig. 4c and 4d, the S p z states near the VBM in the untwisted bilayer have a considerable large energy splitting of about 1.5 eV compared with those in the monolayer, which even shifts the VBM from the S p x -p y states in the monolayer to the S p z states in the untwisted bilayer. Another system that hosts relatively strong moiré potential is twisted bilayer buckled arsenene. As shown in Fig. 4g and 4h, the As p z states near the VBM in the untwisted bilayer arsenene also have a relatively large energy splitting of about 1.0 eV, shifting the VBM from the As p x -p y states to the As p z states. The twisted bilayer arsenene also hosts a large band gap between flat bands as shown in Fig. 4f.
Finally, we discuss the flat bands occurred at the CBM. Generally it is much harder to form flat bands at the CBM. Nevertheless, we found a few exceptions, such as twisted bilayer GaS, 1T-ZrS 2 , HfS 2 and SnS 2 . The band structures of twisted bilayer 1T-HfS 2 and GaS at 7.34 • without including SOC effect are shown in Figs. 5b and 5e, respectively. The figures clearly show that isolated flat bands with similar shape appear at the bottom of the conduction bands in the two systems, with a bandwidth of 55 (52) meV for 1T-HfS 2 (GaS). The flat bands don't change much when the SOC effect is included (see Figs. S18 and S5 in the SI). A similar feature for these systems is that the bottom of the conduction bands of the pristine 1x1 bilayer locates at the M point in the BZ (see Fig. 5c and 5f. So the flat bands at the CBM of these systems are formed by the M point states. On the other hand, when the bottom of the conduction bands locate at the Gamma point, the bands at the CBM remain highly dispersive at large twist angles (see the previously discussed twisted bilayer Bi 2 Se 3 and GaSe for example.).

IV. CONCLUSIONS
To summarize, we report a systematic investigation of band structures to explore flat bands in a series of 2D TBHMs with van der Waals beyond TBG using first-principles methods. Two configurations of moiré superlattice sharing supercells with the same size at twisted angles of 7.34 • and 67.34 • are considered. Our calculations show that the twisted compounds of Bi 2 Se 3 , In 2 Se 3 , Sb 2 Te 3 , Bi 2 Te 3 , GaS, GaSe, InSe, PtS 2 , PtSe 2 , arsenene and antimonene host flat bands at the VBM with the small bandwidths (less than 100 meV). Meanwhile, flat bands also emerge at the CBM for the twisted bilayer materials of GaS, HfS 2 , HfSe 2 , ZrS 2 , SnS 2 and SnSe 2 . Those systems that host flat bands at relatively large twist angles generally have the following characters: (1) The states at the band edges are mainly contributed by the outermost atoms of the layered materials and from those atomic orbitals where the charge density extends towards the stacking interface, such as p z and d z 2 orbitals; (2) the band curvatures at around the band edges are relatively large or even relatively flat band plateaus appearing at the band edge in the Brillouin zone of the primitive cell. Then, by estimating the Hubbard interaction U, we find that the twisted compounds of Bi 2 Se 3 , In 2 Se 3 , GaS, GaSe, InSe and PtS 2 exhibit large Hubbard U over bandwidth W ratios, indicating electron correlation may  be relatively strong in those systems. Note that our study is limited to the systems with relatively large twist angles. When the twist angle further decreases and the system size becomes larger, flat bands may also appear in other systems, although for those systems the bands at the band edges are still quite dispersive at the angles we studied here. Nevertheless, the ultra-flat band systems we discuss here provide ideal candidates for the study of strong correlation effects at large twist angles. Finally, we discuss the characters of systems that host flat bands at large angles and provide a guideline for future exploration of novel 2D MSLs that host flat bands and potential strong electron correlations.

V. ACKNOWLEDGEMENT
This work is supported by the Key-Area Research and Development Program of Guangdong Province of China (Grants No.2020B0101340001). The computational resource is provided by the Platform for Data-Driven Computational Materials Discovery of the Songshan Lake laboratory.