Impact of atomic reconstruction on optical spectra of twisted TMD homobilayers

Twisted bilayers of transition metal dichalcogenides (TMDs) have revealed a rich exciton landscape including hybrid excitons and spatially trapped moir\'e excitons that dominate the optical response of the material. Recent studies have shown that in the low-twist-angle regime, the lattice undergoes a significant relaxation in order to minimize local stacking energies. Here, large domains of low energy stacking configurations emerge, deforming the crystal lattices via strain and consequently impacting the electronic band structure. However, so far the direct impact of atomic reconstruction on the exciton energy landscape and the optical properties has not been well understood. Here, we apply a microscopic and material-specific approach and predict a significant change in the potential depth for moir\'e excitons in a reconstructed lattice, with the most drastic change occurring in naturally stacked TMD homobilayers. We show the appearance of multiple flat bands and a significant change in the position of trapping sites compared to the rigid lattice. Most importantly, we predict a multi-peak structure emerging in optical absorption of WSe$_2$ homobilayers - in contrast to the single peak that dominates the rigid lattice. This finding can be exploited as an unambiguous signature of atomic reconstruction in optical spectra of moir\'e excitons in naturally stacked twisted homobilayers.


I. INTRODUCTION
The recent advances in the field of two-dimensional materials have opened up a new platform to study many-particle physics in nanomaterials [1].The field of twistronics is of particular interest, where the relative twist angle between two vertically stacked 2D materials can be used as an additional external knob to tune the material properties [2][3][4].Consequently, the introduction of a twist angle in these layered systems can lead to intriguing many-particle phenomena, such as unconventional superconductivity and correlated phases [5,6].Transition metal dichalcogenides (TMDs) have emerged as a significant and noteworthy subclass of two-dimensional materials, garnering substantial interest in recent years [7,8].Here, strongly bound electronhole pairs, known as excitons, form due to the twodimensional confinement and the reduced screening of the Coulomb interaction.They dominate the material's optical response even at room temperature [9][10][11][12].Furthermore, long-lived interlayer excitons can form, which exhibit an out-of-plane dipole moment due to the spatial separation of electrons and holes [13][14][15][16][17][18][19][20].Additionally, the large wave function overlap can lead to a significant carrier tunneling between the layers, giving rise to hybrid states of intra-and interlayer excitons [10,19,[21][22][23][24].
When introducing a twist angle between two vertically stacked TMDs, a periodic superlattice emerges.This leads to a moiré pattern, i.e. a periodic superlattice potential that can trap excitons at specific high-symmetry sites [21,[25][26][27][28][29][30][31].Furthermore, the emerging moiré potential can be tuned by changing the twist angle which has been shown to have a large impact on the optical * joakim.hagel@chalmers.seresponse of these materials [25,26,29,[32][33][34].Interestingly, recent studies have revealed that when the twist angle is small, i.e. the period of the moiré pattern is large, the two monolayer crystal lattices do no longer remain rigid, but instead undergo a deformation to reduce the total adhesion energy [35][36][37][38][39][40][41][42][43].Here, some energetically favorable high-symmetry sites grow to form large domains, separated by typically narrow domain walls (cf.Fig. 1).It has been predicted that this atomic reconstruction should have a major impact on the electronic band structure and fundamentally change the potential landscape [38,[44][45][46][47].However, it is still not thoroughly understood how atomic reconstruction can be identified in the optical response of moiré excitons.
In this work, we investigate photoluminescence (PL) and absorption spectra in atomically reconstructed TMD bilayers.Our approach is based on the ab initio derived continuum model for the deformation of the rigid lattice derived in Refs.[45,46] combined with an equation-ofmotion approach for the optical response of moiré excitons [19,21,30,31].This allows us to obtain microscopic insights into the changes in the moiré exciton energy landscape and the optical response in the presence of atomic reconstruction.The deformation of the moiré supercell into large periodic domains induces strain that can result in deeper potentials for moiré excitons.As a direct consequence, we find more trapped states than in the rigid lattice and the stacking region in which excitons are trapped can differ from the rigid lattice system.Furthermore, we predict that the KK interlayer exciton becomes lower in energy than the intralayer KK state due to the strain-induced potentials at low twist angles.The lowest lying momentum-dark KΛ exciton is also predicted to undergo a larger red-shift than in the rigid lattice, visible in PL spectra.Most importantly, we find a qualitative change in the absorption spectrum in naturally stacked bilayers reflecting the bright KK exci-tons, where multiple peaks emerge as a consequence of the strain-induced potentials -in stark contrast to the single peak that is present in the rigid lattice.Consequently, we predict an unambiguous signature of atomic reconstruction in the optical response of moiré excitons in naturally stacked twisted TMD homobilayers.

II. MOIR É POTENTIAL IN ATOMICALLY RECONSTRUCTED TMD BILAYERS
To gain access to the moiré exciton energy landscape in atomically reconstructed bilayers, we first describe the change in geometry in the superlattice.Following the approach in Refs.[45,46] we set up a functional for the stacking energy, taking into account the strain energy from continuum theory and a parameterized equation for the adhesion energy between the layers, where the model parameters have been fitted to data from density functional theory in Refs.[45,47].The functional can then be turned into an optimization problem by expanding the displacement vectors u l (r) as a Fourier expansion u l (r) = n u l n e ignr and minimizing the integral with respect to the Fourier coefficients u l n .Here, l is the layer index, g n are the reciprocal vectors of the moiré lattice and r is the real space coordinate in the supercell.For further details, see the supplemental part (SI) of the paper [48] (see also references [49][50][51][52][53][54][55][56][57] therein).
After obtaining the displacement vectors, we study the geometrical change in the supercell.Here, we find that the energetically more favorable H h h stacking grows into hexagonal domains (Kagome pattern) at small twist angles, separated by a thin boundary wall -in good agreement with experimental studies [35] and previous theoretical works [45].This reconstruction is illustrated in Fig. 1a showing the interlayer distance for H-type stacked WSe 2 homobilayers at the twist angle θ = 0.6 • .Each high-symmetry stacking is reflected by a specific interlayer distance, which can be computed from DFT [58] and thus the interlayer distance can be used as an indicator for a change in the local atomic registry.The blue domains denote the H h h regions, which are significantly larger than in the rigid lattice, cf.Fig. 1.b.The other high-symmetry stackings (H X h and H M h ) consequently shrink in size and form thin boundary networks (cf.Fig. 1.a).We focus here on the examplary H-type stacked WSe 2 homobilayer, but the same approach holds for all vertically stacked TMDs.In R-type stacked homobilayers, we instead obtain large triangular patterns as a result of the atomic reconstruction, which yields a different potential landscape (see SI for further details).
The impact of atomic reconstruction on the exciton energy landscape does not only stem from the geometrical change in the superlattice, but also from the strain that accompanies the lattice deformation.Assuming small displacements, the displacement vector u l (r) can be associated to the linear strain tensor ε ij = 1 2 (u i,j + u j,i ), where u i,j = 1 2 (∂ i u j + ∂ j u i ) and i(j) = (x, y).Due to its D3h symmetry, we can decompose the dominating strain contributions to the moiré potential into two parts (see SI) [59][60][61]: scalar potential composed of inhomogeneous uniaxial strain in both x-and y-direction [62,63] and inhomogeneous vector gauge potential known as piezo potential [45,47,64].Details about the strain potential model are provided in the SI.
In addition to the strain-induced potentials, there is a stacking-dependent alignment shift of the two monolayer band structures, originating from charge polarization between the two layers, that already exists in the nonreconstructed system [19,30].In a homobilayer, these only affect the energy of interlayer excitons [19].Furthermore, it is only present in R-type stacked bilayers due to the lack of inversion symmetry, which allows for local charge transfers between the layers [27].The remaining interlayer coupling mechanism contributing to the total moiré potential is given by the interlayer hybridization.For some electronic states, in particular in the Λ valley in the conduction band, the large wave function overlap between the layers leads to a strong mixing of states in both layers, which can be described as a mixing of interlayer and intralayer excitons into hybrid exciton states [19,21].In H-type stacking, the strongest hybridization occurs around the H h h stacking due its shorter interlayer distance (cf.Fig. 1.a-b).To sum up, we have four different contributions to the moiré potential: 1) scalar strain potential S ξ Lg , 2) piezo potential P ξ Lg , 3) stacking-dependent alignment potential V ξ Lg and 4) tunneling T ξ LL ′ g .The first three act as a periodic renormalization of exciton energies, i.e. they are diagonal with respect to the layer index, and can thus be modeled as an external potential acting on free excitons M ξ Lg = S ξ Lg + P ξ Lg + V ξ Lg , where L = (l e , l h ) is a compound layer index and ξ = (ξ e , ξ h ) is the exciton valley index.
In Fig. 1.c-d, we show the resulting band edge variation for the conduction band (Fig. 1.c) and the valance band (Fig. 1.d) in the bottom layer at θ = 0.6 • .Similarly, the band edge variation for the top layer is shown in Fig. 1.e-f, where the conduction band is given by Fig. 1.e and the valance band is given by Fig. 1.e at θ = 0.6 • .Here, the only potentials that are contributing are the scalar strain potential and the piezo potential.In general, homobilayers mostly reconstruct due to atomic rotation [37], which is associated with the piezo potential and will be concentrated at the H M h /H X h sites.However, there is also a significant scalar strain potential close to the boundary regions of the H h h domains which also significantly contribute to the band edge variation.The combination of scalar potential and piezo potential in turn gives the tilted star-shape of the total band edge variation.Comparing the variation for different layers (Fig. 1.e to Fig. 1.d ) we find that the minima of the conduction band and the maxima of the valance band occur at the same sites, thus yielding a very efficient band gap renormalization for interlayer excitons.In contrast, comparing Fig. 1.c and Fig. 1.d we find that the minima for the conduction band and the maxima for the valance band occur at different sites in the supercell.This means that there is a certain driving force for charge separation.However, the binding energies of the intralayer excitons are very strong (∼ -140 meV for KK) which in turn means that one would need a very deep potential in order to have a complete charge separation.At angles θ > 0.5 • , the formation of excitons is still favored in comparison to the separated unbounded electron/holes.For much smaller supercells, however, the formation of in-plane charge transfer excitons can occur as reported in [38].
The effective exciton potential as shown in Fig. 1.g gives the calculated potential depth for KK (blue) and KΛ (yellow) excitons as a function of twist angle, both for intra-(dashed) and interlayer (solid) excitons.Since the piezo potential in H-type stacking is the same for both the valence and the conduction band in both layers, the only contributing factor is the scalar strain potential.The color map of this effective potential is given in the supplemental part of the paper.Below θ = 2 • we find a significant decrease of the KK interlayer exciton potential (solid blue) due to the efficient band gap renormalization stemming from the scalar strain potential.The KK intralayer (dashed blue) potential has a much shallower potential due to the conduction and valance band shifting in the same direction, although with different rates.Comparing to the effective KΛ (solid yellow) potential, it only reaches half the potential depth of the KK.This can be understood from different orbital compositions of the valleys, which directly impact the change in the orbital overlap with scalar strain [65].Note that these effective exciton potentials are all in stark contrast to the rigid lattice where they do not occur at all for H-type stacking.For R-type stacking, we predict in general much weaker band edge variations, having only a maximum variation of ∼40 meV.This difference comes primarily from the lack of scalar strain in R-type stacking, were the domain walls are formed purely from shear strain -in good qualitative agreement with the experimentally observed values [37].In turn, this means that the effective intralayer exciton potential will be close to 0 meV.There will, however, be a significant addition to the interlayer potential in the form of the piezo potential (see SI).

