T wisted van der Waals structures, such as twisted bilayer graphene [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] (TBG), twisted double bilayer graphene [17][18][19] (TDBG), and twisted transition-metal-dichalcogenides [20][21][22][23][24][25][26][27] are in the vanguard of quantum materials research [28][29][30] . The twist between the layers leads to large-scale periodic perturbations of stacking configurations, called a moiré superlattice. Because atomic layers in van der Waals (vdW) materials are not rigid but instead behave as deformable membranes, moiré suprelattices acquire additional attributes. As two atomic layers with a small relative twist angle come in contact, the atomic positions relax to minimize the total energy. Through the relaxation process domains of lowest-energy configurations form and become separated by domain walls of transitionary configurations [31][32][33] (Fig. 1a). The generalized stacking fault energy function (GSFE), which provides the energetic variations across different stacking configurations, is the fundamental property that describes relaxed vdW interfaces 31,34 . The GSFE is commonly calculated using density functional theory (DFT) 31,34 . Experimental techniques 35 to probe the GSFE are currently restricted to the stable lowestenergy configuration, and are very limited in energy resolution compared to the variability among theoretical descriptions.

Here we show that the generalized stacking fault energy function (GSFE) is encoded in fine details of the relaxed moiré super-lattice patterns at the low twist-angle limit. In particular, the shape of domains and domain walls networks, as well as domain wall width, abide by transitionary configurations beyond the lowest-energy stackings of the domains. More specifically, we distinguish between single and double domain walls (SDW and DDW). SDWs separate two distinct stacking configurations of a moiré superlattice (for instance, ABCA [MM'] and ABAB [MX'] in the TDBG [for twisted H-stacked MoSe 2 /WSe 2 , or T-H-MoSe 2 /WSe 2 for short] example of Fig. 1a). DDWs, formed from the collapse of two SDWs, separate identical phases (ABAB for TDBG and MX′ for T-H-MoSe 2 /WSe 2 in Fig. 1a). The formation and nature of DDWs result from attraction of SDWs as they are brought together (for instance, due to external or relaxation induced strain), and is proven here to provide a reliable read-out of the underlying energetics. In cases of inequivalent two lowestenergy configurations (as in Fig. 1), the SDW develops a finite curvature κ, allowing one to extract the domains energy imbalance with an accuracy outperforming the ~3 meV/atom of the gold standard of quantum chemistry 36,37 .

Moiré metrology, presented here, correlates measurable spatial patterns of the relaxed moiré superlattice (such as shapes of domains, SDWs and formation of DDWs) with modeling based on the GSFE. To do so, we developed a continuous twodimensional relaxation simulation. The model searches for local interlayer displacement fields that minimize the total energy of the multilayer, as a sum of elastic and stacking energy terms (see Supplementary Information S1-2 for more details, also see ref. 33 for an alternative approach). The equations are solved in real space and thus capture subtle experimental details that remained underexplored. Fig. 1b-g is a tour-de-force of moiré metrology combining experimental imaging of different systems, techniques and length-scales (Fig. 1b-d), and their respective modeling (Fig. 1e-g). Fig. 1b-d were acquired with modern scanning probe microscopy (SPM) techniques: scanning tunneling microscopy (constant current mode) and spectroscopy 17 (STM and STS respectively) and mid-infrared range (mid-IR) scanning nearfield optical microscopy 11 (SNOM). These techniques resolve stacking configurations based on local topographic and electronic (STM and STS), as well as mid-IR optical conductivity (mid-IR SNOM) contrasts. In low strain TDBG, the model (Fig. 1e) captures the fine curving of SDWs (Fig. 1b). In cases of higher strain (Fig. 1c and modeling in Fig. 1f) we observe the formation of one dimensional DDW structures (inset of Fig. 1f highlights an example). Similarly, DDW formations and SDW curving were observed (Fig. 1d) and modeled (Fig. 1g) in T-H-MoSe 2 /WSe 2 , with excellent agreement across different length-scales of the image (see Supplementary Information S3 for additional analysis). Next, we will illustrate in detail how moiré super-lattices reveal the energy landscape information using TBG and TDBG as prototypical examples.

