Excitons in mesoscopically reconstructed moiré heterostructures

Moiré effects in vertical stacks of two-dimensional crystals give rise to new quantum materials with rich transport and optical phenomena that originate from modulations of atomic registries within moiré supercells. Due to finite elasticity, however, the superlattices can transform from moiré-type to periodically reconstructed patterns. Here we expand the notion of such nanoscale lattice reconstruction to the mesoscopic scale of laterally extended samples and demonstrate rich consequences in optical studies of excitons in MoSe2–WSe2 heterostructures with parallel and antiparallel alignments. Our results provide a unified perspective on moiré excitons in near-commensurate semiconductor heterostructures with small twist angles by identifying domains with exciton properties of distinct effective dimensionality, and establish mesoscopic reconstruction as a compelling feature of real samples and devices with inherent finite size effects and disorder. Generalized to stacks of other two-dimensional materials, this notion of mesoscale domain formation with emergent topological defects and percolation networks will instructively expand the understanding of fundamental electronic, optical and magnetic properties of van der Waals heterostructures.


Supplementary Note 1: Mesoscopic reconstruction in SEM
To image mesoscopic lattice reconstruction in MoSe 2 -WSe 2 HBLs (as in Fig. 1f,g of the main text) we used the secondary electron (SE) imaging technique in low-energy scanning electron microscopy (SEM) initially pioneered on SiC crystals [1] and recently adopted to twisted WSe 2 homobilayers [2]. In brief, the technique is based on out-lens detection of SEs emitted from atoms near the crystal surface by inelastic scattering of incoming primary electrons. In a TMD bilayer, the top and bottom hexagonal lattices form an effective cavity for electrons, rendering the scattering probability for the incoming beam in the cavity and thus the generation of SE dependent on the relative atomic positions in given stacking configurations. In optimized experimental geometries, the contrast of SE yield under channeling conditions allows to discriminate domains of different stacking configurations [2]. We adopted the technique to image R-and H-type MoSe 2 -WSe 2 HBLs. The samples were prepared by the same method as for optical spectroscopy yet without hBN encapsulation.
For both R-and H-type HBLs, the measurements were performed in a Raith eLine system The images illustrate for nearly-aligned H-and R-type MoSe 2 -WSe 2 HBLs clear departures from canonical moiré patterns or uniformly reconstructed periodic domains. They actually exhibit large variations in domain shape and size with the following common features: from the tips and edges to the cores, the reconstructed pattern evolves from large micron-sized two-dimensional (2D) domains to elongated one-dimensional (1D) stripes which merge into a network of zero-dimensional (0D) arrays with dimensions well below 100 nm with variations in domain shape and size. In addition to line defects such as tips and edges, local strain such as interfacial bubbles in the bottom right panel of Fig. 1 can also give rise to mesoscopic reconstruction. The main difference between H-and R-type HBLs is that only one staking (H h h ) is optimal in the former, whereas two stackings (R X h and R M h ) yield optimally reconstructed domains in the latter.