III. MOIR É EXCITON ENERGY LANDSCAPE IN RECONSTRUCTED TMD BILAYERS
Having a microscopic access to the potential landscape of moiré excitons, we can now set up the complete moiré exciton Hamilton operator.Here, we start from a decoupled monolayer basis in electron-hole picture and add the moiré potential as periodic modifications to these energies [21,30,31].Transforming to an exciton basis, the Hamiltonian then reads where E ξ LQ is the binding energy of the decoupled monolayer as obtained from the generalized Wannier equation [66] and Q is the center-of-mass momentum.The effective masses for electrons and holes are obtained from Ref. [67].Furthermore, X ( †) denotes the annihilation (creation) operator of excitons.Moreover, M ξ Lg is the combined layer-diagonal moiré potential component and T ξ LL ′ g is the tunneling matrix element.The latter gives rise to the hybridization of intra-and interlayer excitons, which corresponds to the off-diagonal component of the moiré potential [58].
The periodicity of the superlattice is taken into account by applying the zone-folding scheme [21,30,31], where we restrict the summation over Q to the first mini-Brillouin zone (mBZ), where states with larger momentum are folded back via mBZ lattice vectors g giving rise to the formation of mini-subbands.The final moiré exciton energies are then obtained by first expanding the exciton operators into a moiré exciton basis , where Y ( †) are the moiré exciton annihilation (creation) operators, η is the new band index, and C ξη * Lg (Q) are the mixing coefficients that mix different sub-bands and layer configurations.Applying this transformation to Eq. 1 we obtain the moiré exciton eigenvalue equation Here, E ξ ηQ are the final moiré exciton energies, which can be obtained by solving Eq. 2 numerically.
In Fig. 2, we show the calculated moiré exciton band structure for H-type stacked hBN-encapsulated WSe 2 bilayers (for R-type stacking see SI) with a reconstructed lattice (Fig. 2.a) and with a rigid lattice (Fig. 2.b).We find a large red-shift in the reconstructed lattice, where the lowest sub-band is located ∼36 meV below the rigid lattice energies.Interestingly, due to the efficient moiré potential for the interlayer exciton, the nature of the lowest lying state has changed from an intralayer exciton to an interlayer exciton.While in the rigid lattice (cf.Fig. 2.b) most bands have a significant curvature, we find multiple ultra-flat bands in the reconstructed lattice (cf.Fig. 2.a), both for the interlayer excitons (black) and the intralayer excitons (blue).The spacing between the first two bands is about 15 meV, which is significantly larger than the small spacing (of about 3 meV) in the rigid lattice indicating a much stronger spatial confinement.This is a direct consequence of the additional strain potentials in the reconstructed lattice.Note that also the hybridization of the higher-lying bands is strongly suppressed as there is a lack of band crossing due to their flatness in the reconstructed lattice, cf. the color gradient in Fig. 2.a-b denoting the degree of hybridization.We need to go to higher energies to see some hybridization in the reconstructed case (red color around 10 meV in Fig. 2.a).
The moiré exciton wave function in real space |Ψ(r)| 2 for the lowest-lying states are shown in Fig. 2c-d.We find that in the reconstructed lattice the wave function is localized close to the edges of the H h h -domains of the of the moiré unit cell (cf.Fig. 2.c), reflecting the potential minima from the strain-induced potential.In contrast, in the rigid lattice with no additional strain potential (cf.Fig. 2.d) there is a much larger probability distribution centered at the H h h stacking configuration (cf.Fig. 1.b).Here, the moiré potential only consists of the tunneling and therefore the exciton localization is determined only by the stacking with the strongest tunneling [19].When increasing the twist angle, the reconstructed lattice starts to revert back to the rigid lattice, making the tunneling the dominating contribution of the moiré potential.This change between different regimes implies the possibility of tuning the localization of moiré excitons, which should have important consequences for exciton transport in these materials [68,69].

