# Fileset

[s41467-022-35213-5.pdf](https://mdr.nims.go.jp/filesets/f6bb338b-18b5-40a4-b1a2-8f19150b5a57/download)

## Creator

Dorri Halbertal, Simon Turkel, Christopher J. Ciccarino, Jonas B. Hauck, Nathan Finney, Valerie Hsieh, [Kenji Watanabe](https://orcid.org/0000-0003-3701-8119), [Takashi Taniguchi](https://orcid.org/0000-0002-1467-3105), James Hone, Cory Dean, Prineha Narang, Abhay N. Pasupathy, Dante M. Kennes, D. N. Basov

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Unconventional non-local relaxation dynamics in a twisted trilayer graphene moiré superlattice](https://mdr.nims.go.jp/datasets/d18ec3ca-c2ea-4090-b6b0-2708a4e27753)

## Fulltext

Unconventional non-local relaxation dynamics in a twisted trilayer graphene moiré superlatticeArticle https://doi.org/10.1038/s41467-022-35213-5Unconventional non-local relaxationdynamics in a twisted trilayer graphenemoiré superlatticeDorri Halbertal 1,8 , Simon Turkel 1,2,8, Christopher J. Ciccarino 3,Jonas B. Hauck 4, Nathan Finney1, Valerie Hsieh1, Kenji Watanabe 5,Takashi Taniguchi 6, James Hone 1, Cory Dean 1, Prineha Narang 3,Abhay N. Pasupathy 1,2, Dante M. Kennes 4,7 & D. N. Basov 1The electronic and structural properties of atomically thin materials can becontrollably tuned by assembling them with an interlayer twist. During thisprocess, constituent layers spontaneously rearrange themselves in search of alowest energy configuration. Such relaxation phenomena can lead to unex-pected and novel material properties. Here, we study twisted double trilayergraphene (TDTG) using nano-optical and tunneling spectroscopy tools. Wereveal a surprising optical and electronic contrast, as well as a stacking energyimbalance emerging between themoiré domains.We attribute this contrast toan unconventional form of lattice relaxation in which an entire graphene layerspontaneously shifts position during assembly, resulting in domains of ABA-BAB and BCBACA stacking. We analyze the energetics of this transition anddemonstrate that it is the result of a non-local relaxation process, in which anenergy gain in one domain of the moiré lattice is paid for by a relaxation thatoccurs in the other.The discovery of superconductivity in rotationally misalignedgraphene bilayers established moiré engineering as a robustway to create strongly correlated phases in van der Waalsheterostructures1,2. Since this initial discovery, a wide range ofmoiré materials have emerged with fascinating electronic proper-ties such as correlated insulators3–8, strange metals9,10, electronicnematics11–14, and Wigner crystals15, among other unconventionalphases. Moiré materials additionally hold promise as platforms fortunable quantum simulation, enabling future applications in mate-rials discovery and design16. With the advent of new and morecomplexmoiré device geometries, with greater than two layers5–8 orgreater than one twist angle between layers17–20, it becomesincreasingly important to consider the effects of lattice relaxation,or the spontaneous rearrangement of atoms in search of a lowerenergy configuration, on the final microscopic crystal structure. Inmirror symmetric twisted trilayer graphene, for instance, it wasrecently observed21 that lattice relaxation leads to the emergence ofmoiré defects that are not observed in the simpler twisted bilayersystem. A fuller understanding of relaxation phenomena in moiréheterostructures thus holds the potential to enable exploitation ofthese effects as a means of engineering novel or otherwise unstablematerial systems.Received: 29 July 2022Accepted: 18 November 2022Check for updates1Department of Physics, Columbia University, New York, NY 10027, USA. 2Condensed Matter Physics and Materials Science Division, Brookhaven NationalLaboratory, Upton, NY 11973, USA. 3Harvard John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA.4Institute for Theory of Statistical Physics, RWTH Aachen University, and JARA Fundamentals of Future Information Technology, 52062 Aachen, Germany.5Research Center for Functional Materials, National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan. 6International Center for MaterialsNanoarchitectonics, National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan. 7Max Planck Institute for the Structure and Dynamics ofMatter, Center for Free Electron Laser Science, Hamburg, Germany. 8These authors contributed equally: Dorri Halbertal, Simon Turkel.e-mail: dorrihal@gmail.comNature Communications |         (2022) 13:7587 11234567890():,;1234567890():,;http://orcid.org/0000-0002-1665-0724http://orcid.org/0000-0002-1665-0724http://orcid.org/0000-0002-1665-0724http://orcid.org/0000-0002-1665-0724http://orcid.org/0000-0002-1665-0724http://orcid.org/0000-0002-5011-7579http://orcid.org/0000-0002-5011-7579http://orcid.org/0000-0002-5011-7579http://orcid.org/0000-0002-5011-7579http://orcid.org/0000-0002-5011-7579http://orcid.org/0000-0002-3195-4777http://orcid.org/0000-0002-3195-4777http://orcid.org/0000-0002-3195-4777http://orcid.org/0000-0002-3195-4777http://orcid.org/0000-0002-3195-4777http://orcid.org/0000-0003-3399-1341http://orcid.org/0000-0003-3399-1341http://orcid.org/0000-0003-3399-1341http://orcid.org/0000-0003-3399-1341http://orcid.org/0000-0003-3399-1341http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-8084-3301http://orcid.org/0000-0002-8084-3301http://orcid.org/0000-0002-8084-3301http://orcid.org/0000-0002-8084-3301http://orcid.org/0000-0002-8084-3301http://orcid.org/0000-0003-2967-5960http://orcid.org/0000-0003-2967-5960http://orcid.org/0000-0003-2967-5960http://orcid.org/0000-0003-2967-5960http://orcid.org/0000-0003-2967-5960http://orcid.org/0000-0003-3956-4594http://orcid.org/0000-0003-3956-4594http://orcid.org/0000-0003-3956-4594http://orcid.org/0000-0003-3956-4594http://orcid.org/0000-0003-3956-4594http://orcid.org/0000-0002-2744-0634http://orcid.org/0000-0002-2744-0634http://orcid.org/0000-0002-2744-0634http://orcid.org/0000-0002-2744-0634http://orcid.org/0000-0002-2744-0634http://orcid.org/0000-0002-9838-6866http://orcid.org/0000-0002-9838-6866http://orcid.org/0000-0002-9838-6866http://orcid.org/0000-0002-9838-6866http://orcid.org/0000-0002-9838-6866http://orcid.org/0000-0001-9785-5387http://orcid.org/0000-0001-9785-5387http://orcid.org/0000-0001-9785-5387http://orcid.org/0000-0001-9785-5387http://orcid.org/0000-0001-9785-5387http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-022-35213-5&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-022-35213-5&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-022-35213-5&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-022-35213-5&domain=pdfmailto:dorrihal@gmail.comTwisted double trilayer graphene (TDTG), a moiré material thathas not yet been experimentally investigated, is a natural next step inextending the moiré paradigm to more complex structures, in whichlattice relaxation can lead to unexpected atomic configurations. TDTGis formed in a manner analogous to twisted bilayer (TBG) and twisteddouble bilayer graphene (TDBG), namely by stacking two Bernal tri-layers from the same source crystal with a small relative twist. In thelow twist angle limit, moiré patterns in few-layer graphene (FLG) sys-tems generate large domain structures with distinct crystallographicstackings22,23. If rigidly stacked in a manner that preserves the ABAstacking of each trilayer (referred to here as the rigid scenario), TDTGforms domains of mixed rhombohedral and Bernal character withABABCB and BCBABA stackings (Fig. 1a), where each domain containsa unit of three rhombohedrally stacked layers (3R).It is conceivable, however, that stresses and strains applied to asample during the fabricationprocess canact as aneffective annealing,allowing the system to explore other stacking configurations beforereaching the lowest energy equilibrium. In the case of TDTG, applyinga simple translation to the second layer (Fig. 1b) results in one domain(ABABAB) with pure Bernal stacking and another domain (BCBACA)with a unit of four rhombohedrally stacked layers (4R). This layer slidescenario provides an interesting test case from an energy standpoint.While the pure Bernal phase is energetically favorable, its realizationcomes at the price of increasing the local stacking energy of themixedrhombohedral domain by transforming it from a 3R to a 4R config-uration. Prior studies of relaxation effects in twisted van der Waalsheterostructures22,24–27 have focused on processes in which latticerelaxation acts to uniformly reduce the stacking energy at everylocation in space, generally by minimizing the area of higher energydomains and maximizing the area of their lower energy counterparts,as occurs in minimally twisted TDBG. Consideration of the scenariospresented in Fig. 1a, b for TDTG raises the question of whetherrelaxation processes can likewise act to offset an increase of thestacking energy in one area of a devicewith a decrease in an area that isas far as several microns away. Such relaxation at a distance wouldoffer a means of stabilizing elusive atomic configurations with distinctelectronic properties by structurally coupling them to low-energyphases through moiré patterning.ResultsIn this work, we utilize mid-infrared scanning near-field opticalmicroscopy (SNOM) and scanning tunneling microscopy (STM) andspectroscopy (STS) to characterize the optical and electronicproperties of TDTG samples. The use of SNOM and STM/S onidentical samples allows us to characterize the electronic propertiesof a device over both large (microns) and small (nanometers) areas,giving direct experimental access to the interplay of length scalesthat is crucial to non-local relaxation dynamics. Figure 1c, d showsphase resolved SNOM images28–30 taken over a ~10 μm region of theTDTG device shown in the inset of Fig. 1c. In the rigid scenario forTDTG (Fig. 1a), we expect the local stacking configuration in eachdomain to be ABABCB and BCBABA, which are related to each otherby inversion across the x–y plane (Cz2). In the absence of external Cz2symmetry breaking mechanisms, such as displacement field or anasymmetric dielectric environment, these two configurations aretherefore expected to possess identical electronic and structuralproperties. Contrary to this expectation, a moiré superlattice withclear optical contrast between adjacent domains is observed both inthe amplitude (Fig. 1c) and phase (Fig. 1d) of the near-field signal.Furthermore, there is a clear imbalance in the stacking energy of thetwo domains, as evidenced by the convexity (concavity) of thebright (dark) triangles in Fig. 1d. This difference in stacking energy isparticularly clear in regions of large heterostrain, such as in the topleft of Fig. 1d, where a linear pattern is observed, corresponding toFig. 1 | Moiré superlattice of twisted double trilayer graphene (TDTG). a Thelattice structure of the two lowest energy stacking configurations of TDTG whenassuming both top and bottom trilayer graphene (TLG) are Bernal stacked. Theinterfacial atoms are colored to highlight the AB/BA stacking configurations. b TheTDTG system after a global slide of the middle layer of one of the TLG sheets. Thedirection of slide is marked by red arrows indicating the transition from an ABA(dashed red honeycomb lattice) to a BAB configuration. Such a slide, if realized,results in the formation of the unstable BCBACA phase. c, d Mid-IR near-fieldamplitude (c) and phase (d) of the TDTG stack. Bright (convex) and dark (concave)triangular domains are observed across the sample. Near-field signal in c, d is the4th harmonic demodulation of the probe tapping frequency (laser illumination at1000 cm−1, see “Methods” for additional experimental details). Inset of c: Opticalimage of the TDTG/hBN stack. The bottom and top TLGs are highlighted by purpleand green frames, respectively. The red rectangle marks the scan area of (c, d).Article https://doi.org/10.1038/s41467-022-35213-5Nature Communications |         (2022) 13:7587 2double domain walls (DDWs) that emerge due to the collapse ofunstable domains22.Figure 2 examines the TDTGmoiré superlattice at the sub-micronlength scale. The complex near-field signal is presented in Fig. 2a–das afunction of tapping-probe demodulation harmonic (focusing on thegreen square in Fig. 1c, d). These finer scans further demonstrate theoptical contrast between adjacent domains aswell as the curvingof thedomain walls (DWs). Additional details emerge in these higher-resolution images as well, including clear bright features along theDWs and bright spots at the DW intersections in Fig. 2d. Furthermore,we find that the complex near-field optical contrast between thedomains is a monotonically increasing function of demodulation har-monic (Fig. 2e), indicating that the signal is coming from the graphenedevice rather than the substrate. We examine the microscopic elec-tronic structure in greater detail by performing STM/S measurementsover regions of the same TDTG sample. The local density of states(LDOS) measured near the Fermi level is displayed in Fig. 2f, whichshows a striking contrast in tunneling conductivity between differentdomains of the moiré lattice. Characteristic spectra acquired on eachof the two domains are plotted in Fig. 2g, demonstrating that theobserved LDOS contrast is caused by a large spectral peak near zeroenergy that is present in only the convex domain. The spectral shape inthe concave domain, on the other hand, is largely featureless at lowenergy and possesses no comparable peak. Moreover, measurementsof the tunneling spectrum on the untwisted region (Fig. 2h) indicatethat the source crystal is Bernal trilayer, confirming that the moirécontrast is a property of the six-layer system. This suite of measure-ments unambiguously demonstrates a significant difference inelectronic structure between each of the observed moiré superlatticedomains. In the Supplementary Information (section 3) we considerand rule out alternative sources of Cz2 symmetry breaking, includingthe effect of the hBN substrate and the possibility of atomic-scale near-field tomography, i.e., the breaking ofCz2 symmetry by the sharpprobeinteracting with individual atomic layers. The naive expectation of therigid scenario (Fig. 1a) is therefore clearly not realized in our TDTGdevice.Applying a global translation to one of the six layers in a TDTGheterostructure has the potential to lower the overall stacking energyof a device even as itmight raise the energy density in certain confinedregions.Oncewe accept that the energy barrier to such a transition canbe overcome (see below), there is in principle no reason to restrict ouranalysis to the particular layer slide scenario depicted in Fig. 1b. In aneffort to match the experimental observations, we have thereforeperformed DFT calculations of the band structure, density of states,and stacking energy density of all 25 possible TDTG stacking config-urations (eachof the five-layer interfaces can be stacked as either ABorBA). Eight of these configurations describe, in the minimally twistedregime, moiré pairs that are related by Cz2 symmetry (see Supple-mentary Information section 4), which we have already ruledout above.Figure 3 explores the remaining twenty-four TDTG configurationswith crystallographically distinctmoiré domains. These can be dividedinto four groups (corresponding to the four columns of Fig. 3) basedon their symmetries. Moiré pairs in Fig. 3a–d are connected by blackhorizontal lines (C2 symmetry pairs are connected by colored lines asindicated). Each possible moiré domain is characterized by itsFig. 2 | Imaging of the TDTG moiré. a–d Mid-IR near-field imaging (over greensquare in Fig. 1a) showing bright (convex) and dark (concave) domains: normalizedamplitude (i) and phase (ii) at different probe demodulation harmonics from 2nd(a) to 5th (d). The measurement was done at room temperature and a laser fre-quency of 1000 cm−1. a–d share a color-bar and scale-bar. Dashed arcs in (d) with aradius of 850 nm highlight the domain shape. e Amplitude (top) and phase (bot-tom) contrasts between the two domains of (a–d) as a function of demodulationharmonics, revealing a monotonic trend. f Differential conductance measured bySTS, acquired at a sample bias of −20meV.gTunneling spectrummeasured in eachof the two domains (marked by correspondingly colored dots in (f)). The dashedline marks the energy at which (f) was acquired. Each curve is normalized to itsvalue at 500 mV. h Tunneling spectrum measured on exposed TLG section (greendot in Fig. 1a) indicative of a Bernal stacked TLG source crystal.Article https://doi.org/10.1038/s41467-022-35213-5Nature Communications |         (2022) 13:7587 3calculated Fermi level spectral weight and stacking energy density(reflected in Fig. 3a–d by vertex color and size, respectively). Moirépairs with a large difference in stacking energy result in curveddomains, similar to those seen experimentally, as confirmed by atomicrelaxation calculations (Fig. 3e–h). In seeking a match with ourexperimental observations, we therefore require a pair of moirédomainswith a largedifference in both stacking energy andFermi levelspectral weight tomatch the domain curvature and electronic contrastrevealed by SNOMand STM.While three groups ofmoiré pairs possesssufficient domain curvature (Fig. 3f–h), only the ABABAB/BCBACAconfiguration displays calculated spectra that are consistent with ourSTS results (Fig. 3k). The concave domain of this pair (BCBACA) showsa peak at low energy where the convex domain (ABABAB) remainsfeatureless, as in experiment. In addition, the sharp steps at ~± 350mVin the BCBACA domain, which are associated with the edges of elec-tronic bands, align quantitatively with similar steps observed experi-mentally in the concave domain (compare Fig. 2g). This excellentmatch of the calculated electronic structurewith themeasured densityof states spectrum points conclusively to ABABAB/BCBACA as thestacking configurations spontaneously realized in our TDTG device.DiscussionConstructing an ABABAB/BCBACA stacking configuration from aminimally twisted Bernal trilayer source crystal requires a global slid-ing of themiddle layer of one of the two twisted trilayers (illustrated inFig. 3k, inset). The energetics of sucha global translation are seeminglycounter-intuitive because, first, it requires a large energy input to thesystem to realize a universal layer shift, and second, that shift results ina local increase in the stacking energy density of one of the twodomains. Sliding of a graphene layer can be viewed as continuouslytraversing the stacking energy landscape shown in Fig. 4a. To trans-form adjacent layers from AB to BA stacking, as required to realize theexperimentally observed configuration, the system must cross a for-midable energy barrier of ∼6× 104 eVμm2 set by the saddle-point (SP) ofthe stacking energy function (indicated by SP in Fig. 4a). This cannotconceivably be overcome by thermal excitation alone. The only step inour experiment during which the sample is subjected to forces ofsufficient magnitude to induce a sliding transformation is the stackingprocess, which involves pressing together each constituent trilayer ofthe TDTG device before peeling them away from the exfoliation sub-strate (Fig. 4b).When stacking induced sliding transitions like this havebeen observed in the past31, they have as a rule been from ametastable(rhombohedral) to a stable (Bernal) phase. In our case, however, thetransition from ABABCB/BCBABA (3R/3R) to ABABAB/BCBACA (0R/4R) acts to decrease the thermodynamic stability of nearly half of thedevice area.Lattice relaxation in TDTG therefore takes an unconventionalform in which an energy gain in one half of the crystal is paid for by arelaxation process that occurs in the other. The energy justification forthis phenomenon is studied in Fig. 4c, where each curve represents thetotal energy density (stacking and elastic energy) of a givenmoiré pair,after atomic relaxation, as a function of moiré wavelength λ. All curvesFig. 3 | Exploration of candidates for TDTG moiré superlattice structures.a–d Each panel addresses a group of configurations. Each moiré pair is connectedby a horizontal black line. Configurations that are C2 symmetry pairs are connectedby blue (Cx2), green (Cz2), andmagenta (Cx2Cz2) lines, respectively. Each configurationis marked by a circle whose color indicates the low energy spectral weight, andwhose size indicates the configuration’s stacking energy density (see legend).e–h The domain formation is reflected by the stacking energy density for a twistangle of 0.04∘ for each moiré pair. Each panel is the result of the atomic relaxationcalculation (see “Methods”) using the DFT calculated energy imbalance of thecorrespondingmoirépair (triangles indicate thephaseswithmatching colors to theconfiguration text in (a–d)). i–l DFT calculated electronic densities of states fordifferentmoiré pairs (see “Methods”). The inset shows the lattice structure for eachconfiguration (with consistent colors as in (a–d)). The red arrow indicates therequired global layer sliding in order to realize the particular moiré superlatticefrom the rigid ABABCB/BCBABAconfiguration. The configurationswhere themoirépairs are also C2 symmetry pairs were omitted here for brevity, as they could notproduce an energy imbalance (these missing configurations are explored in Sup-plementary Information section 4).Article https://doi.org/10.1038/s41467-022-35213-5Nature Communications |         (2022) 13:7587 4are referenced to the energy of the rigid configuration (ABABCB/BCBABA), marked by Eno−slide. For small λ (and therefore weak relaxa-tion) the ABABCB/BCBABA and ABABAB/BCBACA configurations areenergetically equivalent, indicating that a layer slide transition withoutadditional relaxation cannot reduce the global energy. As the moiréperiod increases, the relaxation strengthens, and the energeticimbalance between ABABAB and BCBACA domains (absent in theABABCB/BCBACA configuration) drives the expansion of the Bernaldomain, thus reducing the global energy of the ABABAB/BCBACAconfiguration relative to its rigid counterpart. If provided sufficientshear forces to overcome the SP energy barrier, the atomic relaxationprocess can therefore create and stabilize the formation of theotherwise unstable BCBACA phase simply by reducing its relativevolume fraction.We visualize the dynamics of this non-local relaxation in Fig. 4d,where we plot the instantaneous solution to the energy minimizationproblem for each stacking considered in Fig. 4c as a function ofiteration number within the steepest-descent optimization algorithm,which can be interpreted as an effective time coordinate, t. At t =0,before any relaxation has taken place, all configurations have similarenergies. As the systemsflowdown their respective energy landscapes,large domains of uniform stacking are formed, separated by straightDWs. In this intermediate regime, the rigid (ABABCB/BCBABA) andlayer slide (ABABAB/BCBACA) scenarios have equivalent energies.Only when the relaxation process has reached a point where the DWsbegin to curve, with the lower energy Bernal phase pushing into themetastable BCBACA configuration, does a layer slide transitionbecome energetically favorable (see crossing of the black and orangelines in Fig. 4d). Minimally twisted TDTG thus spontaneously seeks asolution inwhich a transition to a locallymetastable phase (BCBACA) isenabled by a shared phase boundary with a proximate stable structure(ABABAB).The development of new and increasingly complex moiré het-erostructures demands a revisiting of some of the basic assumptionsof van der Waals engineering. It is sometimes convenient to think oflayered materials as immutable building blocks that can be exfoliatedand stacked at will. In reality, however, van derWaalsmaterials inhabita complicated energy landscape thatmustbe carefully navigatedwhendesigning new device architectures. Our measurements of minimallytwisted TDTG reveal a surprising crystallographic transformation thatoccurs during the stacking process. The mechanism underlying thistransition involves a non-local energy balancing that enables the for-mation of rhombohedral domains by their coupling to a simulta-neously formed relaxed Bernal structure. This has immediateFig. 4 | Energy analysis for atomic relaxation driven formation of unstablephases by sliding layers. a Stacking energy density at the interface between twographene sheets for different configurations. The plot presents cuts along thezigzag and armchair directions. Inset: Stacking energy density in the two-dimensional configuration-space reflecting the stacking energy density as a func-tion of translation of one sheet relative to the other (green honeycomb lattice).b Scenarios for generation of a new phase during stacking. Top: New phase(orange) generation during contact of the two TLGs (blue). Bottom: New phasegeneration (orange) as the TDTG stack is picked up from the surface. c Comparingtotal energy density as a function of moiré periodicity for different relaxed moirésuperlattices (see “Methods” for details on atomic relaxationmodels). All curves arereferenced to the rigid scenario, ABABCB/BCBABA. ABABAB/BCBACA becomesincreasingly energetically favorable as themoiré period increases. d Simulating theABABCB/BCBABA to ABABAB/BCBACA phase stability inversion through theatomic relaxation process. Each curve shows the total energy density for a parti-cular phase at different optimization steps of the gradient-descent process. Above:instantaneous stacking energy densities of the lowest energy moiré system atrepresentative steps. A sharp transition between favorable ABABCB/BCBABA(black) to favorable ABABAB/BCBACA (orange) is observed as the relaxationprogresses.Article https://doi.org/10.1038/s41467-022-35213-5Nature Communications |         (2022) 13:7587 5implications both for research andpractical applications, as it reveals aheretofore unexplored factor in determining the final structure offabricated devices. As device geometries become increasingly com-plex, non-local relaxation phenomena are likely to play an importantrole in device design. Detailed knowledge of this and similar relaxationprocesses holds the potential to utilize the power of lattice relaxationfor engineering novel and otherwise unstable material systems.MethodsSamples preparationExfoliation. Graphene and hBN flakes were mechanically exfoliatedfrom the bulk single crystals onto SiO2/Si (285 nm oxide thickness)chips using the tape-assisted exfoliation technique (the tape used wasScotchMagic Tape). The Si chips were treated with O2 plasma (using abenchtop radio frequency oxygen plasma cleaner of Plasma Etch Inc.,PE-50 XL, 100 W at a chamber pressure of 215 mTorr) for 20 s forgraphene and no O2 plasma treatment for hBN. The chips were thenmatched with respective exfoliation tape. In the graphene case, thechip+tape assembly were heated at 100 ∘C for 60 s and cooled to roomtemperature prior to removing the tape. Such thermal treatment wasnot done for hBN.Stack preparation. The heterostructure was assembled using stan-dard dry-transfer techniques32 with a polypropylene carbonate (PPC)film mounted on a transparent-tape-covered polydimethylsiloxane(PDMS) stamp. The transparent tape layer was added to the stamp tomold the PDMS into a hemispherical shape which provides precisecontrol of the PPC contact area during assembly33. Theheterostructurewas made by first picking up the hBN (20-nm thick). Prior to pick-up,mechanically exfoliated TLG flakes on Si/SiO2 were separately pat-terned with anodic-oxidation lithography34 to facilitate the “cut-and-stack” technique35. Next, the PPC film with the heterostructure on topis mechanically removed from the transparent-tape-covered PDMSstamp and placed onto a Si/SiO2 substrate such that the final pick-uplayer is the top layer. Then theunderlying PPCwas removedby vacuumannealing at 350 ∘C.Near-field imaging techniquesThe mid-IR near-field scans in this work were acquired with a phase-resolved scattering type scanning optical microscope imaging (s-SNOM) with a commercial system (Neaspec), using a mid-IR quantumcascade laser (HedgehogbyDaylight Solutions) tunedbetween8.7 and10.2μm. The laser light was focused to a diffraction limited spot at theapex of a metallic tip, while raster scanning the sample at tappingmode.We collect the scattered light (power of 3–5mW) by a cryogenicHgCdTedetector (KolmarTechnologies). Thenear-field amplitude andphase were extracted as harmonic components of the tapping fre-quency using an interferometric detection method, the pseudo-heterodyne scheme, by interfering the scattered light with a modu-lated reference arm at the detector28. The near-field scans of Figs. 1 and2 were taken at 983 and 1000 cm−1, respectively.Scanning tunneling microscopy and spectroscopySTM/S measurements were conducted in a home-built STM underultra-high vacuum at 7 K. The tungsten tunneling tip was electro-chemically etched and calibrated against the Au(111) surface state priorto each sample approach. Spectroscopy was measured using a lock-inamplifier to record the differential conductance with a bias modula-tion between 1 and 7mVat 927Hz, a set point voltage of 250mV, and aset point current of 120 pA.Electronic structure theory calculations of generalized stackingfault energy function (GSFE) and DOSIn order to capture the generalized stacking fault energy function(GSFE) and electron density of states (DOS), we rely on densityfunctional theory calculations. The different six-layer graphene stack-ing configurations were captured within a hexagonal unit cell with anin-plane lattice constant of a = b = 2.459 Å. We describe the systemusing a 24 × 24 × 1 k-point mesh within the plane-wave code JDFTx36.Fermi smearing of width 0.01 Hartree is applied to the electronicoccupations. We use ultrasoft pseudopotentials37 and the PBEsolexchange-correlation functional38. In order to remove any artificialinteractions between periodic images in the out-of-plane direction, weuse a Coulomb truncation technique39. The plane-wave cutoff used is40Hartrees.We use a stringent charge density cutoff of 1000Hartreesin order to densely sample the z direction of the unit cell, which isimportant for accurately describing the energetics of the differentstacking configurations and therefore for comparison among them.Van der Waals interactions between the graphene layers are modeledusing the DFT-D3 scheme40.Using these calculations as starting points, we can then capture theelectronic density of states. We describe the electronic states of oursystems using a real-spaceWannier representation based onmaximally-localized Wannier functions41. This allows us to sample the electronicenergies at arbitrary wave vectors. In our DOS calculations, we sample5.76 × 107 wave vectors to accurately converge the DOS. We use a Lor-entzian with a broadening of width 4.3 meV to smooth the results.Atomic relaxation calculationsModeling of the atomic relaxation of TDTG structures was performedwithin a continuity model22,42. In this model, the total energy of thesystem is taken as the sumof elastic energy and a stacking energy term.The total energy wasminimized in search for the inter-layer real spacedisplacement field corresponding to the relaxed structure. Thestacking configuration at the interface was imposed to be AA at thefour corners of the moiré unit-cell. The mechanical relaxation para-meters (bulk and shear moduli) for TLG as well as the generalizedstacking fault energy function (GSFE) for the TLG/TLG interface werecalculated using DFT (see DFT section in “Methods” for details). Theresulting mechanical coefficients for TLG (in meV per unit-cell) are:bulk modulus—K = 210,971, shear modulus—G = 151,580.The GSFE coefficients were extracted from a 7 × 7 sampling of theconfiguration between two ABA-TLG, with the vertical positions of theatoms relaxed at each configuration. The Fourier components of theresulting energies were then extracted to create a convenient func-tional form for the GSFE used to describe the stacking energy term inthe atomic relaxation calculations. For simplicity, the GSFE for con-figurations other than the nominal case (ABABCB/BCBABA) used theextracted GSFE for the nominal case while imposing the stackingenergies at the lowest energy configurations as calculated by DFT(values are detailed in Supplementary Table 1). The GSFE for a givenconfiguration was taken as the closest curve (L2 norm) with the samefunctional structure, that satisfies the imposed lowest energy config-uration. This approach circumvented the need to calculate the fullstacking landscape for all systems. The validity of this approach wasassessed comparing the resulting GSFE with the full GSFE calculationfor BCABCA/CABABC yielding similar results. The resulting GSFEcoefficients are detailed in Supplementary Table 2.Data availabilityThe data used in this study are available in the Harvard Dataversedatabase [https://doi.org/10.7910/DVN/VDKZMA].References1. Cao, Y. et al. Correlated insulator behaviour at half-filling in magic-angle graphene superlattices. Nature 556, 80–84 (2018).2. Cao, Y. et al. Unconventional superconductivity in magic-anglegraphene superlattices. Nature 556, 43–50 (2018).3. Xie, Y. et al. Fractional Chern insulators in magic-angle twistedbilayer graphene. Nature 600, 439–443 (2021).Article https://doi.org/10.1038/s41467-022-35213-5Nature Communications |         (2022) 13:7587 6https://doi.org/10.7910/DVN/VDKZMA4. Nuckolls, K. P. et al. Strongly correlated Chern insulators in magic-angle twisted bilayer graphene. Nature 588, 610–615 (2020).5. Liu, X. et al. Tunable spin-polarized correlated states in twisteddouble bilayer graphene. Nature 583, 221–225 (2020).6. Cao, Y. et al. Tunable correlated states and spin-polarized phases intwisted bilayer–bilayer graphene. Nature 583, 215–220 (2020).7. He,M. et al. Symmetry breaking in twisteddouble bilayer graphene.Nat. Phys. 17, 26–30 (2021).8. Burg, G. W. et al. Correlated insulating states in twisted doublebilayer graphene. Phys. Rev. Lett. 123, 197702 (2019).9. Ghiotto, A. et al. Quantum criticality in twisted transition metaldichalcogenides. Nature 597, 345–349 (2021).10. Jaoui, A. et al. Quantum critical behaviour in magic-angle twistedbilayer graphene. Nat. Phys. https://doi.org/10.1038/s41567-022-01556-5 (2022).11. Rubio-Verdú, C. et al.Moiré nematic phase in twisteddouble bilayergraphene. Nat. Phys. 18, 196–202 (2022).12. Samajdar, R. et al. Electric-field-tunable electronic nematic order intwisted double-bilayer graphene. 2D Mater. 8, 034005 (2021).13. Kerelsky, A. et al. Maximized electron interactions at the magicangle in twisted bilayer graphene. Nature 572, 95–100 (2019).14. Jiang, Y. et al. Charge order and broken rotational symmetry inmagic-angle twisted bilayer graphene. Nature 573, 91–95 (2019).15. Li, H. et al. Imaging two-dimensional generalized Wigner crystals.Nature 597, 650–654 (2021).16. Kennes, D. M. et al. Moiré heterostructures as a condensed-matterquantum simulator. Nat. Phys. 17, 155–163 (2021).17. Hao, Z. et al. Electric field–tunable superconductivity in alternating-twist magic-angle trilayer graphene. Science 371, 1133–1138(2021).18. Park, J. M., Cao, Y., Watanabe, K., Taniguchi, T. & Jarillo-Herrero, P.Tunable strongly coupled superconductivity in magic-angle twis-ted trilayer graphene. Nature 590, 249–255 (2021).19. Siriviboon, P. et al. A new flavor of correlation and super-conductivity in small twist-angle trilayer graphene. Preprint athttps://arxiv.org/abs/2112.07127 (2021).20. Lin, J.-X. et al. Zero-field superconducting diode effect in small-twist-angle trilayer graphene. Nat. Phys. 18, 1221–1227 (2022).21. Turkel, S. et al. Orderly disorder in magic-angle twisted trilayergraphene. Science 376, 193–199 (2022).22. Halbertal, D. et al. Moirémetrology of energy landscapes in van derWaals heterostructures. Nat. Commun. 12, 1–8 (2021).23. Kerelsky, A. et al. Erratum: Moiréless correlations in ABCA gra-phene. Proc. Natl Acad. Sci. USA 118, e2017366118 (2021).24. Nam, N. N. T. & Koshino, M. Lattice relaxation and energy bandmodulation in twisted bilayer graphene. Phys. Rev. B 96,075311 (2017).25. Haddadi, F., Wu, Q., Kruchkov, A. J. & Yazyev, O. V. Moiré flat bandsin twisted double bilayer graphene. Nano Lett. 20,2410–2415 (2020).26. Guinea, F. & Walet, N. R. Continuum models for twisted bilayergraphene: effect of lattice deformation and hopping parameters.Phys. Rev. B 99, 205134 (2019).27. Moore, S. L. et al. Nanoscale lattice dynamics in hexagonal boronnitride moiré superlattices. Nat. Commun. 12, 5741 (2021).28. Sunku, S. S. et al. Photonic crystals for nano-light inmoiré graphenesuperlattices. Science 362, 1153–1156 (2018).29. Hesp, N. C. et al. Observation of interband collective excitations intwisted bilayer graphene. Nat. Phys. 17, 1162–1168 (2021).30. Jiang, L. et al. Soliton-dependent plasmon reflection at bilayergraphene domain walls. Nat. Mater. 15, 840–844 (2016).31. Yang, Y. & Al, E. Stacking order in graphite films controlled by vander Waals technology. Nano Lett. 19, 8526–8532 (2019).32. Wang, L. et al. One-dimensional electrical contact to a two-dimensional material. Science 342, 614–618 (2013).33. Kim, K. et al. Van der Waals heterostructures with high accuracyrotational alignment. Nano Lett. 16, 1989–1995 (2016).34. Li, H. et al. Electrode-free anodic oxidation nanolithography of low-dimensional materials. Nano Lett. 18, 8011–8015 (2018).35. Saito, Y., Ge, J., Watanabe, K., Taniguchi, T. & Young, A. F. Inde-pendent superconductors and correlated insulators in twistedbilayer graphene. Nat. Phys. 16, 926–930 (2020).36. Sundararaman, R. et al. JDFTx: software for joint density-functionaltheory. SoftwareX 6, 278–284 (2017).37. Garrity, K. F., Bennett, J. W., Rabe, K. M. & Vanderbilt, D. Pseudo-potentials for high-throughput DFT calculations. Computat. Mater.Sci. 81, 446–452 (2014).38. Perdew, J. P. et al. Restoring the density-gradient expansion forexchange in solids and surfaces. Phys. Rev. Lett. 100,136406 (2008).39. Sundararaman, R. & Arias, T. A. Regularization of the Coulombsingularity in exact exchange by Wigner-Seitz truncated interac-tions: towards chemical accuracy in nontrivial systems. Phys. Rev. B87, 165122 (2013).40. Grimme, S., Antony, J., Ehrlich, S. & Krieg, H. A consistent andaccurate ab initio parametrization of density functional dispersioncorrection (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 132,154104 (2010).41. Marzari, N. & Vanderbilt, D. Maximally localized generalized Wan-nier functions for composite energy bands. Phys. Rev. B 56,12847–12865 (1997).42. Carr, S. et al. Relaxation and domain formation in incommensuratetwo-dimensional heterostructures. Phys. Rev. B 98, 224102 (2018).AcknowledgementsNano-imaging research at Columbia is supported by DOE-BES grant DE-SC0018426. STMmeasurements were supported by the Office of BasicEnergy Sciences, Materials Sciences and Engineering Division, U.S.Department of Energy (DOE) under Contract No. DE-SC0012704. ANPacknowledges salary support from the National Science Foundation viagrant DMR-2004691. The development of nano-optical methods issupported as part of Programmable Quantum Materials, an EnergyFrontier Research Center funded by the U.S. Department of Energy(DOE), Office of Science, Basic Energy Sciences (BES), under award DE-SC0019443. Research on atomic relaxation is supported byW911NF2120147.Work byC.J.C. and P.N. was primarily supported by theDepartment of Energy, Photonics at Thermodynamic Limits EnergyFrontier Research Center, under Grant No. DE-SC0019140. Weacknowledge funding by the Deutsche Forschungsgemeinschaft (DFG,GermanResearch Foundation) under RTG 1995 andRTG2247,within thePriority Program SPP 2244 “2DMP”, under Germany’s Excellence Strat-egy - Cluster of Excellence Matter and Light for Quantum Computing(ML4Q) EXC 2004/1 - 390534769 and - Cluster of Excellence andAdvanced Imaging of Matter (AIM) EXC 2056 - 390715994. Weacknowledge computational resources provided by the Max PlanckComputing andData Facility andRWTHAachenUniversity under projectnumber rwth0811. Thisworkwas supported by theMax Planck-NewYorkCity Center for Nonequilibrium Quantum Phenomena. P.N. acknowl-edges support as a Moore Inventor Fellow through Grant No.GBMF8048 and gratefully acknowledges support from the Gordon andBetty Moore Foundation. D.N.B. is Moore Investigator in QuantumMaterials EPIQS GBMF9455. D.H. was supported by a grant from theSimons Foundation (579913).Author contributionsD.H. conducted the SNOM experiments with supervision by D.N.B. S.T.conducted and analyzed the STM and STS experimentswith supervisionby A.N.P. N.R.F. and V.S. fabricated the studied TDTG samples withsupervision by J.Hone and C.D. J.Hauck performed tight binding calcu-lationswith supervision byD.M.K. C.J.C. performed ab initio calculationsArticle https://doi.org/10.1038/s41467-022-35213-5Nature Communications |         (2022) 13:7587 7https://doi.org/10.1038/s41567-022-01556-5https://doi.org/10.1038/s41567-022-01556-5https://arxiv.org/abs/2112.07127with supervision by P.N. D.H. developed the atomic relaxation code andperformed related calculation and analysis as well as the near-fieldtomography modeling in the Supplementary Information. K.W. and T.T.grew the hBN crystals. S.T. and D.H. wrote the manuscript with con-tributions from C.J.C., J.Hauck, D.M.K., and D.N.B. D.M.K. and D.N.B.supervised the project. All authors contributed to discussions andreviewed the manuscript.Competing interestsThe authors declare no competing interests.Additional informationSupplementary information The online version containssupplementary material available athttps://doi.org/10.1038/s41467-022-35213-5.Correspondence and requests for materials should be addressed toDorri Halbertal.Reprints and permissions information is available athttp://www.nature.com/reprintsPublisher’s note Springer Nature remains neutral with regard to jur-isdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative CommonsAttribution 4.0 International License, which permits use, sharing,adaptation, distribution and reproduction in any medium or format, aslong as you give appropriate credit to the original author(s) and thesource, provide a link to the Creative Commons license, and indicate ifchanges were made. The images or other third party material in thisarticle are included in the article’s Creative Commons license, unlessindicated otherwise in a credit line to the material. If material is notincluded in the article’s Creative Commons license and your intendeduse is not permitted by statutory regulation or exceeds the permitteduse, you will need to obtain permission directly from the copyrightholder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.© The Author(s) 2022Article https://doi.org/10.1038/s41467-022-35213-5Nature Communications |         (2022) 13:7587 8https://doi.org/10.1038/s41467-022-35213-5http://www.nature.com/reprintshttp://creativecommons.org/licenses/by/4.0/http://creativecommons.org/licenses/by/4.0/ Unconventional non-local relaxation dynamics in a twisted trilayer graphene moiré superlattice Results Discussion Methods Samples preparation Exfoliation Stack preparation Near-field imaging techniques Scanning tunneling microscopy and spectroscopy Electronic structure theory calculations of generalized stacking fault energy function (GSFE) and DOS Atomic relaxation calculations Data availability References Acknowledgements Author contributions Competing interests Additional information