Moiré metrology of twisted bilayer graphene. To study the energy landscape of TBG, we focus on the interplay between SDW and DDW formations. Fig. 2a presents a non-local nanophotocurrent map of TBG in the minimal twist limit <0.1°. Bright spots in the photo-current map highlight the AA sites (indicating higher absorption-see "Methods" and Supplementary Information S4). The AA sites are connected by domain walls separating AB and BA domains. The resultant moiré super-lattice is clearly affected by strain, inferred from the distorted triangular pattern, especially near the edges of the stack. There, we observe the merging of two SDWs into a single DDW (selected locations are marked in Fig. 2a). We successfully account for the observed network within a model addressing a competition between SDWs and DDWs. To grasp the essential physics, we first assume a characteristic energy of forming a segment of DDW and SDW. We define a dimensionless domain-wall formation ratio as the ratio of DDW and SDW line energies, β ¼ γ DDW =γ SDW . In addition to β, the model input includes the AA sites of the moiré pattern as the fixed vertices of the triangles forming the network. We explore the SDW vs. DDW structures that emerge for a given value of the single tuning parameter β. The case of β ¼ 2 implies there is no benefit in forming a DDW, and the optimal structure would simply be straight SDWs connecting the AA sites. For β < 2 the two SDWs attract each other favoring the emergence of DDW segments (see Supplementary Information S4 and Supplementary Video 1 for details). Our modeling captures the overall shape of the experimental map for β ¼ 1:90 (Fig. 2a). The agreement is remarkable considering the minimal modeling we employ. We conclude that in order for a TBG model to reproduce the experimental picture, two SDWs have to sufficiently attract one another as quantified by the fitted β. In that sense, as we show more rigorously below, moiré metrology puts constraints the GSFE.

To quantify how the observed moiré networks constrain the stacking energy landscape, we span all realistic GSFE's satisfying the symmetry of TBG over a 2D unit-less parameter space (ζ, τ) (as illustrated in Fig. 2b and discussed at Supplementary Information S4), such that each point on the (ζ, τ) plane represents one GSFE candidate. We solve a set of 1D relaxation problems describing the profiles of SDWs and DDWs (see Supplementary Information S1 for more details on relaxation codes), and extract the domain wall formation ratio β. This allows us to define a band in (ζ, τ) plane of GSFE's that comply with the

β ¼ 1:90 (see Supplementary Fig. S3g). Fig. 2b compares one GSFE moiré constrained candidate with the β ¼ 1:90 band (magenta) with the well-accepted choice of GSFE of ref. 31 (blue), which notably falls outside of the band with β ¼ 1:98. The moiré metrology analysis indicates that SDWs implied by GSFE in ref. 31 insufficiently attract one another (blue curve in Fig. 2c) to account for the observed network, as indeed revealed in Fig. 2d. In contrast, the moiré-constrained candidate (magenta in Fig. 2b) with a flatter saddle point promotes stronger SDWs attraction across a broad range of domain wall orientations (magenta curve in Fig. 2c; see Supplementary Information S1 and S5 for additional details), and yields excellent agreement with the data (Fig. 2e). Regardless of the good agreement, the relatively flat saddle point comes as a surprise, and may in fact correct for an unknown effect unrelated to interlayer energy.