IV. IMPACT OF RECONSTRUCTION ON OPTICAL SPECTRA
Having access to the moiré exciton energy landscape, we can now determine the optical response of TMD bilayers in presence of atomic reconstruction.In WSe 2 homobilayers, the lowest-lying state in R-type stacking is the strongly hybridized KΛ exciton, while in H-type stacking KΛ and K ′ Λ ′ excitons are the lowest ones.Note that the latter two are degenerate in H-type stacking, but they have a reversed layer configuration, as a consequence of the reversed spin-orbit coupling the layers [19][20][21].We first focus on PL spectra, where these dark states are expected to be visible via indirect phonon sidebands [70].Here, the exciton can not directly emit a photon due to the large momentum mismatch between the valleys.Instead, it decays via a two-step process, first scattering to a virtual bright state within the light cone via phonon emission followed by the emission of a photon [70].
Figure 3.a-b shows the PL spectra for WSe 2 with Htype (cf.Fig. 3.a) and R-type stacking (cf.Fig. 3.b) at 4K at different twist angles.The reference frame for the energies are the same as in Fig. 2.a-b.We observe sidebands from momentum-dark states, where the smaller (larger) peak stems from optical (acoustic) phonon modes [19,21].The energetic position of the KΛ(K ′ Λ ′ ) phonon replicas will primarily depend on the interlayer hybridization, stemming from the strong interlayer tunneling of electrons around the Λ(Λ ′ ) valley.This gives rise to a twist-angle-dependent dehybridization, in turn yielding a blue shift when increasing the twist angle (gray peaks in Fig. 3), i.e. the larger the twist-angle the smaller is the wave function overlap [21,29].Since this component of the moiré potential is already present in the rigid lattice we find a very similar behavior in both the rigid and the reconstructed lattice when changing the twist angle.There is, however, a clearly larger red-shift for the reconstructed lattice at lower twist angles due to the increasing strain-induced potentials.This is especially noticeable for R-type stacking where the piezo potential contributes to the effective exciton potential.Consequently, the change in the peak position closely resembles the change in the effective exciton potential for KΛ (SI Fig. 3).In Fig. 3.b we can also see that the peak position of the reconstructed lattice is converging towards the peak position of the rigid lattice when increasing the twist angle, indicating that the lattice becomes more rigid at larger angles.So far we have been focusing on photoluminescence that is dominated by phonon sidebands of the energetically lowest dark KΛ excitons.These states are strongly hybridized [19,21], and are thus subject to a pronounced moiré potential already in the rigid lattice.The situation is qualitatively different for bright KK excitons that are accessible in absorption spectra.Here, the reconstructed lattice gives rise to a deep strain-induced potential, whereas the rigid lattice only experiences a small moiré potential due to the weak periodic tunneling strength around the K-valley.In Fig. 2.a we have already demonstrated the significant change occurring in the band structure of KK excitons, where multiple flatbands have emerged in the reconstructed lattice for both the intralayer and the interlayer excitons.
In Fig. 4.a-b, we have calculated the absorption spectrum of H-type stacked WSe 2 bilayers as a continuous function of the twist angle for both the rigid lattice (Fig. 4.a) and the reconstructed lattice (Fig. 4.b).Here, we immediately observe a drastic change in the spectrum, where below θ ≈ 1.0 • multiple peaks appear in a reconstructed lattice -in stark contrast to the single-peak spectra in the rigid lattice that are mostly angle-independent.Furthermore, at smaller angles we find a noticeable splitting around ∼15 meV between the two brightest peaks in the reconstructed lattice -in good agreement with experimental observations [44].We also find a small splitting in the rigid lattice due to the tunneling of holes, which consequently traps the lowest state at H h h (cf.Fig. 2.d).This state, however, is very close to the next non-trapped state (cf.Fig. 2.b), and thus the trapping is not very efficient.In the reconstructed lattice (cf.Fig. 4.b), we instead have a rich optical response, consisting of multiple peaks.Here, the lowest branches stem from excitons trapped in the deep strain potential pockets (cf.Fig. 2.a/c).Consequently, these peaks undergo a non-linear blue-shift as we increase the twist angle, following the non-linear twist angle dependence of the strain-induced moiré potential (cf.Fig. 1.g).Furthermore, in Fig. 2.a, we find that the lowest lying KK states are of pure interlayer character (black lines in Fig. 2.a).These are not visible in absorption spectra due to their by orders of magnitude smaller overlap between the electron and the hole wave function compared to intralayer excitons [71].When increasing the twist angle, the number of peaks is reduced until it converges to the rigid lattice case around θ = 1.2 • that is characterized by only one resonance.
In R-type stacking (shown in the SI), we do not observe the same effect.Here, the strain-induced potential for intralayer KK excitons is negligible and the impact on the spectra only stems from the weak tunneling of carriers.We predict this drastic change in optical absorption due to atomic reconstruction to be similar for all naturally stacked TMD homobilayers.Note that the multi-peak structure is also well-known from TMD heterostructures, such as MoSe 2 -WSe 2 , already in the rigid lattice case due to deep moiré potentials [30].Here, the impact of lattice reconstruction is only of quantitative nature with changes for the peak separation and amplitude.The situation is qualitatively different in naturally stacked TMD homobilayers, where the strain-induced potential in the reconstructed lattice is the dominant mechanism and leads to the characteristic multi-peak structure for the bright A exciton.The latter can thus be considered as an unambiguous optical signature for the presence of lattice reconstruction.