Supplementary Note 2: Modelling of mesoscopic reconstruction
Reconstruction of MoSe 2 -WSe 2 HBL is driven by the interplay of interlayer adhesion energy and strain. The sum of the intralayer and interlayer energies of the lattice is given by the integral over the HBL area S as [3]: where U (r) and W s (r) are the respective intralayer and interlayer strain energy densities which depend on the in-plane displacements in MoSe 2 and WSe 2 layers. Assuming equal elastic properties and lattice constants in both layers [4], the displacement fields are given by: The intralayer strain energy density is given by: with the strain tensor: where the first and second Lamé parameters λ and µ are obtained by averaging over the respective parameters in MoSe 2 and WSe 2 MLs.
The interlayer term W s (r) = V s (r) − εZ 2 s (r) quantifies the adhesion energy density defined via the expressions [5]: {A 1 e −Qd 0 cos φ n (r) + A 2 e −Gd 0 sin [φ n (r) + ϕ s ]}, Z s (r) = 1 2ε In the next step, we assume an equilateral triangle for the shape of the HBL area S.
Supplementary Fig. 3a shows the triangle baseline placed on the x-axis as a borderline to the moiré region in the HBL with twist angle θ and without reconstruction. At the borderline, the displacement field is zero, u x (x, 0) = u y (x, 0) = 0, whereas the remaining two sides of the triangle are free from boundary conditions. We assume that the y-axis bisects the angle θ and divides the equilateral triangle into two right-angled triangles, yielding two constrains: u y (x, y) = −u y (−x, y).
Finally, we discretize the displacement field u(r) with a square mesh, replace the integration in Eq. (1) by the summation over the triangle area, and express the derivatives in Eq. (4) through the finite differences method. The density of the square mesh (we choose 288 moiré unit cells in all simulations of equilateral triangles) determines the number of discretization parameters (i.e. if the mesh divides the height of triangle h into N equal sectors, the displacement field is defined by ∼ (N + 1) 2 / √ 3 parameters). To find the parameter set that minimizes the total energy in Eq. (1), we use the trust-region algorithm implemented in the Optimization Toolbox TM of MATLAB ® to determine the local energy minimum of the final displacement field u f as a function of the initial displacement field u i set prior to reconstruction. In all calculations, we monitor the convergence of numerical results.
In Fig. 2a,d of the main text we used the following initial sets of displacement fields. The ideal moiré case was obtained with the initial set u i = 0 (twisted HBL with θ = 0.4 • ) and without reconstruction. The periodic case illustrates the final set after the reconstruction with the initial field u i = 0. All remaining cases show the final sets with the rotated initial displacement: which realizes the rotation of layers by ±θ/2 around the point (0, αh) and leads to untwisted HBL regions when the stacking expands from the rotation point to the entire area of the triangle. To denote these initial displacement fields we use the dimensionless y-coordinate of the rotation point, α. The top panel shows periodic reconstruction, and the central panel illustrates reconstruction with one rotated point at α = 0. The maps in the bottom panels were calculated for two rotation points at α = 0 and α = 1/12 with the initial displacement field: This initial displacement field realizes three areas upon reconstruction: two untwisted domains around points with α = 0 and α = 1/12 with zero initial displacement in between. The resulting patterns exhibit regions of both periodically reconstructed and rotated patterns.
This procedure can be expanded to model realistic samples with complex reconstruction patterns as in Supplementary Fig. 1 and 2 by taking into account multiple rotation points for a given sample geometry and strain distribution.      The PL spectra were recorded with an excitation power of 2 µW; for 0D arrays, additional spectra at 0.01 µW excitation power are shown (light and dark purple for R-type, light and dark brown for H-type) and scaled for better visibility (factors of 10 and 100 for R-type, 100 and 1000 for H-type).      Supplementary Figure 11. Characteristics of sample 4 without DR data.  Figure 12. Power-dependent interlayer PL of different domain types in nearlyaligned R-type HBL. a -c, Normalized PL spectra at selected excitation powers for a 2D domain (a), a region with 1D stripes (b), and domains of 0D arrays (c). All data were acquired on sample 1 with excitation at 725 nm.  Zeeman splitting between σ + and σ − polarized peaks in a magnetic field of 9 T under linearly polarized excitation (π). Bottom panels: Valley Zeeman splitting ∆E as a function of magnetic field. The solid lines are linear fits to the data with g-factors and error bars obtained from leastsquare best fits. The negative P c and the positive g-factor assign the PL to singlet R X h interlayer exciton. All data were recorded with excitation powers of 20 and 0.01 µW for 2D and 0D regions.     Supplementary Fig. 20, where the amplitude of the primary component dominates at high excitation powers over the amplitudes of the two complementary channels with power-independent ratio.
With this assignment of primary and secondary PL decay components to non-radiative population loss processes (including Auger) and radiative decay, respectively, we attribute the third channel of H h h triplet exciton PL decay on 440 ns timescale and R X h singlet exciton state with 45 ns decay time to population feeding from energetically higher-lying long-lived reservoirs, presumably constituted by momentum-dark configurations of valence band holes in WSe 2 and conduction band electrons in MoSe 2 . Note that the absence of such reservoirs could explain the lack of the third channel in the PL decay of the H h h singlet state that is the top-most energy level according to our calculations summarized in Table 1