Moiré metrology of twisted double bilayer graphene. Compared to TBG, the TDBG system makes an even more interesting casestudy due to the small yet finite imbalance between the two lowest-energy phases: Bernal (ABAB) and rhombohedral (ABCA) stackings 17 . This imbalance results in an energy cost per-unitarea (σ) for rhombohedral relative to Bernal stackings, leading to characteristic curved domains 17,38 (see Fig. 1a,b). Exploring large areas of TDBG reveals a rich distribution of rhombohedral domain shapes (see Fig. 1c and other TDBG images in this work). Fig. 3a summarizes this distribution as a histogram of inverse curvature values (κ -1 ), extracted from images as in Fig. 1b (see Supplementary Information S6 for more examples). The histogram reveals a distinct clustering about a value of κ -1 = 440 ± 120 nm, which we use to assess the accuracy of several variants of the GSFE from available DFT functionals (Fig. 3b and see "Methods"). All reported GSFE variants are qualitatively similar to the TBG case, peaking at the BAAC configuration, and having a saddle point barrier between ABAB and ABCA. A closer inspection (inset) reveals a profound difference between the GSFEs for the ABCA relative to the ABAB that governs domain curvature. We model the domain curvature and structure by a continuous 2D relaxation code (Supplementary Information S1). Fig. 3c,d show two representative cases, with disparate outcomes. In Fig. 3c (resembling the experimental case of Fig. 1b) the energy is minimized by slight bending of the SDW into the ABCA region. As the twist angle decreases (or as strain increases as in Fig. 1c), at some point it becomes energetically beneficial to form DDWs (Fig. 3d). As the twist angle further decreases, the shape of the ABCA domains remains unchanged. Similarly, solving for the domain formation for all DFT approaches and across a wide twist angle range, we compare the extracted κ -1 . Interestingly, κ is independent of the twist angle for all GSFE variants (with values indicated by colored lines over Fig. 3a), which is not generally the case (see discussion in Supplementary Information S7). The domain structures are further captured by the 2D "soap-bubble" model, as seen in turquoise dashed lines in the representative cases of Fig. 3c, d and more generally in Supplementary Videos 2-5. This model approximates the total energy as a sum of a domain area term and two line-energy terms as

where γ 1,2 are the lineenergies of SDW and DDW as a function of the domain-wall orientation respectively, the integrations are along the domain walls, and S is the area of the domain (see Supplementary Information S7). All model parameters require only the GSFE (and elastic properties) to describe domain shapes, with no additional tuning parameters (Supplementary Information S7). One approach, DFT-D2, remarkably reproduces the experimental cluster (Fig. 3a), due to relatively high σ and comparable lineenergies to other approaches (see Supplementary Information S2). It is noteworthy that the DFT-D approach has not been previously considered as the leading approach when theoretically benchmarked against the Quantum Monte-Carlo method for the binding energy of AA and AB stacking of bilayer graphene 39 .

The rhombohedral domains represented in the histogram of Fig. 3a exemplify a well-defined electrostatic environment near charge neutrality point (CNP) in the absence of the interlayer bias. As shown recently 38 , upon charging and biasing the balance between the rhombohedral and Bernal phases can shift. An extreme demonstration of malleability of TDBG moiré patterns under a non-uniform distribution of charges and high strain conditions is presented in Fig. 3e. Three holes (marked by blue circles) punctured one of the bilayers. This procedure prompts a highly strained moiré pattern, most strongly manifested in the densely packed parallel DDWs structures connecting the two bottom holes. The stack shows strong defect-induced doping (see discussion in Supplementary Information S5), apparent in the enhanced nearfield contrast between the ABCA (dark) and ABAB (bright) phases (compare to contrast of Fig. 1c). Further support for the high non-uniformity of charge distribution is an observed region of flipped balance, where the ABAB phase becomes unstable relative to ABCA across a sharp (~50 nm) interface (to the left of the top hole). Attempting to model the moiré superlattice with the DFT-D2 GSFE at CNP (red in Fig. 3b) fails to capture the observed structure of excessively curved SDWs (color-map of Fig. 3f). However, when introducing a doping level of 8 × 10 12 cm -2 the resulting GSFE (dashed light-green in Fig. 3b) better captures the observed structure (green dots in Fig. 3f tracking the domain walls in the calculation). The difference between the two models becomes more pronounced for regions of higher strain, as highlighted in the inset of Fig. 3f (compare green and red dots in respect to the experimental map). Therefore, minuscule energy differences between models of order 0.1 meV/u.c. (inset of Fig. 3b) result in measurable spatial features of the relaxed moiré patterns. To put this figure in context, the theoretical method which is widely considered as the goldstandard of ab-initio quantum chemistry 36,37 yields an accuracy as low as 3 meV/atom 37 .