V. CONCLUSION
Overall, our work provides new microscopic insights into many-particle mechanisms dominating the optical response in twisted TMD homobilayers.The atomic reconstruction has a large impact on the moiré exciton energy landscape where strain-induced potentials lead to significant red-shifts and multiple new flat bands.Furthermore, we have calculated moiré exciton wave functions and predict that excitons are trapped at different sites in the reconstructed lattice reflecting the maxima of the strain potentials.Finally, we find a significant change in optical spectra of naturally stacked TMD homobilayers due to reconstruction.In particular, we predict the emergence of a multi-peak structure in absorption spectra -in stark contrast to the single-peak dominated spectra in a rigid lattice.Our work can help guide future experiments in the growing field of 2D material twistronics.

VI. ACKNOWLEDGMENTS
In the main part of the paper the total potential landscape and its corresponding KK moiré exciton band structure is discussed for H-type stacking.Here, we show strain fields for H-type stacking and the effective moiré exciton potential stemming from the strain.In addition, we also discuss the corresponding KΛ band structure.

A. Strain field
In Fig. 1.a-b we show the calculated strain field for the shear strain in the form of the maximum shear strain value γ max = u max − u min , where u max,min is given by Here, ε ij is the linear strain tensor (see section III B for further details).At lower angles (cf.Fig. 1.a) we can see that the magnitude of the shear strain increases while the H h h region also grows in size.At larger twist angles (cf.Fig. 1.b) we can see that the magnitude of strain decreases and instead forms a more triangular pattern around H X h , in very good agreement with previous works [1].Furthermore, in Fig. 1.c-f we show the calculated scalar strain for different angles in both the top layer and the button layer.Here, we predict that at smaller angles (cf.Fig. 1.c/e) a pronounced scalar strain will be present along the edges of the H h h domains with a range of ± 2%.Note that the shear component of the strain is still the dominate component, but that the scalar strain adds a significant contribution to strain profile.At larger angles (cf.Fig. 1.d/f) we find that the scalar strain decreases in magnitude to range of approximately ± 0.4%, in good agreement with previous theoretical predictions [1].

B. Effective moiré exciton potential
The total band edge variation shown in Fig. 1.c-f in the main part of the paper is the result from the combined scalar strain and piezo potential.In Fig. 2.a-d we show the effective exciton potential that arises from these band edge variations for the KK exciton.Here, the only contributing component is the scalar strain since the piezo potential is a local rigid shift of the entire band structure.First thing to note is that the potential minimum is localized around the edges of the H h h domains with varying sign, stemming from the strong scalar strain close arXiv:2308.14633v2[cond-mat.mes-hall]4 Mar 2024 to the edges.Furthermore, the interlayer exciton potential Fig. 2.b/d has a much deeper potential than the intralayer exciton potential Fig. 2.a/c.This comes as a consequence of the scalar strain shifting the valance and conduction band in opposite direction for the different layers due to one layer exhibiting compressive strain and the other tensile strain.In turn, this yields a very efficient moiré potential for the interlayer excitons which shifts them below the intralayer excitons in the case of KK excitons (Fig. 2 in main part).The intralayer excitons also exhibits a significant moiré potential, although not as deep as the interlayer potential.This stems from the difference in strain gauge factors for the bands as can be seen from Table II.