Domain
Stacking Spin configuration  Supplementary Fig. 19 using tri-exponential decay with characteristic lifetime (τ i ) and amplitude (A i ) of each decay channel. For large 2D domains, exemplified in Supplementary Fig. 22 a and b, the interlayer PL is bright and characterized by a single peak at ∼1.40 eV of the triplet exciton state in H h h stacking. Consistently, the corresponding DR spectra feature single-peak resonances of intralayer excitons in MoSe 2 and WSe 2 . In these domains, the degree of circular polarization is high and the degree of linear polarization is zero, as indicated by numbers in Supplementary Fig. 23 a and b. For regions with elongated 1D domains, exemplified in Supplementary Fig. 22 c and d, the interlayer exciton PL is detected around 1.40 eV or with slightly red-shifted additional peaks. In contrast to 1D stripes in R-type HBLs, the PL preserves a sizable degree of circular polarization, whereas the degree of linear polarization is either negligible or very small where the density of 1D stripes is high (Supplementary Fig. 23 c and d) For 0D array regions, finally, exemplified in Supplementary Fig. 22 e and f, the interlayer exciton PL is very low and red-shifted by several tens of meV compared to 2D domains. This red-shift is size-dependent, with smaller domains exhibiting larger red-shifts. Moreover, for excitation powers below 100 nW, the PL spectra develop into quantum dot like spectrally narrow peaks. As shown in Supplementary Fig. 23 e and f,    Supplementary Figure 23. af, PL spectra recorded in circular (left panels, with σ + and σ − polarized detection shown in red and blue, respectively) and linear (right panels, with two orthogonal orientations of linear polarization for maximum and minimum PL intensity shown in purple and green, respectively) basis from the same regions as in Supplementary Fig. 22 af. The numbers in each panel denote the respective degrees of circular and linear polarization, P c and P l . In the same manner, we performed secondary electron modulated SEM and PL hyperspectral imaging on the R-type sample and co-aligned the two images ( Supplementary Fig. 24).
From the direct correlation between SEM and PL images, we again obtained the results consistent with the analysis in the main text on reconstruction patterns and their optical characteristics of R-type samples. For 2D domains of R X h stacking, exemplified in Supplementary Fig. 25 a and b, their interlayer exciton PL is bright and the spectra feature a single peak at ∼1.33 eV from the singlet exciton state. Consistently, their DR spectra show single-peak resonances of intralayer excitons in MoSe 2 and WSe 2 . The interlayer exciton PL from these 2D domains has negative degree of circular polarization and zero degree of linear polarization ( Supplementary Fig. 26 a and b). For 1D stripe regions, exemplified in Supplementary Fig. 25 c and d, the interlayer exciton PL is blue-shifted and broadened compared to that from 2D domains of R X h stacking. This PL features almost zero degree of circular polarization but high degree of linear polarization along the elongated direction of the stripes. Finally, for 0D array regions, exemplified in Supplementary Fig. 25 e and f, the interlayer exciton PL is reduced and blue-shifted by several tens of meV compared to 2D domains. In contrast to H-type, the smaller 0D domains show larger blue-shifts because of the quantum confinement. When excitation powers are below 100 nW, the PL spectra of 0D domains have the characteristic of quantum dots, showing narrow peaks with negative degree of circular polarization and zero degree of linear polarization (Supplementary Fig. 26 e and f ). The DR spectra of 0D regions in R-type exhibit for both MoSe 2 and WSe 2 intralayer exciton resonances characteristic splittings and broadenings with smaller domains featuring larger energy splittings, which is similar to H-type.  Supplementary Figure 26. af, PL spectra recorded in circular (left panels, with σ + and σ − polarized detection shown in red and blue, respectively) and linear (right panels, with two orthogonal orientations of linear polarization for maximum and minimum PL intensity shown in purple and green, respectively) basis from the same regions as in Supplementary Fig. 25 af. The numbers in each panel denote the respective degrees of circular and linear polarization, P c and P l .