To understand the enhanced sensitivity of moiré metrology under strain (as seen in Fig. 3f), we propose an alternative description of the moiré superlattice in terms of a geometric interference pattern of the lattices of the two layers (see Supplementary Information S1). At minute twist angles, the relaxed moiré patterns are essentially a projection of the detailed energy landscape over space, accumulated over large regions compared to the atomic scale. The introduction of strain between the layers, whether naturally occurring or externally controlled, alters the interference pattern (Supplementary Information S8). As strain pushes the domain walls in Figs. 123together, it also promotes their interaction; both effects are reflected in the relaxed moiré pattern (also see Supplementary Video 6).

Moiré metrology, introduced here, correlates first principle calculations of the stacking energy function with measurable spatial features of twisted vdW systems. The stacking energy function is widely used for modeling twisted multilayers across a broad range of twist angles and strain conditions, and has direct implications for the electronic band-structure 40 . Moiré metrology is not restricted to the discussed material systems and can be universally applied with other systems of great current interest [41][42][43][44][45][46][47][48] . Therefore, by providing a reliable account of the stacking energy function, moiré metrology has a broad impact across the field of vdW heterostructures. Furthermore, the moiré metrology tools can also be used for modeling and designing non-uniform strain fields in realistic devices. Finally, due to its outstanding stacking energy sensitivity, we propose moiré metrology as a concrete experimental path to provide much needed benchmarks for first-principle theoretical approaches 49 .

Samples preparation Source crystals. The bulk crystals of MoSe 2 and WSe 2 were grown by the self-flux method 50 : single crystals of MoSe 2 were prepared by combining Mo powder (99.997%, Alfa Aesar 39686) and Se shot (99.999%, Alfa Aesar 10603) in a ratio of 1:50 (Mo:Se) in a quartz ampoule. The ampoules were subsequently sealed under vacuum (∼5 × 10 -6 Torr). The reagents were then heated to 1000 °C within 24 h and dwelled at this temperature for 8 weeks before being cooled to 350 °C over 4 weeks. At 350 °C the Se flux was decanted through alumina wool (Zircar D9202) in a centrifuge and the ampoules were quenched in air. The subsequently obtained MoSe 2 single crystals were annealed at 275 °C with a 200 °C gradient for 48 h to remove any residual Se. A similar process was also used to synthesize single crystals of WSe 2 with 1:15 (W:Se) as the starting ratio and using W powder (99.999%, Alfa Aesar 12973). Kish graphite source crystals were purchased from Graphene Supermarket.

Exfoliation. The MoSe 2 and WSe 2 monolayers, and Graphene and hBN flakes were mechanically exfoliated from the bulk single crystals onto SiO 2 /Si (285 nm oxide thickness) chips using the tape-assisted exfoliation technique (the tape used was Scotch Magic Tape). The exfoliation followed ref. 51 , such that the Si chips were treated with O 2 plasma (using a benchtop radio frequency oxygen plasma cleaner of Plasma Etch Inc., PE-50 XL, 100 W at a chamber pressure of ~215 mTorr) for 20 s for graphene, for 10 s for MoSe 2 and WSe 2 , and no O 2 plasma treatment for hBN. The chips were then matched with respective exfoliation tape. In the graphene case the chip+tape assembly were heated at 100 °C for 60 s and cooled to room temperature prior to removing the tape. Such thermal treatment was not done for other materials. The MoSe 2 and WSe 2 monolayer relative crystallographic orientation was obtained by linear-polarization-resolved second-harmonic generation (SHG).