C. KΛ exciton band structure
In Fig. 3.a-b we show the calculated moiré exciton band structure for the low-lying KΛ exciton in H-type stacked WSe 2 in the reconstructed lattice (cf.Fig. 3.a) and the rigid lattice (cf.Fig. 3.b).Note that the reference point for the energies is given with respect to the A exciton in the reconstructed and the rigid lattice, respectively.Due to the reversed spin-orbit coupling in one of the layers, the energetic distance between the interlayer and the intralayer KΛ exciton is increased.This consequently makes the KΛ band structure in H-type stacking less hybridized than in R-type (cf.Fig. 8.a-b).However, the tunneling still remains the dominant component of the moiré potential as can be perceived from the colorgradient of the bands.
When considering the reconstructed lattice we also observe a decreased distance between the KΛ exciton and its respective KK exciton compared to the rigid lattice.This comes as a consequence of the valley dependence of the scalar strain potential that was discussed in Fig. 1.g in the main part of the manuscript.In Fig. 3.c we show the KΛ exciton wave function for the reconstructed lattice, which is, in contrast to the rigid lattice, more spread out and the probability distribution is pushed towards the edges of the H h h domains.This comes as a consequence of the competing potentials between electron tunneling and scalar strain, where the most efficient tunneling occurs at the H h h sites and the scalar strain is localized at the boundary of this domains.In the rigid lattice, the additional strain potential is not present, which means that KΛ exciton wave function is only determined by the tunneling strength, thus yielding a strong localization at the H h h site.

II. ATOMIC RECONSTRUCTION IN R-TYPE STACKING
In this section, we cover the impact of atomic reconstruction on the exciton energy landscape and its optical properties in R-type stacked WSe 2 .Similarly to Fig. 1  In contrast to H-type stacking, which exhibits a hexagonal pattern, here we instead have large triangular domains forming.Since the R M h and R X h stacking are both equivalent in terms of stacking energy, these regions will grow to be the same size.Meanwhile, the less energetically favorable R h h stacking will drastically shrink in size.This is in stark contrast to the rigid lattice (cf.Fig. 4.b) where it is clearly visible.

in
The resulting band edge variation for the conduction band is shown in Fig. 4.c for the bottom layer and Fig. 4.e for the top layer.Similarly, the variation for the valance band is shown in Fig. 4.d for the bottom layer and Fig. 4.f for the top layer.Here, we have two contributing potentials.The first is the alignment shift that is already present in the rigid lattice, but now with a triangular geometry due to the change in local stacking configuration.
The second is the piezo potential that stems from the shear strain induced by atomic reconstruction.Note that the scalar potential present in H-type stacking does not appear in R-type stacking for homobilayers.This comes as a consequence from the domain wall formation, which is nearly purely formed from shear strain [2], in good agreement with experimental observations [1].Comparing Fig. 4.c with Fig. 4.d we see that the conduction and valance band is shifting in the same direction in the same layer.Since there is no scalar potential present, the bands also shift the same amount, making for a vanishing moiré potential for intralayer excitons.In contrast, by comparing Fig. 4.c with Fig. 4.f we find that the conduction band is shifted down where the valance band is shifted up, resulting in a very efficient moiré potential in the case of interlayer excitons.This is already expected in the rigid lattice.Here, we have, however, a significant contribution to the depth from the piezo potential, which stands in contrast to H-type stacking, where the piezo potential is the same for both layers, making it negligible when considering the exciton potential.
InFig.4.g the effective exciton potential depth is shown for interlayer KK and interlayer KΛ excitons as a function of the twist angle.Here, we find that the KK and KΛ states have the same dependence on θ up until larger angles where the size of the mBZ increases such that the exciton wave function overlap starts to suppress the potentials, thus yielding a difference between KK and KΛ at larger angles.The potential minimum is around θ = 1.2 • and then it starts to increase again due to the increased size of the supercell -in good agreement with previous theoretical predictions [3].The spatial map of the KK interlayer potential is shown in Fig. 5.a-b.Here, we can clearly see the formation of triangular potential wells, stemming from the alignment shift and piezo potential, and consequently trapping interlayer excitons in real space.Furthermore, the sign of the potential between the different layer configurations is different, but with the same depth.Consequently, this means that the interlayer excitons are degenerate, but trapped in different regions.