Supplementary Note 6: Analysis of intralayer exciton absorption spectra in reconstructed domains
In the following, we analyze the DR spectra of intralayer excitons in Fig. 3a and 4a (main text) recorded on different positions of reconstructed R-and H-stacks. First, we note that the DR spectra exhibit asymmetric lineshapes due to multiple interferences at the interfaces of samples containing stacks of hBN, TMD, SiO 2 and Si. To simplify the lineshape analysis, we transform the DR spectra into absorption spectra by adapting the numerical method developed in Ref. [6].
We begin with a simple case of two interfaces provided by a thin TMD layer on top of SiO 2 in vacuum. Since the thickness of TMD layer d is much smaller than the wavelength of light λ, DR (defined as ∆R/R) can be written as [7]: where ǫ 1,2,3 is the dielectric function of vacuum, TMD layer and SiO 2 . In linear response approximation, ǫ 1,2,3 = 1 + χ 1,2,3 with dielectric susceptibility of each layer χ 1,2,3 . With vanishing absorption in vacuum and SiO 2 at the wavelengths of relevance in the spectral window of intralayer exciton transitions in MoSe 2 and WSe 2 , Im(ǫ 1 ) = Im(ǫ 3 ) = 0, and thus DR in Eq. 12 is simply proportional to the imaginary part of the optical susceptibility χ ′′ 2 of the thin TMD layer, which in turn corresponds to its absorption.
To model the DR response of our samples with two hBN layers that embed the TMD layer on the top and bottom sides and ∼ 500 µm Si substrate underneath the SiO 2 layer, we account for additional interfaces between the top hBN layer and the TMD layer, the TMD layer and the bottom hBN layer, its interface to SiO 2 , and the SiO 2 -Si interface by a phase factor e iα [8] in the effective susceptibilityχ 2 = e −iα 0 χ 2 , where we approximate α by a constant α 0 in the relevant wavelength range. By decomposingχ 2 into real and imaginary parts asχ ′ 2 + iχ ′′ 2 , the imaginary part of the effective susceptibilityχ ′′ 2 can thus be expressed as χ ′′ 2 = cos(α 0 )χ ′′ 2 + sin(α 0 )χ ′ 2 . Finally, as DR is proportional toχ ′′ 2 , using Kramers-Kronig relations we obtain: where ω denotes the angular frequency. Using Eq. 13, we compute χ ′′ 2 (ω) of the exciton transitions in the regions of MoSe 2 and WSe 2 monolayers from the respective DR spectra in the charge-neutral regime for different values of the phase shift α 0 from −π to π. Since the exciton absorption of monolayer TMD is expected to have a Lorentzian lineshape in the absence of free charge carriers [9], we choose the value of α 0 such that the obtained absorption spectra χ ′′ 2 (ω) can be accurately fitted by a Lorentzian. Successively, we use this value of α 0 to compute the absorption spectra of intralayer excitons in HBL regions as the transition energies are similar to monolayer exciton transitions. The respectively obtained absorption spectra χ ′′ 2 (ω) in Supplementary Fig. 27 and 28 are predominantly positive, providing confidence in our approach.
In the following analysis of absorption spectra, we assume intralayer exciton transition energies are sensitive to the local stacking and domain size. In regions of two proximal stackings A and B (R X h and R M h , and H h h and H X h in the case of 0D reconstructed domains of R-and H-stacks, respectively) with the corresponding energy minima, coupling of the intralayer exciton states via tunneling gives rise to a doublet peak structure. By introducing tunnel-hopping between the domains of different registries, we model the eigenenergies of the coupled two-level system with the following Hamiltonian: where E A and E B are the intralayer exciton energies for the respective commensurate stackings, t and β are the hopping parameters, and a m is the domain size (or moiré period).
Using t = 20 meV and β = 0.07 nm −1 estimated from the calculations of Yu et al. [10], we compute the energies (in meV) in Supplementary Tab. 2 by fitting the peak evolution in Supplementary Fig. 27

Supplementary Note 7: Density functional theory calculations
We used density functional theory (DFT) to calculate exciton g-factors and oscillator strengths of interlayer excitons in R-and H-type MoSe 2 -WSe 2 HBLs as described in detail in Refs. [11,12]. DFT calculations were performed with the PBEsol exchange-correlation functional [13] as implemented in the Vienna ab-initio simulation package (VASP) [14].
Van der Waals interactions were included with the DFT-D3 method [15] and Becke-Johnson damping [16]. Moreover, spin-orbit interactions were included at all stages. Elementary cells with thickness of 35 Å in the z-direction were used to minimize the interactions between periodic images. The atomic positions were relaxed with a cutoff energy of 400 eV until the total energy change was less than 10 −6 eV. Calculations were performed for high-symmetry points of HBL moiré patterns in R-and H-type stackings on the Γ-centered k grid of 6 × 6 divisions with 600 bands and the cutoff energy of 300 eV.
The theoretical g-factors and oscillator strengths of interlayer excitons obtained from our DFT calculations for R-and H-type MoSe 2 -WSe 2 HBLs are listed in Table. 1 of the main text. For completeness, we list in Supplementary Table 3 interlayer exciton g-factors in different stackings and momentum-dark spin-valley configurations. Supplementary Fig. 29 gives a graphical representation of the energetic ordering of singlet and triplet interlayer exciton states with respective dipolar selection rules and oscillator strengths, labelled by the according spin configuration and atomic registry of high symmetry stackings in R-and Htype HBLs. The insets of Supplementary Fig. 29 a and  and valence band holes at K in WSe 2 or at Γ in the hybrid band of MoSe 2 -WSe 2 according to the electronic band structure [17]. For each state, the absolute g-factor values of both singlet and triplet spin configurations are listed for distinct atomic registries. Note that due to different symmetry in R-and H-type MoSe 2 -WSe 2 HBLs, the spin-valley configurations are distinct [18]: In R-type (H-type), KK (K ′ K) interlayer exciton states are momentum-direct with lowest energy states in singlet (triplet) configuration, whereas K ′ K (KK) states are momentum-indirect with triplet (singlet) lowest-energy states.