The role of stacking on the electronic structure of MoSe2 at small twist angles

We consider two high symmetry stackings AA and AB and examine the changes induced in the electronic structure by considering small angles of rotation of 3.48° from both these stackings. In both cases we largely recover the low energy electronic structure of the untwisted limit. We additionally find flat bands emerging above the dispersing bands. Surprisingly, while the rotation from the AA end leads to one flat band above the highest occupied band at Γ, one finds two flat bands emerging for small rotations from the AB end. Examining the real space localisation of the flat bands allows us to discuss the origin of the flat bands in terms of quantum well states and qualitatively understand the dependence of the number of flat bands found on the twist angle.


Introduction
Advances in experimental techniques have allowed us to assemble van der Waals heterostructures from various layered two dimensional materials.The weak van der Waals interactions that exists between the layers has allowed for an enormous freedom in the construction of the heterostructures.The layers could be displaced or rotated with respect to each other.The most surprising aspect that has emerged is that these assembled heterostructures show unique characteristics which are not found in their component materials.This has led to the huge interest in these materials [1,2,3,4,5,6,7,8].The twisted bilayers were first realised in graphene [1,2,3,4] and were accompanied by the formation of flat bands and strongly correlated phases.Later on, this has been achieved in layered Mo and W based transition metal dichalcogenides [5,6,8].At certain twist angles, these materials exhibit several unusual properties like the formation of flat bands, zero resistive behaviour, electric field modulated metalinsulator transitions etc. which are not found in their untwisted limit [5,6].
Rotating one layer with respect to the other breaks the periodicity of the primitive cell.As a result, a large moiré cell appears.This leads to a smaller Brillouin zone.Consequently, bands get folded to a significantly small Brillouin zone, emerging as weakly dispersing flat bands.This has been an explanation offered for the formation of flat bands and the associated strong correlation phenomena in these twisted bilayers [15,16].Besides, there have been theoretical reports that shows formation of flat bands at certain limits using continuum model [17,18,19].However, there is no discussion given regarding the mechanism responsible for its formation.
Each layer is coupled to the other by weak van der Waals forces implying weak interactions between the layers.Consequently, twisting one layer with respect to the other is not expected to introduce large changes in the electronic structure compared to the untwisted limit.Therefore, we have labelled the bands of the moiré cell with the kpoints of the untwisted unit cell.In order to do this, we have projected the eigenfunctions of the moiré bands onto those of the untwisted limit.In case of a perfect supercell this would give us the primitive cell bandstructure.However, for twisted bilayers the unit cells are not perfect supercells.Hence, one would expect that projecting the moiré cell bandstructure onto the untwisted limit, would provide us with a measure of the extent of perturbation.As the perturbation is expected to be weak, we should largely recover the band structure of the untwisted limit.
While the unfolding of the electronic structure for large unit cells has been used widely in the context of alloys [9,10], doped unit cells [9,11], our earlier work [12,13,14] established that the electronic structure in these van der Waals systems should be viewed as a perturbation to the untwisted limit.In a recent work [13], we have examined the electronic structure of two large and similar sized moiré cells for twist angles 19.03 • and 3.48 • for MoSe 2 .By projecting out the eigenfunctions of the bands of the large moiré cell onto the primitive cell eigenfunctions for a twist angle of 19.03 • , dispersing bands were found similar to the untwisted limit with the electronic structure being hardly perturbed.However, carrying out the same analysis for a twist angle of 3.48 • , we found that in addition to the dispersing bands which can be identified with the untwisted limit, a weakly dispersing band emerges above these bands.This has been associated with the presence of regions of strong perturbation in the moire cell.
However the effect of rotations upon starting with different stackings has not been adequately examined.We explore this by considering twist angles near 0 • and 60 • for AA stacked MoSe 2 .Here 60 • rotation from the AA stacking leads to another high symmetry stacking, AB stacking.Carrying out a similar analysis, we found a pair of split off bands at the top of the valence band.This split-off band has a vanishingly small bandwidth and is even more flat than what we had for 3.48 • .This has been associated with stronger localisation found here.In each of the instances, the localization of the wavefunction associated with the flat band is in a region in which one has an enhanced interlayer hopping interaction.This leads to the antibonding states to be pushed to higher energies compared to the states in the region in which the interaction strength is less.We can associate the flat band with the antibonding states arising from these interactions.The regions characterized by enhanced interlayer interactions can be considered to be associated with a potential well of width L. This mapping allows us to discuss the emergence of the flat bands for a twist angle of 56.52 • with the solutions of a potential well far separated from its neighbour, allowing us to have two bound states (the flat bands).However, for a twist angle of 3.48 • , the solutions are for two potential wells of similar depth and width L, separated by a small region.In the limit where the separation between them vanishes, we can approximate the well width to 2L.As the energy of the solutions varies as 1 (2L) 2 , the shallower energies here allows us to have only one bound state (one flat band).