A. Moiré exciton energy landscape
By solving the moiré exciton eigenvalue equation in Eq. 30 we obtain the moiré exciton energy landscape in the presence of atomic reconstruction.In Fig. 6.ab we show the moiré exciton band structure for R-type stacked WSe 2 .Similarly to the band structure for H-type shown in the main part of the paper, we find multiple strain-induced flat bands for the interlayer excitons (cf.Fig. 6.b).However, in contrast to H-type stacking, the intralayer excitons are not subject to any moiré potential and will still exhibit a very pronounced curvature.In the rigid lattice, the lowest lying KK state is an intralayer state trapped at the R h h sites due to the weak tunneling of carriers at these sites.In the reconstructed lattice, the R h h stacking configuration will drastically shrink and thus make the tunneling even less efficient, in turn making the lowest lying KK intralayer state more curved than in the rigid lattice and slightly blue shifted.Furthermore, we also find that the interlayer exciton state becomes the lowest lying in R-type stacking, similarly to H-type stacking, due to the very efficient moié potential for interlayer excitons.
In Fig. 6.c we show the exciton wave function in the reconstructed lattice for θ = 0.6 • .Here, the wave function is localized at the R M h sites with a degenerate exciton at the R X h sites due to the symmetry of the potentials (cf.trapped at R X h , implying tunability of localization via an external electrical field [4].Additionally, the wave function is following the geometry of the moiré potential (cf.Fig. 4.c).In similarity with H-type stacking, the wave function trapping site is changed in comparison with the rigid lattice.In the rigid lattice the trapping of the wave function is instead determined by the very weak tunneling around the K valley (note that this trapping is very inefficient, which can be observed from the spacing of the bands in Fig. 6.b)In R-type stacking, the only high symmetry stacking which allows for tunneling around the K valley is R h h , which consequently traps the excitons in the rigid lattice at this site.This restriction on the tunneling stems from the C 3 symmetry of the d-orbitals around the metal atom which mainly composes the orbitals of the K valley [5,6].

B. Optical response
With the moiré exciton band structure for R-type stacked WSe 2 (cf.Fig. 6.a-b) in the reconstructed and rigid lattice regime, we can now calculate the optical response of the material.In Fig. 7 we show the optical absorption for the reconstructed (blue peaks) and the rigid lattice (gray peak).In contrast to H-type stacking where we predicted a very stark difference between the rigid lattice and the reconstructed lattice, in R-type stacking the difference is negligible.This comes from the lack of a strain-induced moiré potential in R-type stacking.We observe a slight blue shift of the reconstructed peak at smaller twist angles stemming from the less efficient tunneling of carriers when the R h h regions shrink.This also removes the small splitting that otherwise can be found in the rigid lattice (see gray peaks).When increasing the twist angle the peaks align themselves, indicating that the reconstructed lattice is now more rigid.Furthermore, the KK interlayer exciton that does experience a significant moiré potential in the reconstructed lattice is not visible in the absorption spectra due to the small oscillator strength.

C. KΛ exciton band structure
In the main part of the paper we discussed the optical response in the form of PL.There, the low-lying KΛ exciton is dominating the optical response.Here, we show the corresponding moiré exciton band structure for the KΛ exciton in R-type stacked WSe 2 for both the reconstructed lattice (Fig. 8.a) and the rigid lattice (Fig. 8.b).Note that the energies are given with respect to their corresponding A exciton.We notice first the strongly hybridized nature of the KΛ exciton in both the rigid lattice and the reconstructed lattice [5], stemming from the efficient electron tunneling around the Λ valley.This consequently means that the dominating component of the moiré potential is the tunneling, which is present in both the rigid and reconstructed lattice.However, there are some significant changes between the reconstructed case and the rigid case.Notably, we find far more trapped states in the reconstructed lattice (cf.Fig. 8.a) in comparison to the rigid lattice (cf.Fig. 8.b).This is a direct consequence of the additional strain-induced potentials in the reconstructed lattice.In Fig. 8.c-d we show the Kλ exciton wave functions for the reconstructed lattice (Fig. 8.c) and the rigid lattice (Fig. 8.d).In contrast to the KK exciton wave function, the KΛ is already trapped at R M h (R X h ) in the rigid lattice due to the strong tunneling of electrons at these sites.The increased size of these regions, however, increases the size of the wave functions as well, making the KΛ exciton more delocalized in the reconstructed lattice than in the rigid lattice.

III. THEORETICAL MODEL A. Atomic reconstruction
To obtain the moiré exciton landscape in the regime of atomic reconstruction we first need to access the change in relative atomic positions with the twist angle, i.e. the atom deviation from the rigid lattice site.In a continuum model, this deviation is captured by a displacement vector u(r), which translates a point in space r with some distance u(r).Assuming only small displacements, the displacement vector is then related to the linear strain tensor as ε ij = 1 2 (u i,j +u j,i ), where u i,j = 1 2 (∂ i u j +∂ j u i ), where i(j) = (x, y).Following the theory of elasticity, the total elastic energy per unit volume can then be written as [7] where λ and µ are the material specific Lamé parameters and l is the layer index.We now consider the presence of another layer with a marginal twist angle.Following the approach laid out in Refs.[2,3] we use a parameterized form of the adhesion energy and fit it to data from density functional theory (DFT) where R/H denotes R-type and H-type stacking respectively, and the phase is given by γ R = π/2(γ H = 0).
Here, r 0 = θẑ × r + u t (r) − u b (r), which means that G n r 0 ≈ g n r + G n ∆u(r), where ∆u(r) = u t (r) − u b (r) and g n is the reciprocal vector of the mini-Brillouin zone (mBZ).Furthermore, Z R/H (r 0 ) is the deviation of the interlayer distance from d 0 = 0.69 nm [3], given by The parameters κ, a 1 , a 2 and A are all fitted from DFT simulations and are obtained from [3].The parameters are also given in Table I.The total energy can then by described with the following functional integral [3,8] where A M is the moiré unit cell area and r is the real space coordinate in the moiré lattice.To find the relevant displacement vectors u l (r), we expand them as a Fourier series u l (r) = n u l n e ignr and turn the problem into an optimization problem, which can then be solved numerically for the Fourier coefficients [3,8].This is done for the WSe 2 homobilayer up to the 12th moiré shell of the mBZ lattice vectors.

B. Strain in the reconstructed regime
The displacement vectors obtained from the optimization problem described in the previous section can be TABLE I. Fitted parameters for the adhesion energy and the Lamé parameters of bilayer WSe2 obtained from Ref. [3].related for small displacements to the strain in the material via the linear strain tensor In our considered material, this is a second-rank tensor belonging to the D3h symmetry group.Consequently, in 2D u l i,j transforms according to scalar component u l x,x + u l y,y and a vector component (u l x,x − u l y,y , −2u l x,y ) [9][10][11].Since we are interested in the effect reconstruction has on the exciton energy landscape, we are mainly interested in what impact the strain has on electronic bands.For this purpose, we can associate the scalar component as uniaxial strain in each direction and the vector component as a vector gauge potential [12], also known as piezo potential [2,13].
Scalar strain potential: The band shifts due to the scalar component of the strain are obtained by exploiting the linear dependence of the band on the uniaxial strain [14] where ξ is the valley index, λ = (c, v) the band index, and g l ξλ the valley-specific gauge factor obtained from the DFT calculations preformed in Ref. [14] and given in Table II.
Piezo potential: The relationship between the polarization P l i and the linear strain tensor is given by [15] where e l ijk is the piezoelectric tensor component.Similarly to the strain tensor, the piezoelectric tensor belongs to D3h and will thus only have certain components contributing.These are [16] where we have dropped the third index and only use e l 11 as the piezo coefficient from now on.The resulting polarization is the vector component of the linear strain tensor established in earlier multiplied by the piezo coefficient By aligning the y-axis of the polarization with the vector from the metal atom to the chalcogen atom we have e 0 11 = −e 1  11 > 0 for H-type configuration and e 0 11 = e 1 11 for Rtype configuration [2].The value of the piezo coefficient e 11 = 2.03 • 10 −10 C/m is obtained from Ref. [17].Note that there is no polarization in uniformly strained lattices without a shear component.Via Gauss law one can write the charge distribution with respect to the polarization By including the out-of-plane direction it reads In addition to the piezo charge density, there is also the screening-induced charge density given by the Poisson equation where ϕ(r, z) is the electrostatic potential produced by the piezo charges and α l ∥ = d 0 ϵ 0 (1−ϵ l ∥ ) is the in-plane polarizability [2].Here, ϵ l ∥ is the in-plane dielectric screening obtained from Ref. [18].The total charge density ρ tot = ρ ind + ρ piezo and the piezoelectric potential can then be found by solving the full Poisson equation, which is done by expanding the piezo potential as a Fourier expansion ϕ(r, z) = n φn (z)e ign•r and solving for matching boundary conditions of the two dielectric slabs The piezo-induced energy shifts of a charge carrier P λ l (r) is then directly obtained from the piezo potential for each layer ϕ l (r).

C. Derivation of the moiré exciton Hamiltonian
We obtain microscopic access to the moiré exciton energy landscape by first considering a Hamiltonian in second quantization, starting in a decoupled monolayer basis [4][5][6]19].Here, we take into account the strong Coulomb interaction by solving the generalized Wannier equation [20] and add the moiré potential as periodic modifications to these energies.In the regime of atomic reconstruction we will have the two strain induced potentials including the scalar strain potential S(r) (see Eq. 7) and the vector strain potential, which we will denote as piezo potential P (r).
In addition to the strain induced potentials, we will also have the components that already exist in the rigid lattice.These are the stacking-dependent alignment shift V (r) [6,19] and the tunneling of charge carriers T (r) [4-6, 21, 22].The stacking dependent alignment shift V (r) stems from a charge transfer-induced polarization between the layers [23], which consequently give rise to an electrostatic potential that will be dependent on the relative atomic alignment between the layers.Since this spontaneous charge polarization can only happen if we lack inversion symmetry, it will only be present in Rtype stacked configurations [23].The interlayer tunneling T (r) stems from the wave function overlap between the layers giving rise to a mixing of monolayer eigenstates, which can be interpreted as a tunneling of electrons and holes between the layers [5,6,21,24,25].Since the wave function overlap is largely dependent on the interlayer distance, the tunneling rate is also periodic when introducing a twist angle, due to the variation in interlayer distance throughout the supercell [6,21,22,26].Thus, in the reconstructed regime we have four different components of the moiré potential.The corresponding moiré Hamiltonian in real space then reads where h (′) = (l, ξ) is a compound index, taking into account the layer index l and the valley index ξ.Furthermore, Ψ ( †) (r) are the annihilation (creation) operators in real space.In the effective mass approximation, we can approximate the wave functions close to the high symmetry points in the Brillouin zone (BZ) as a plane wave expansion Ψ λ † h (r) = k e ik•r λ † h,k , where k is the momentum of the electron/hole.In momentum space, the Hamiltonian will thus have the following form where the matrix elements are given by the Fourier coefficients in the expansion of their real space constituents Here, p λ hh (g) is obtained by solving the integral for the Fourier coefficients The area of the moiré supercell is denoted as A M .Note that we have used the piezo potential here, but the same method to extract the coefficients holds for all four contributions to the moiré potential.
In the case of the tunneling and the alignment shift, we can obtain the real space map of the potentials by smoothly interpolating between different high symmetry stackings [4,6,19] where v λ h , A λ h and B λ h are parameters which are fitted to data from density functional theory Ref. [6].From section III A we have G n r 0 ≈ g n r + G n ∆u(r), where ∆u(r) = 0 for the smooth interpolation of the rigid lattice.We then take into account the deformation of the lattice, by allowing the displacement vectors to deform the geometry, in turn yielding After obtaining the full moiré Hamiltonian in electronhole basis, we can now transform it into exciton basis by considering the pair operator expansion.Under the assumption of low densities, we can expand the single particle operators in pairs D † hk,h ′ k ′ = c † hk v h ′ k ′ .Consequently, we obtain the following expression [4,5,19,27] We can now transform into exciton basis by expanding the pair operators in a set of orthogonal basis functions where α hh ′ (β hh ′ ) = m c(v) h(h ′ ) /(m c h + m v h ′ ), ξ h is the valley index of the compound h and µ is the exciton quantum number, which in this work is restricted to be 1s.Furthermore, Ψ µ hh ′ (k) is the exciton wave function that solves the Wannier equation [20].Moreover, X ( †) are the exciton annihilation (creation) operators.Performing the above transformation now allows us to write the complete moiré exciton Hamiltonian Here, L = (l e , l h ) is a compound layer index, Q is the center-of-mass momentum and ξ = (ξ e , ξ h ) is the compound exciton valley index.Additionally, we have introduced the notation M ξ Lg for the total layer-diagonal moiré potential given by M ξ Lg = S ξ Lg + P ξ Lg + V ξ Lg , which includes the scalar strain potential S ξ L (g), the piezo potential P ξ L (g) and the alignment shift V ξ L (g).Furthermore, E ξ LQ is the exciton dispersion given by where E b ξ are the exciton binding energies as obtained from the generalized Wannier equation [20] and ε λ ξ λ is the valley splitting.Effective electron (hole) masses and valley splittings are obtained from Ref. [28].The matrix elements of the Hamiltonian Eq. 23 then read where m λ l λ (g) = s λ l λ + p λ l λ + v λ l λ is the short notation for the layer-diagonal Fourier coefficients, while t λξ λ l λ l ′ λ describes the tunneling Fourier coefficients, both obtained from Eq. 18.Here, F ξ LL ′ (q) are the form factors given by F ξ LL ′ (q) = k Ψ ξ * L (k)Ψ ξ L ′ (k + q).