Stack preparation. All heterostructures were assembled using standard dry-transfer techniques 52 with a polypropylene carbonate (PPC) film mounted on a transparent-tape-covered polydimethylsiloxane (PDMS) stamp. The transparent tape layer was added to the stamp to mold the PDMS into a hemispherical shape which provides precise control of the PPC contact area during assembly 53 .

All graphene heterostructures were made by first picking up the bottom-layer boron nitride (h-BN) (>25 nm thick), followed in the case of the TDBG samples by a graphite bottom gate-layer (>5 nm thick), then a dielectric BN layer (>25 nm thick). Prior to pick-up, mechanically exfoliated graphene flakes on Si/SiO 2 were separately patterned with anodic-oxidation lithography 54 to facilitate the "cut-andstack" technique 55 where it was used (all samples but those used for STM, which used the established tearing method). In the case of the TDBG samples, additional anodic-oxidation lithography was used prior to pick-up to provide additional texture to the strain landscape, e.g., cut holes inside the bulk of one of the graphene layers or non-rectilinear edge geometries.

In the case of MoSe 2 /WSe 2 the PPC was used to pick up a thin layer of exfoliated h-BN and a few layers of graphene. Then a monolayer WSe 2 was picked up and using the SHG data MoSe 2 monolayer was lifted on a rotation stage with ~1°twist angle. The stack was flipped over a Si/SiO 2 (285 nm) chip at 120 °C. In the last step, the sample was thermally annealed in a high vacuum chamber to remove the PPC at 250 °C for 1 h. Some of the presented measurements (TDBG stacks of Fig. 3e, Supplementary Figs. S8b, d-f and red bins in histogram of Supplementary Fig. S8a) were taken at this point while the stack was on a PDMS/transparent-tape/PPC structure. This provided access to the meta-stable large rhombohedral domain before they were suppressed by thermal annealing (see Supplementary Information S9).

After optional mid-assembly scanning probe measurement and/or optional encapsulation of the twisted-graphene layers, the PPC film with the heterostructure on top is mechanically removed from the transparent-tape-covered PDMS stamp and placed onto a Si/SiO 2 substrate such that the final pick-up layer is the top layer.

In the case of the TBG device presented in this work, the underlying PPC was removed by vacuum annealing at T = 350 °C. Standard plasma etching and metal deposition techniques 52 were then used to shape and make contact to the samples. In the case of the TDBG devices for STM (Fig. 1b) the stack was made using the established tearing method, using PPC as a polymer to sequentially pick up hBN, half of a piece of graphene followed by the second half with a twist angle.

For all samples for STM imaging, standard metal deposition techniques were avoided in order to maintain a pristine surface, therefore, direct contact was made to the stack by micro-soldering with Field's metal 56 , keeping temperatures below 80 C during the entire process.

Nearfield imaging techniques. In this work, we used two nano-optical imaging techniques: cryogenic nano-photocurrent imaging (used for TBG imaging) and phase-resolved scattering type scanning optical microscope imaging (used for TDBG imaging) (s-SNOM). Cryogenic photocurrent imaging 57 was done with a home-built cryogenic SNOM, and s-SNOM nearfield imaging was done with a commercial (Neaspec) SNOM 11 . In both cases using mid-IR light (continuous wave CO 2 gas laser [Access Laser] at a wavelength of 10.6 µm) focused to a diffraction limited spot at the apex of a metallic tip, while raster scanning the sample at tapping mode. Fig. 2a was acquired at a temperature of 100 K and while tuning the silicon back-gate to a relatively high doping of 3 × 10 12 cm -2 . In such a case, we observe, for the first time, a non-local photo-current generation regime. In this regime, the light induced temperature profile is broad (relative to system size) and the photo-current generation is located at a distant interface (a monolayer twisted bilayer interface in this case, clearly visible as the bright region with plasmonic fringes at the top right section of Fig. 2a). The signal contrast in such a case results from absorption contrast between different stacking configuration, and not from thermo-electric properties. This unique approach provides a high-resolution image, not limited by thermal length-scales (see Supplementary Information S4 for more details).