Methodology
Ab initio electronic structure calculations have been carried out using density functional theory (DFT) as implemented in the Vienna Ab initio Simulation Package (VASP) [20,21].Projector augmented plane wave potentials (PAW) have been used to describe the ion-electron interactions [22].The generalised gradient approximation [23,24] has been used for the exchange-correlation functional.The lattice constants of MoSe 2 has been kept fixed at the experimental value of 3.289 Å. Internal positions of the atoms have been optimised through a total energy minimization.A vacuum of 20 Å has been introduced along the z-direction for all cases to minimise the interaction between the periodic images which are inevitably present because of the periodic implementation of DFT that we use.Van der Waals interactions have been introduced between the two layers using the DFT-D2 method of Grimme [25].The electronic structures for twisted bilayers have been solved self-consistently at Γ point.Starting from AA stacking, we rotated the upper layer with respect to the lower one anticlockwise by an angle θ.In order to examine the twist induced perturbation, the electronic structures calculated for the large moiré cell have been projected onto the high-symmetry directions of the unrotated primitive cell using the unfolding technique [26].
A primitive cell wave vector k pr gets folded to a supercell wave vector K sc if there exists a supercell reciprocal lattice translation vector G which satisfies the condition k pr = K sc + G. Now using Blöch's theorem, any energy eigenstate corresponding to the supercell wave vector K sc could be written as where u Kscm (r) is the Blöch function and has the periodicity of the supercell.This can be expanded in terms of the reciprocal lattice vectors G of the supercell as So Eqn. 1 can be written as We determined these Fourier coefficients using the open source interface Wavetrans [27] for VASP.Then, the contribution of the m-th eigenstate of supercell k-point K sc at the primitive cell k-point k pr is calculated as where g is the translation vector of the reciprocal lattice of the primitive cell.Only those G coefficients which are equal to g + k pr − K sc will contribute.