D. Diagonalization of moiré exciton Hamiltonian
Once the complete moiré exciton Hamiltonian in the reconstructed regime (see Eq. 23) is developed, we now find a diagonal form of this Hamiltonian.We do this by applying the well-known zone-folding scheme [4,5,19].Here, we restrict the summation over the center-of-mass momentum Q to the first mBZ and then fold it back in with the mBZ lattice vectors g where Q ∈mBZ and the notation E ξ LQ (g) = E ξ LQ+g .By changing to the zone-folding operator basis F ξ LQg = X ξ L,Q+g we get the following expression Here, we have introduced the abbreviation M ξ L (g, g ′ ) = M ξ L (g ′ − g).
We now introduce the hybrid moiré exciton basis, where we expand the zone-folded operators with a set of orthogonal basis functions called mixing coefficients [4][5][6]19] where η is the new moiré exciton band index and C ξη Lg (Q) are the mixing coefficients that fulfill the following conditions Applying the above transformation to Eq. 27 we reveal the moiré exciton eigenvalue equation where E ξ ηQ are the final moiré exciton energies, which can be obtained by numerically solving the equation above.Consequently, we obtain the diagonal form of the moiré exciton Hamiltonian

E. Optical response
Our approach for the optical response of the material is done in the exact way as described in the supplementary of Ref. [4].Here, we use the photoluminescence (PL) formula for phonon-assisted PL first derived in Refs.[5,29], thus allowing us to obtain the optical response of both bright and momentum dark excitons.