In the s-SNOM case, we collect the scattered light (power of 3 mW) by a cryogenic HgCdTe detector (Kolmar Technologies). The far-field contribution to the signal can be eliminated from the signal by locking to a high harmonic (here we used the 3rd harmonic of the tapping). The phase of the backscattered signal was extracted using an interferometric detection method, the pseudo-heterodyne scheme, by interfering the scattered light with a modulated reference arm at the detector. Fig. 1c was acquired with a level of 1 × 10 12 cm 2 p-doping applied with a Si backgate. Such a level of doping has negligible effect on domain curvature (see Supplementary Information S5 for more details).

STM and STS imaging. STM and STS measurements were carried out in a homebuilt STM under ultra-high vacuum conditions. MoSe 2 /WSe 2 measurements were performed at 300 K while TDBG measurements were performed at 5.7 K. The setpoints of the STM (constant current mode) imaging were V = -1.8 V and I = 100 pA for MoSe 2 /WSe 2 measurements and V = 0.5 V and I = 50 pA for TDBG measurements. The tip-sample bias for the STS measurement of Fig. 1b was -57.5 meV. STM tips were prepared and calibrated for atomic sharpness and electronic integrity on freshly prepared Au (111) crystals. Samples were measured with multiple tips to ensure consistency of results.

Parameter DFT calculations. The DFT calculations of the GSFE parameters were performed with the Vienna Ab initio Simulation Package (VASP) 58 . Plane wave basis sets were employed with energy cutoff of 1200 eV and 500 eV for the calculations for TDBG and MoSe 2 /WSe 2 , respectively. Pseudopotentials were constructed with the projector augmented wave method (PAW) 59,60 . 60 × 60 × 1 and 11 × 11 × 1 Γ-centered k-point grids were used in the calculations for TDBG and MoSe 2 /WSe 2 , respectively. Vacuum spacing larger than 15 Å was added along the z direction to eliminate the artificial interaction between periodic slab images in all the calculations. In the calculations of the GSFE, the x-y coordinates of the atoms in all the 2D layers were fixed and the z coordinates were allowed to relax until forces in the z-direction are less than 1 meV/Å for TDBG and 20 meV/Å for MoSe 2 /WSe 2 . The van der Waals interactions are important in evaluating the energetics of the 2D layer structures. We tested this effect for TDBG by considering four approaches: (1) employing PAW pseudopotentials with the exchangecorrelation functionals treated at the local density approximation (LDA) level. (2) Employing PAW pseudopotentials with the exchange-correlation functionals treated at the generalized gradient approximation (GGA) level 61 and additional van der Waals corrections are applied with the DFT-D2 method of Grimme 62 . (3) Employing PAW pseudopotentials with the GGA functionals and additional van der Waals corrections with the Tkatchenko-Scheffler method 63 . (4) Employing PAW pseudopotentials with a non-local correlation functional (optB88-vdW) [64][65][66] that approximately accounts for dispersion interactions. In the calculations for MoSe 2 /WSe 2 , we employed PAW pseudopotentials with the GGA PBE functionals with the DFT-D3 van der Waals corrections 67 . The bulk modulus and shear modulus for each material were calculated by applying isotropic or uniaxial strain to a monolayer lattice and then performing a quadratic fit to the strain-dependent energies. For TDBG, we assume the elastic coefficients are twice the values extracted for monolayer graphene. The DFT ground state energy is evaluated on a regular grid of different interlayer configurations. The Fourier components of the resulting energies are then extracted to create a convenient functional form for the GSFE.