Results and Discussions
We have considered the primitive cell of MoSe 2 which has AA stacking.While the Mo atoms form a hexagonal lattice, the presence of the Se atoms leads to C 3 rotational symmetry instead of C 6 rotational symmetry as found in graphene.This makes the rotation about 0 • and 60 • inequivalent for TMDCs.For instance, a rotation of 60 • starting from AA stacking leads to another high symmetry stacking AB, using the notation of identifying primitive cells of Ref. [28].Therefore, we examined the twisted bilayers for twist angles near 0 • and 60 • separately.For this purpose, we have considered a twist angle of 3.48 • as well as its conjugate angle 56.52 • as 56.52 • =(60 • -3.48 • ).This could also be viewed as twisting by an angle of 3.48 • clockwise starting from the high symmetry stacking AB.Unit cells of both twist angles are shown in Fig. 1(a) and (b).
In each of the cases, the length of each lattice vector of the unit cell is 54.14 Å, with the moiré cell containing 1626 atoms.For these small twist angles near zero and sixty degrees, one finds large patches of various high symmetry stackings within the moiré cell.These patches are not found for large twist angles.Even if at a particular point, an atom on atom condition is satisfied, the nearest neighbour atoms are substantially rotated away.This therefore does not allow one to have large regions of any high symmetry stackings for large twist angles.This is the primary distinction between small and large angles of rotation.
Considering the twist angle of 3.48 • (Fig. 1(a)), one finds patches of AA stacking at the corners marked by orange semicircular arcs and patches of AB ′ stacking in the middle marked by the green circles.The side view and the top view of these primitive cell stackings are shown in Fig. 1(b).Depending on the relative position of the atoms of the upper layer with respect to those of the lower layer, each stacking has different interlayer separations.For instance, one finds that in AA stacking, Se atoms of the upper layer are sitting exactly on top of the Se atoms of the lower layer.Hence, the electrons on the Se atoms of the upper layer experience larger Coulomb repulsions from the electrons on the Se atoms of the lower layer.Consequently, the vertical interlayer separation increases in the region of AA stacking and is found to be 3.74 Å. Considering the primitive cell for AA stacking one finds that the equilibrium interlayer separation is similar and is found to be 3.79 Å.In AB ′ stacking the upper layer is staggered with respect to the lower layer.This is also clear from the top view of the primitive cell AB ′ stacking (Fig. 1(b)).The Se atoms of the upper layer sit on top of the centre of the hexagon of the lower layer.Consequently, the electrons on the Se atoms of the upper layer feel a smaller Coulomb repulsion than in the previous case and the layer separation is found to be reduced to 3.19 Å in the region of AB ′ stacking.Considering the primitive cell associated with this stacking, one finds the interlayer separation to be 3.20 Å similar to the patches discussed above.These different stackings can be associated with the The twist angles near sixty degrees show large patches of three other high symmetry stackings making it inequivalent to the twist angles near zero degrees.The unit cell for twist angle 56.52 • is shown in Fig. 1(c).Here, at the corners, one finds large regions of AB stacking marked with black solid semicircles.In the middle, patches of AA' (marked with red dashed cirle) and patches of A ′ B stackings (marked with green solid circle) are found.The side view and the top view of the primitive cells of these three stackings are shown in Fig. 1(d).Depending on the orientation of the upper layer with respect to the lower layer, one finds that the that the interlayer separation is substantially different in each case.Considering the primitive cells in each of these three stackings, the equilibrium interlayer separations are found to be 3.25 Å, 3.20 Å and 3.75 Å for AB, AA'(2H) and A'B stackings respectively.This variation in the interlayer separation has also been reflected in the relaxed moiré cell of twist angle 56.52 • .The interlayer separation reaches a maximum value of 3.74 Å in the region of A ′ B stacking and reaches a minimum value of 3.20 Å in the region of AA ′ (2H) stacking.
The electronic structure has been calculated for the moiré cell of twist angle 3.48 • as well as 56.52 • .However, as a consequence of the large moiré cells, the electronic structure gets folded into very small Brillouin zone resulting in flat dispersionless bands irrespective of the angle of twist.The folded electronic structure plotted in the small Brillouin zone associated with the moiré along the Γ to K sc direction, for twist angles 3.48 • , 56.52 • as well as 19.03 • have been shown in Fig. 2. In each of the cases, we find almost dispersionless bands in the entire energy window, which is contrary to experiments.We have discussed this aspect in Ref. [12,13,14], explaining why the interesting phenomena are seen in the small angle regime.
As the moiré cell electronic structure can be viewed as a weak perturbation to the electronic structure of the untwisted limit, each eigenfunction of the moiré cell has been projected onto the eigenfunction of the primitive cell along the Γ to K direction, as described in the methodology section.The projected band structure along the Γ to K of the primitive cell is shown in Fig. 3 for twist angle 3.48 • .The colour bar on the right side represents the percentage of the weight of the energy eigenstate at that primitive cell kpoint.In contrast to the band dispersions expected for a moiré supercell, where one expects almost dispersionless bands, here, we found a set of dispersing bands.These dispersing bands are similar to what is seen in the untwisted limit.The untwisted limit band structure is superposed in Fig. 3 for comparision.One finds that the projected eigenfunctions which have large weights, follow the band dispersions of the untwisted limit.However in addition to that, a weakly dispersing band emerges at the top of the valence band at Γ.This split-off band is well separated from the valence bands below and has a bandwidth of 19 meV.The highest occupied band at K is found to be energetically much lower than that at Γ with the separation being 84 meV.Including spin orbit interactions reduces this Γ-K separation to 5 meV with the valence band maximum remaining at Γ.An expanded view of the near K region of the band structure close to the fermi energy, calculated including spin-orbit interactions is shown in the inset of Fig. 3.This indicates that there are flat bands at K.
In a similar manner, the projected band structure for twist angle 56.52 • along the Γ to K direction of the primitive cell has been plotted in Fig. 4 with the untwisted limit band structure superposed.Here also, we found a set of dispersing bands associated with the valence bands, like the untwisted limit.In addition to those, a pair of split off bands are found above the dispersing valence bands.These split-off bands have a vanishingly small bandwidth.Here also, the highest occupied valence band at K is energetically much lower than that at Γ with the Γ-K separation being 74 meV.Inclusion of spin orbit interactions reduces this Γ-K separation to 13 meV with Γ being the valence band maxima.An expanded view of the near K region of the band dispersions close to the fermi energy, calculated including spin-orbit interactions has been shown in Fig. 4. One finds flat bands emerging at K also.However, these split-off bands are not found for large twist angles [13] and the band structure is similar to the untwisted limit.It immediately raises the question what is the origin of the perturbation in the electronic structure of small twist angles near 0 • and 60 • .
In these materials, the layers are coupled to each other by weak van der Waals forces.In a previous work [30], we have examined the origin of the changes in the electronic structure with layers.This has been carried out by mapping the ab-initio electronic structure of a bilayer primitive cell onto a tight binding model.Within this model, the interlayer hopping interactions have switched off.This gives us the monolayer's bandstructure.This implies that the changes in the electronic structure is due to the interlayer hopping interactions only.Therefore, in order to understand the origin of the perturbation in case of the small twist angles near zero and sixty degrees, one needs to examine the interlayer hopping interactions in these moire cells.As these scale with the interlayer Se-Se distance, we plotted the spatial profile of the interlayer Se-Se bond lengths which are less than 3.8 Å in each of the cases.This is shown in 5 (a) and (b) for twist angles 3.48 • and 56.52 • respectively.The regions where these shortest bond lengths are concentrated is around the region of AB ′ stacking for twist angle 3.48 • as shown in Fig. 5 (a) and around the region of 2H stackings for twist angle 56.52 • as shown in Fig. 5 (b).In these regions, the Se atoms of one layer have more number of Se neighbours from the other layer with these bond lengths.Hence these regions give rise to a strong perturbation to the electronic structure of the untwisted limit in each case.
Examining the charge density associated with the highest occupied band, the untwisted limit has the charge density almost uniformly distributed on both layers.However, at these small twist angles the spatial distribution of the interlayer hopping interaction strengths modifies the charge density localization.The unit cell as well as the localization of the charge density for a twist angle of 3.48 • is shown in Fig. 5(a).It follows the spatial profile of the enhanced interlayer hopping interactions and is localized in the same region.A similar localization trend is also seen for the twist angle of 56.52 • .
The hopping interaction between the orbitals on the two layers leads to a set of bonding and antibonding states.One can view the charge density plotted in Figs.6(a) and (b) as being associated with the antibonding state emerging from the region with enhanced interaction between the two layers.The stronger interactions also pushes the antibonding state to higher energies compared to the state in the region with weaker interactions.One can therefore associate a quantum well of spatial width L in the region associated with each of the two lobes of the charge density in Fig. 5(a), with a small barrier in between.A similar potential well of width L can be associated with the region where the charge density is localized in Fig. 6(b).As the interlayer separations in the regions of the well should be similar, one expects potential wells with similar depths, say V 0 , in both cases.Associating the flat bands that we find with the bound states of the potential well problem, we can try to qualitatively understand the existence of flat bands at 3.48 • and 56.52 • .Let the well of width L have two bound states for the potential V 0 .Approximating the barrier width to zero in the scenario where the two potential wells merge, their width is now 2L.As the energy of each state varies as 1 (2L) 2 , the lowest energy bound state is now shallower, and so we have only one flat band here.This qualitative analysis helps us understand the origin of one flat band for 3.48 • , and the existence of two flat bands for 56.52 • , with the energetically deeper flat band being associated with the excited state of the potential well problem.This is evident from the charge density shown in Fig. 7. Examining the dispersional width of the band, one finds it to be 19 meV for a twist angle of 3.48 • , while it is vanishingly small for 56.52 • .In order to understand this, we have plotted the area of localisation of the charge density of this band in each case.This has been determined by calculating the integrated charge density around each Mo atom by considering a sphere of a radius of 1.262 Å (half of Mo-Se bond length).We have marked the area of the charge density localisation up to where the integrated charge density over a sphere around each Mo atom falls to 10 % of its maximum value for both the twist angles 3.48 • and 56.52 • in Fig. 8 (a) and (b) respectively.
It is evident from Fig. 8 that the area of charge density localisation of the split off In conclusion, we have examined the electronic structure of small rotations about AA and AB stacking.This leads to a low energy electronic structure involving dispersive bands, quite similar to the untwisted limit.In addition we see almost flat bands just above these dispersing bands.We find that the number of flat bands forming, and their localization can be discussed in terms of the spatial extent of the perturbing potential.

Figure 1 .
Figure 1.Unit cell for twist angle (a) 3.48 • and (c) 56.52 • .The Mo/Se atoms are shown by violet/green spheres.Regions with high symmetry stackings have been indicated.The primitive cell stacking AA and AB ′

Figure 3 .
Figure 3. Electronic structure (without SOC) for twist angle 3.48 • projected along the primitive cell Γ to K direction, with the colour representing the weight at that k-point.The green line shows the primitive cell band structure for AA stacking at an interlayer separation of 3.49 Å.An expanded view of the highest occupied bands near K including SOC is shown as an inset.

Figure 4 .
Figure 4. Electronic structure for twist angles 56.52 • projected along the primitive cell Γ to K direction, with the colour representing the weight at that k-point.The green line shows the primitive cell band structure for AA stacking at an interlayer separation of 3.49 Å.An expanded view of the highest occupied bands near K including SOC is shown as an inset.

Figure 6 .
Figure 6.Charge density associated with the highest occupied split off band at Γ for twist angle (a) 3.48 • and (b) 56.52 •

Figure 7 .igure 8 .
Figure 7. Charge density associated with the second highest split off band at Γ for twist angle 56.52 •