FIG. 1 .
FIG. 1. Spatial map of the interlayer distance in (a) a reconstructed lattice and (b) a rigid lattice at θ = 0.6• for H-type stacked WSe2 homobilayers.Here, the interlayer distance is associated with a high-symmetry stacking.The band edge variation for the conduction band in (c) the bottom layer and (e) the top layer at the K-point, and the valance band edge variation in (d) the bottom layer and (f ) the top layer at the K-point.Note the reversed sign between the layers and the difference in the potential depth between valance and conduction band.Since this is a purely strain-induced potential, it is not present in the rigid lattice.The shape of the potential stems from the linear combination of the scalar strain potential and the piezo potential.(g) Maximum exciton potential depth ∆Emax in the supercell is shown as a function of twist angle for both intralayer excitons (dashed) and interlayer excitons (solid).

FIG. 2 .
FIG. 2. Band structure and wave function of the bright KK exciton in H-type stacked WSe2 bilayer at θ = 0.6 • in (a,c) the reconstructed lattice and (b,d) the rigid lattice.Energies are normalized to the lowest-lying KK exciton in the rigid lattice.Note that the reconstruction drastically changes the localization site of excitons and that the efficient moiré potential for interlayer excitons shifts them below the intralayer KK states.The band index η = 0 indicates which wave function is shown.

FIG. 3 .
FIG. 3.Normalized photoluminescence spectrum at 4K at different twist angles, displaying phonon replicas of the energetically lowest KΛ exciton in (a) H-type and (b) Rtype WSe2 homobilyers.The spectra from the reconstructed and the rigid lattice are shown in blue and gray, respectively.A signature of lattice reconstruction is a more pronounced red-shift of the exciton resonance for smaller twist angles.

FIG. 4 .
FIG. 4. Twist-angle dependence of the normalized absorption in H-type WSe2 bilayers for (a) the rigid lattice and (b) the reconstructed lattice.The emergence of multiple peaks in the latter is due to the strain-induced intralayer potential and is an unambiguous signature of lattice reconstruction.

FIG. 3 .
FIG. 3. Band structure for the low-lying dark KΛ exciton in H-type stacked WSe2 at θ = 0.6 • in (a) the reconstructed lattice and (b) the rigid lattice.Energies are normalized to the lowest lying A exciton in their respective lattice.The corresponding KΛ exciton wave function for the lowest lying state is shown for (c) the reconstruction lattice and (d) the rigid lattice.

FIG. 4 .
FIG. 4. Spatial map of the interlayer distance in (a) a reconstructed lattice and (b) a rigid lattice at θ = 0.6 • for R-type stacked WSe2 homobilayers.Here, the interlayer distance is associated with a high-symmetry stacking.The band edge variation for the conduction band in (c) the bottom layer and (e) the top layer at the K-point, and the valance band edge variation in (d) the bottom layer and (f ) the top layer at the K-point.Note the reversed sign between the layers.The shape of the potential stems from the linear combination of the the piezo potential and the alignment shift.(g) Maximum exciton potential depth ∆Emax in the supercell is shown as a function of twist angle for interlayer excitons.

FIG. 5 .
FIG. 5. Effective moiré exciton potential depth at θ = 0.6 • for the KK interlayer exciton for the electron in (a) the top layer and (b) for the electron in the bottom layer.Note that the sign is reversed when changing the layer configuration.
FIG. 6. Band structure and wave function of the bright KK exciton in R-type stacked WSe2 bilayer at θ = 0.6 • in (a,c) the reconstructed lattice and (b,d) the rigid lattice.Energies are normalized to the lowest-lying KK exciton in the rigid lattice.Note that the reconstruction drastically changes the localization site of excitons and that the efficient moiré potential for interlayer excitons shifts them below the intralayer KK.

FIG. 7 .
FIG. 7. Absorption spectrum showing the the bright KK exciton for different twist angles in the R-type WSe2 homobilayer.Here, the blue peaks stem from the reconstructed lattice and the gray peaks from the rigid lattice.

FIG. 8 .
FIG. 8. Band structure of the low-lying dark KΛ exciton in R-type stacked WSe2 at θ = 0.6 • in (a) the reconstructed lattice and (b) the rigid lattice.Energies are normalized to the lowest lying A exciton in their respective lattice.The corresponding KΛ exciton wave function for the lowest lying state is shown for (c) the reconstructed and (d) the rigid lattice.

TABLE II .
Linear strain gauge factors for WSe2 for singleparticle bands, obtained from DFT calculations Ref.[14].