# Fileset

[s41467-023-42275-6.pdf](https://mdr.nims.go.jp/filesets/ac3dc4b8-542f-4946-ae8a-ed044b46755b/download)

## Creator

Jiachen Yu, Benjamin A. Foutty, Yves H. Kwan, Mark E. Barber, [Kenji Watanabe](https://orcid.org/0000-0003-3701-8119), [Takashi Taniguchi](https://orcid.org/0000-0002-1467-3105), Zhi-Xun Shen, Siddharth A. Parameswaran, Benjamin E. Feldman

## Rights

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

## Other metadata

[Spin skyrmion gaps as signatures of strong-coupling insulators in magic-angle twisted bilayer graphene](https://mdr.nims.go.jp/datasets/6459d56d-961f-4bb6-9756-c54bc64189a3)

## Fulltext

Spin skyrmion gaps as signatures of strong-coupling insulators in magic-angle twisted bilayer grapheneArticle https://doi.org/10.1038/s41467-023-42275-6Spin skyrmion gaps as signatures of strong-coupling insulators in magic-angle twistedbilayer grapheneJiachen Yu 1,2,8, Benjamin A. Foutty2,3,8, Yves H. Kwan4,8, Mark E. Barber 1,2,Kenji Watanabe 5, Takashi Taniguchi 6, Zhi-Xun Shen 1,2,3,7,Siddharth A. Parameswaran 4 & Benjamin E. Feldman 2,3,7The flat electronic bands in magic-angle twisted bilayer graphene (MATBG)host a variety of correlated insulating ground states, many of which are pre-dicted to support charged excitations with topologically non-trivial spin and/or valley skyrmion textures. However, it has remained challenging to experi-mentally address their ground state order and excitations, both because someof the proposed states do not couple directly to experimental probes, andbecause they are highly sensitive to spatial inhomogeneities in real samples.Here, using a scanning single-electron transistor, we observe thermodynamicgaps at even integer moiré filling factors at low magnetic fields. We find evi-dence of a field-tuned crossover from charged spin skyrmions to bare particle-like excitations, suggesting that the underlying ground state belongs to themanifold of strong-coupling insulators. From the spatial dependence of thesestates and the chemical potential variation within the flat bands, we infer a linkbetween the stability of the correlated ground states and local twist angle andstrain. Our work advances the microscopic understanding of the correlatedinsulators in MATBG and their unconventional excitations.Magic-angle twisted bilayer graphene (MATBG) has emerged as aremarkably rich platform to investigate correlated ground states ofstrongly interacting electrons1–6. The observation of insulating elec-trical transport when an integer number ν of electrons/holes occupyeach moire unit cell1,3,6–10 was one of the earliest indications of strongelectronic correlations in MATBG and remains among their most sali-ent manifestations. However, the nature of these correlated insulatorsand the novel excitations they host is still a major puzzle under activeinvestigation. In contrast to the abundance of transport signatures ofinsulating ground states in samples without hBN alignment, there islittle evidence of truly gapped ground states in the fewmeasurementsthat directly probe the thermodynamic density of states11–16 (the onlyreported example is13). Instead, such thermodynamic measurementsconsistently reveal a robust cascade of spin/valley flavor phase tran-sitions which are present up to temperatures T ≈ 100K13. These tran-sitions occur between gapless broken symmetry states that may beviewed as the parent correlated states fromwhich the flavor-polarizedinsulating ground states emerge at low temperatures17. It is thereforeimperative to corroborate the existence of true gapped ground statesat commensurate filling from a thermodynamic perspective.The spin, valley, and sublattice degrees of freedom in MATBGconspire with the narrow moiré bands to form a near-degenerateReceived: 26 May 2023Accepted: 5 October 2023Check for updates1Department of Applied Physics, Stanford University, Stanford, CA 94305, USA. 2Geballe Laboratory of Advanced Materials, Stanford, CA 94305, USA.3Department of Physics, Stanford University, Stanford, CA 94305, USA. 4Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3PU,UK. 5ResearchCenter for Electronic andOpticalMaterials, National Institute forMaterials Science, 1-1 Namiki, Tsukuba305-0044, Japan. 6ResearchCenter forMaterials Nanoarchitectonics, National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan. 7Stanford Institute for Materials and EnergySciences, SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA. 8These authors contributed equally: Jiachen Yu, Benjamin A. Foutty, Yves H.Kwan. e-mail: bef@stanford.eduNature Communications |         (2023) 14:6679 11234567890():,;1234567890():,;http://orcid.org/0000-0002-9363-5501http://orcid.org/0000-0002-9363-5501http://orcid.org/0000-0002-9363-5501http://orcid.org/0000-0002-9363-5501http://orcid.org/0000-0002-9363-5501http://orcid.org/0000-0002-3674-1607http://orcid.org/0000-0002-3674-1607http://orcid.org/0000-0002-3674-1607http://orcid.org/0000-0002-3674-1607http://orcid.org/0000-0002-3674-1607http://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-1454-0281http://orcid.org/0000-0002-1454-0281http://orcid.org/0000-0002-1454-0281http://orcid.org/0000-0002-1454-0281http://orcid.org/0000-0002-1454-0281http://orcid.org/0000-0002-5055-5528http://orcid.org/0000-0002-5055-5528http://orcid.org/0000-0002-5055-5528http://orcid.org/0000-0002-5055-5528http://orcid.org/0000-0002-5055-5528http://orcid.org/0000-0002-4962-0548http://orcid.org/0000-0002-4962-0548http://orcid.org/0000-0002-4962-0548http://orcid.org/0000-0002-4962-0548http://orcid.org/0000-0002-4962-0548http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-023-42275-6&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-023-42275-6&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-023-42275-6&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-023-42275-6&domain=pdfmailto:bef@stanford.edumanifold of competing ground states18–24. In the absence of extrinsicfactors such as strain or substrate alignment, theoretical considera-tions single out a specific “Kramers intervalley-coherent” (KIVC) insu-lator as the energetically favored candidate atmoiréfilling factors ν = 0and ± 218,19,22. However, the KIVC state is part of a larger family of“strong-coupling” insulators, driven by interaction strength that ismuch larger than bandwidth. These insulators exhibit similar topolo-gical structure: they are formed by a coherent superposition of flavor-polarized flat Chern bands, and their energetic hierarchy is sensitive todetails of the modeling. Intriguing real-space patterns as well asunconventional collective excitations are predicted for these many-body ground states17,18. Furthermore, realistic departures from thisidealized limit such as heterostrain between the two graphene sheetsand local variations in twist angle can qualitatively alter the phasediagram, and have been theoretically shown to stabilize distinct butclosely-competing alternatives to the strong-coupling states23–26.Establishing the precise nature of the correlated states that emergefrom this delicate competition is thus a challenging experimentalquestion, whose answer can vary even within a single sample.Studying the spectrumof excitations can help identify the groundstate, as it can carry an imprint of the underlying broken symmetryinsulator. Specifically, it has been predicted17,27–29 that the combinationof strong correlations and nontrivial band topology that stabilizes agiven strong-coupling insulator can also trigger a local deformation ofthe flavor degrees of freedomwhen charges are added or removed, sothat the lowest-energy charged excitations are spin/pseudospin sky-rmions rather than single electrons or holes. This could also haveimportant implications for superconductivity in MATBG: it has beenproposed that charge-e pseudospin skyrmions of the KIVC state couldpair into charge-2e bosons that condense to drive the system into asuperconducting state17. While there is some numerical support forthis scenario28, experimental evidence has remained elusive, in partdue to the challenge of directly accessing the pseudospin degrees offreedom.ResultsIn this work, we use high-resolution local electronic compressibilitymeasurements conducted with a scanning single-electron transistor(SET) to demonstrate the existence of thermodynamically gappedstates at ν = ± 2. These states are present at zero magnetic field B andtheir gap sizes are constant or increasing at low fields, but decreaselinearly at moderate fields. Using Hartree-Fock calculations and sym-metry analysis, we show that this nontrivial gap dependence can berationalized in terms of a ground state with spin skyrmions as thelowest-lying charge excitations at small fields, yielding to conventionalelectrons and holes at larger B. The experimental observation of spinskyrmions at ν = ±2 is most naturally explained by a ground statewhose pseudospin order corresponds to that of a strong-couplinginsulator, among which the KIVC insulator is found to be the mostenergetically favorable candidate. However, since our experimentdoes not directly couple to the flavor degree of freedom, we cannotrule out the possibility of a related strong-coupling state which canalso exhibit spin skyrmions. In the same regions where we observe athermodynamic gap at ν = ±2, we also find evidence for a thermo-dynamic gap at the charge neutrality point (CNP, ν =0) as well as asmaller total chemical potential change across theflatbands relative toareas in which the gaps are absent. These observations are consistentwith a theoretical picture in which the strong-coupling state is desta-bilized as strain and/or twist angle variations broaden the flat bandsand suppress the effect of interactions. Our work thus providesexperimental evidence for spin skyrmions in MATBG, and givesmicroscopic insight into the role of spatial inhomogeneities in tuningthe delicate balance between its competing ground states.We first discuss the signatures of correlated insulating groundstates at zero magnetic field. Figure 1a shows a line cut of the inversecompressibility dμ/dn as a function of spatial position within thesample. We observe a sawtooth-like pattern, which has been exten-sively reported in MATBG11–16,30, throughout the measured region.Superimposedon this background, thedata also reveal incompressiblepeaks at densities corresponding to ν = ±2 in certain locations (Fig. 1b).These features represent thermodynamically gapped states, and thelocation at which each respective state is most robust is distinct. Theincompressible peak at ν = 2 is most pronounced at θ ≈ 1.14°, and isobserved over 300nm along the measured trajectory. Its gap size as afunction of spatial position is shown in Fig. 1c (see Methods). We alsoindicate the correspondingθ in Fig. 1c, but note that strain,which is notdirectly probed in our experiment,may also play a role in the observedspatial dependence. A weaker incompressible state occurs at ν = −2 atangles near θ = 1.18°. The gaps are observed only within a relativelysmall area and were absent in other previously measured low-disorderareas within the sample16. This suggests that stabilizing these fragilestates requires fine-tuning parameters such as twist angle and strain.To clarify the nature of the correlated ground states, we investi-gate their dependence on the perpendicular magnetic field.Figure 2a–c shows dμ/dn as a function of ν and B, measured at loca-tions with θ = 1.14° and θ = 1.18° (distinct from that in Fig. 1a, b),respectively. The gapped phase at ν = ±2 survives up to B ≈ 6 T and isnot sloped in the ν - B plane, indicating it has Chern number C =0. Thehigh-field Hofstadter (Chern) insulators in Fig. 2a (see SupplementaryInformation Section 3) follow the same universal pattern reportedpreviously, independent of the existence of the ν = ±2 insulators. Thefielddependenceof the thermodynamic gapsΔν to the lowest availableba c-4 -3 -2 -1 0 1 2 3 402468101.18°1.14°2 3-0.500.511.5-6 -4 -2 0 2 4 6Vg (V)100200300400500x(nm)1.191.181.141.09-1 0 40(°)θ100 150 200 250 300 350 40000.20.40.60.811.21.17 1.14(°)θx (nm)dμ/dn (10-11 meV cm2)Δ 2 (meV)dμ/dn (10-11meV cm2 )νFig. 1 | Thermodynamically gapped states at zeromagnetic field inmagic-angletwisted bilayer graphene (MATBG). a Inverse electronic compressibility dμ/dn asa function of gate voltage Vg measured along a linear trajectory. Dashed linesdenote gate voltages corresponding to filling factors ν = ±4.b Individual line cuts ofdμ/dn as a function of ν measured at fixed locations in (a) with respective twistangles of θ = 1.14° (orange) and θ = 1.18° (blue). Inset, zoom in showing theincompressible peak near ν = 2 at θ = 1.14°. c Extracted thermodynamic gap Δ2 ofthe ν = 2 correlated insulator as a function of position and twist angle along thespatial line cut in (a).Article https://doi.org/10.1038/s41467-023-42275-6Nature Communications |         (2023) 14:6679 2charged excitations at filling factor ν is plotted in Fig. 2d (see Methodsand Supplementary Fig. S2). In the θ = 1.14° location, Δ2 decreases lin-earlywhenB ≳ 1.8 T,with a slope corresponding to an effective g-factorg* = 4.85 ± 0.03 (Fig. 2c, the error primarily comes from uncertainty indetermining the lower bound of the linear fit). However, the magni-tude of the gap at low fields is smaller than predicted by extrapolatingthis linear trajectory. Evenmore strikingly, in the θ = 1.18° location, Δ±2slightly increases with B up to about 4–4.5 T, and then decreases lin-early with a similar g* as in the 1.14° location (Fig. 2d). Qualitativelysimilar gap dependence is also observed at ν = 2, suggesting similarexcitations are realized (see Supplementary Information Section 2).The large g* extracted suggests the presence of an electron orbitalmoment that couples to the Zeeman field, even in the absence ofsubstrate-induced C2Z symmetry breaking and without any zero-fieldquantum anomalous Hall signatures at odd integer fillings. This orbitalmoment has its origin in the Berry curvature carried by the underlyingbands31,32, and contributes to an additional orbital g-factor gc: g* = gs +gc (gs = 2 for spin Zeeman coupling).The gap evolution at ν = ±2 presented above can be understoodwithin the framework of strong-coupling states. In the absence ofstrain, strong-coupling theory predicts18–24 that the ground state ateven integer fillings is part of a family of correlated insulators, whosetopological structure at ν = ±2 is reminiscent of a quantum Hall fer-romagnet consisting of two filled Landau levels with opposite C(Fig. 2e, f). In the following, we assume spin ferromagnetism appro-priate for a ferromagnetic Hund’s coupling for simplicity, but ourconclusions are similar for the spin-valley-locked state favored by anantiferromagnetic coupling (see Supplementary Information Sec-tion 5). We find that the gap at B = 0 is controlled by spin skyrmions(Fig. 2e–g, Supplementary Fig. S3) that dominantly involve a singleChern sector (see Supplementary Information Section 4) and are thelowest-energy charged excitations at low fields. The dominant effectof an applied B field is to impose a large Zeeman penalty and reducethe skyrmion size, which also raises its Coulomb energy cost. Hencethe gap initially increases at small fields, until the skyrmion excita-tions become energetically more costly than a single particle-holepair. This occurs when the dashed skyrmion line in Fig. 2g crosses thelowest available particle-hole excitation, whose nature depends onquantitative details such as the energy difference δ between filledbands of the renormalized band structure. Above this critical field,the thermodynamic gap is controlled by the lowest-energy particle-hole excitation, which closes linearly because of the coupling to theorbital magnetization of the Chern bands (Fig. 2g, Supplemen-tary Fig. S3).To confirm the above picture, we have performed Hartree-Fockcalculations on the Bistritzer-MacDonald model that demonstrate thestability of spin skyrmions in MATBG and their fingerprint on the non-monotonic gap (see Supplementary Information Sections 4 and 6,which also includes a discussion of alternativemechanisms). Note thatwithin ourmodeling, the ground state at even integer filling is the KIVCstate, but we expect similar behavior from other strong-couplingcandidates. We find that the lowest-energy hole excitation at ν = 2 is aspatially localized skyrmion in one Chern sector with a topologicallynon-trivial spin texture (Fig. 3a, b) and a non-topological pseudospin- 3 - 2 - 1 0 1 2 3 4024681011aθ = 1.14°-1 0 4dμ/dn (10-11 meV cm2)B (T)ce f gb1.8 2 2.223456-1011.5θ = 1.18° θ = 1.18°dμ/dn (10-11  meV cm2 )ν ν νΔsz = ↑sz = ↓sz = ↑sz = ↓J, λΔphΔskg* = -gc+gsg* = -gcg* = gsg* = 0skyrmionν = 2δδ < 0C = 1 C = -1 KK'C = 1C = -1KK'J, λ(T)B±2(meV)Δ0 2 4 600.511.5 1.14°1.18° ν = 2Fit1.18° ν = -2d-2024-1011.5dμ/dn (10-11  meV cm2 )-2.4Fig. 2 |Magneticfield dependence of the strong-coupling insulator. adμ/dn as afunction of ν and perpendicular magnetic field B at the θ = 1.14° location in Fig. 1a.Data is truncated beyond the limits of the colorbar. b, c dμ/dn measured at aθ = 1.18° location (distinct from that in Fig. 1a, b). d Thermodynamic gap Δ±2 atν = ±2 as a function of B at the respective spatial locations in (a–c). A linear fit of theθ = 1.14° gap (red line) for B ≳ 1.8 T. e Schematic of ν = 2 strong-coupling bandssymbolizedby twopseudo-spinors in the valley Blochspheresof eachChern sector.Dispersion and deviations from the chiral ratio manifest as couplings J, λ whichselect the strong-coupling ground state. Note that the spin structure has not beenindicated. f Schematic mean-field band diagram of a strong-coupling insulatorground state, assuming for simplicity a spin ferromagnetic configuration. The spinup/down (↑/↓) bands in the C = +1/−1 sectors are colored blue, green, red, andyellow, respectively. Δph(Δsk) denotes particle-hole (skyrmion) gap. Other types ofparticle-hole gap involving different pairs of bands are not indicated.g Energy gapsas a function of B for different possible charged excitations at ν = 2 correspondingto different spin and orbital Zeeman couplings. The gap Δ2 is set by the lowestenergy excitation. The slope of each particle-hole gap (solid lines) is given by theeffective g* from the relation: E(B) = E(B =0) + g*μBB, and depends on the directionof spin/orbitalmoments relative to the perpendicularmagnetic field.We assume gc> gs = 2 (gc, orbital g-factor of the Chern bands; gs, spin g-factor; μB, Bohr magne-ton). AB =0 energy offset δ <0 (see Supplementary Fig. S3 for the case ofδ >0, andSupplementary Information Section 4) controls the relative band alignment at lowfields. The field-dependence of the skyrmion gap (dashed lines) is governed by theZeeman suppression of spin flips.Article https://doi.org/10.1038/s41467-023-42275-6Nature Communications |         (2023) 14:6679 3modulation (see Supplementary Information Section 4). Skyrmions donot form on electron-doping due to the larger bandwidth of the bandsabove the chemical potential28,33. The initial increase of Δ2 uponincluding spin Zeeman coupling at low fields is consistent with thediscussion above (blue squares in Fig. 3c), and the gap closes at higherfields due to the orbital coupling to themagnetic field (red line). Whileprecise details of the microscopic modeling and fluctuations beyondmean-field will quantitatively affect the results, our calculations pro-vide a proof-of-principle demonstration that spin skyrmions remainthe lowest-energy hole excitations of the ν = 2 strong-coupling stateeven in realistic settings where the flat bands deviate significantly fromidealized Landau level-like topological structure.In the presence of strain, the only gapped state energeticallycompetitive with the strong-coupling insulators is the “incommensu-rate Kekulé spiral” (IKS) state23. This is a spin-singlet state at ν = ±2whose gap strictly decreases with increasing B, and is hence difficult toreconcile with the observed field dependence of the gap at θ = 1.18°(see Supplementary Information Section 4). The observed CNP gap(see below), which is absent for strains large enough to stabilize the IKSstate, is further evidence against this possibility for both the θ = 1.14°and 1.18° regions.Under the low-strain conditions for which strong-coupling orderis stabilized at ν = ±2, theory also predicts strong-coupling order atν =0 and hence an interaction-driven charge gap is expected at theCNP18–20,22–24, even in the absence of substrate alignment. So far, onlyone transport experiment6 has reported anactivation gapat theCNP innon hBN-aligned devices. This could be due to the sensitive depen-dence of the low-energy physics on fine-tuning parameters like twistangle, local disorder, and strain, which complicates the interpretationof global measurements. In contrast, the scanning SET can measurelocally in widely separated areas, which allows us to probe and com-pare multiple locations whose microscopic parameter configurationsdiffer substantially. We find that the slope of μ(ν) in the vicinity of CNPis significantly larger in the areas where ν = ±2 gaps are observed(denoted as Area 1) compared to a far-separated region of the samedevice in which no gaps are present at B = 0 (denoted as Area 2 andcharacterized in detail in ref. 16), suggestive of a spontaneous gap atthe CNP (Supplementary Information Section 7). In Fig. 4a, we com-pare μ(ν) at the θ = 1.14° location (Area 1) in Fig. 1b and the θ = 1.06°location (Area 2). The shape of μ(ν) near the CNP inArea 2 is consistentwith that expected from a massless Dirac dispersion on the electronside16, whereas in Area 1, it changes significantly more rapidly (Fig. 4a,a bcB (T)Skyrmion gapParticle-hole gap (g* = 4)0 2 4 6 8 1039.039.540.040.50-1 1δn (total)δn (C = 1) δn (C = -1)(e/AUC)-2 0(e/AUC)-2 0(e/AUC)-2 0sx/|s| (C = 1) sy/|s| (C = 1) sz/|s| (C = 1)sz/|s| (C = 1)Δ 2 (meV)0-1 1 0-1 1 0-1 1Fig. 3 | Hartree-Fock calculations of spin skyrmions of the strong-couplingstate at ν = 2. a Real-space charge density of a hole skyrmion, measured relative tothe insulator, in the periodic simulation area of 13 × 13 moiré unit cells. The parentinsulator in Hartree-Fock is the KIVC state. The charge modulation is localized inone Chern sector (here we show it in C = 1; a degenerate excitation can be found inthe opposite Chern sector by time-reversal symmetry). b Real-space spin structureof the skyrmion in the C = 1 sector. c Δ2 as a function of B, reflecting two possibleexcitations. Blue squares are numerically computed values of the skyrmion gap inthe presence of the spin Zeeman effect. Red line indicates evolution of the particle-hole gap from zero field assuming an additional effective orbital coupling g* = 4.This leads to a crossover from the skyrmion gap at small fields to the particle-holegap at a critical field B ≈ 5 T. The Hartree-Fock calculations here are expected tooverestimate the gap. wAA = 30meV, εr = 6. Inset, spin projection of the skyrmionalong the Zeeman axis in the C = 1 sector. In-plane projection is denoted witharrows.Article https://doi.org/10.1038/s41467-023-42275-6Nature Communications |         (2023) 14:6679 4inset), despite the significantly smaller total change in chemicalpotential across the flat bands, Δμtot = μ(ν = 4) - μ(ν = -4).While it is experimentally challenging to distinguish a small ther-modynamic gap from a vanishing density of states at the Dirac point atB =0, the evolution with magnetic field can provide additional evi-dence to help differentiate them. Figure 4b shows the extracted ν = 0zLL gap Δ0 as a function of B. In Area 2, Δ0 follows a linear trajectorythat extrapolates to 0 at zero field, whereas in Area 1, it saturates atapproximately 2meV below about 2 T, exhibiting close quantitativeagreement in the 1.14° and 1.18° locations (Fig. 4b). We note that themagnitude of Δ0 in nonzero B is also much larger for Area 1, whereasthe other zLL gaps show a similar magnitude (Supplementary Fig. S9)and are consistent with a phenomenological moiré valley splittingmodel16. Taken together, thedata strongly suggest that theCNP inArea1 is gapped (See Supplementary Section 7 for more detail, includingmicrowave impedance microscopy measurements). Previously repor-ted single-particle gaps induced by hBN alignment in MATBG were7–10meV30,32, much larger than the extrapolated gap size in this work,suggesting that the CNP gap may form spontaneously as a result ofinteractions rather than substrate alignment. This is consistent withtheoretical predictions as noted above.Themechanismthatgives rise to twodistinct types of behaviors inthe same device is unknown; we speculate that the ground states arehighly sensitive to twist angle and strain23–26,34, and hence postulate aminimal explanation in terms of local heterostrain variations. Whereasfor small strain the strong-coupling insulator is the energeticallyfavored ground state at ν = 0 and ±2, even modest heterostrains areknown todegrade the ν =0gap in favor of a gapless semimetal24–26, andto strongly suppress the strong-coupling gap at ν = ±2. This inter-pretation is also consistent with the Δμtot of the low energy bands,which is much smaller in Area 1 than Area 2 (Fig. 4c). Our calculationsgenerally predict a trendof decreasingΔμtot with decreasing strain andwith increasing twist angle (Fig. 4d), with quantitative values depen-dent onmicroscopic parameters. However, the largemagnitude of thedifference in Δμtot between Area 1 and Area 2 suggests that twist anglealone cannot explain the difference and strain is lower in Area 1(Supplementary Information Section 8). Although a gap is expected toeventually open as the IKS state emerges for sufficiently large strains,mean-field theory predicts a steep suppression of the strong-couplinggapbefore this occurs.Our data are consistentwith a scenario inwhichArea 2 has an intermediate level of strain, such that a gap is absent orbelow the resolution of our measurement.In conclusion, our measurements provide evidence of thermo-dynamically gapped ground states at ν =0, ±2. The magnetic fielddependence of the ν = ±2 gapped state is consistent with a strong-coupling insulator whose low field charged excitations are spin sky-rmions. Spatially resolved imaging shows that these states are stabi-lized when the change in chemical potential across the flat bands issmall, indicating relatively low strain. In a broader context, our workestablishes measurement of thermodynamic gaps as a powerful diag-nostic tool for probing unconventional charged excitations. Theseresults furthermotivate direct imaging of flavor ordered ground statesand their topological excitations35–37, especially efforts to experimen-tally determinewhether certain types of strong-coupling ground statessupport pseudospin skyrmion excitations, and their relation tosuperconductivity.MethodsSample fabricationThe MATBG device is fabricated using standard dry-transfer techni-ques, followed by standard lithographic patterning, as detailed inref. 16. During the device fabrication process, we have avoided hBNsubstrate alignment to the graphene layers, which has been confirmedwith measurements presented above and also in ref. 16 at an inde-pendent location in the same device.SET measurementAn SET sensor with a diameter of 50–80nm was brought to < 50nmabove the MATBG sample surface. The resulting spatial resolution isabout 100 nm. We modulate the MATBG and gate voltages, V2d,ac andVg,ac, respectively, and measure inverse compressibility dμ/dn ∝ Ig, ac/I2d,ac, where the corresponding currents Ig, ac, I2d,ac are demodulatedfrom the current through the SETprobe at their respectivemodulationfrequencies using standard lock-in techniques. A d.c. offset voltageV2d,dc is further applied to the sample tomaintainmaximum sensitivityof the SET and minimize tip-induced doping. All data presented aretaken at 330mK in an Unisoku USM1300 system with a customizedmicroscope head.Gate capacitance and twist angle determinationThe capacitance between gate and sample (i.e., the conversionbetween density and applied gate voltage) was calibrated by measur-ing the slope of the Landau levels emanating from charge neutrality,which yielded a value consistent with geometric considerations. Thetwist angle is then calculated using θ(r) =αffiffiffiffiffiffiffiffiffiffiffiffiffiffi3pns rð Þ8qwhere a= 0.246 nmis graphene’s lattice constant, and ns(r) is the carrier density at fullfilling.Extraction of Δ2To extract the thermodynamic gaps, we identify either the pointssurrounding an incompressible peak at which dμ/dn crosses zero,i.e. the local extremum of μ(n), or the local minimum of dμ/dn if itnever crosses zero. We numerically integrate dμ/dn between thefilling factors corresponding to these two densities, and the result-ing step in chemical potential Δμ is taken to be the gap size. How-ever, if dμ/dn exhibits a minimum both sides of the incompressiblepeak, which arises if there is a large negative background from thesawtooth, then prior to integration, we shift the entire dμ/dn curveso that the less negative local minimum in dμ/dn is zero. This pro-cedure will slightly change (< 30 μeV) the bare value of the extractedgap but has negligible effect on the extraction of g. An exampleFig. 4 | Signature of spontaneous gap at charge neutrality and spatial depen-dence. a Chemical potential μ atB =0 as a function of ν at the θ = 1.14° location andthe θ = 1.06° location reported in ref. 16. Inset, zoom in of the dispersion near ν =0.b Thermodynamic gap Δ0 at charge neutrality as a function of B in the samelocations from (a) as well as that from Fig. 2b. c Change in chemical potential Δμfrom ν = −4 to 4, extracted from four distinct spatial line cuts taken in Area 1 andArea 2. d Numerical calculation of Δμ as a function of twist angle and strain for anexample set of parameters wAA = 60meV, εr = 6, showing that Δμ rises withdecreasing θ and increasing strain.Article https://doi.org/10.1038/s41467-023-42275-6Nature Communications |         (2023) 14:6679 5demonstrating the gap extraction procedure is shown in Supple-mentary Fig. S2.Hartree-Fock calculationsThe properties and energetics of skyrmions of the strong couplinginsulator at ν = 2 were theoretically investigated using translationsymmetry-breaking Hartree-Fock calculations of the Bistritzer-MacDonald model on a periodic system of size 13 × 13. In agreementwith prior theoretical work, the numerical ground state at ν = 2 isidentified with the KIVC insulator. Upon hole-doping, the spin/pseu-dospin structure of the numerically obtained skyrmions is consistentwith a non-linear sigma model analysis that captures the mainsymmetry-breaking terms of the Hamiltonian. The dependence of thethermodynamic gap Δ2 on the external Zeeman field was estimated bycalculating the minimum energy Hartree-Fock solutions (includingskyrmions) for total electron numbers Nν=2, Nν=2 ± 1, where Nν=2 cor-responds to filling factor ν = 2. The chemical potential rangeΔμtot = μ(ν = 4) - μ(ν = -4) was computed by evaluating the total energyfor a band insulator with electron numbers Nν=-4 and Nν=4.Microwave impedance microscopy measurementsMicrowave impedance microscopy (MIM) measurements were per-formed in a 3He cryostat with a custom scanner incorporating Atto-cube nanopositioner. An etched tungsten wire was used as the MIMprobe and was attached to a quartz tuning fork for topographic sen-sing. During measurement the tip was held approximately 20 nmabove the sample’s surface. The measurements reported here werecarried out at 6.8 GHz. At GHz frequencies the tip and sample arestrongly capacitively coupled, enabling sub-surface sensing withoutrequiring additional electrical contacts on the sample.MIMoperates inthe near-field limit and the spatial resolution is dictated by the tipdiameter, ~100 nm, rather than the microwave wavelength.Data availabilityThe data that support the findings of this study are available from thecorresponding authors upon reasonable request.Code availabilityThe codes that support the findings of this study are available from thecorresponding authors upon reasonable request.References1. Cao, Y. et al. Correlated insulator behaviour at half-filling in magic-angle Graphene Superlattices. Nature 556, 80 (2018).2. Cao, Y. et al. Unconventional superconductivity in magic-anglegraphene superlattices. Nature 556, 43 (2018).3. Yankowitz, M. et al. Tuning superconductivity in twisted bilayergraphene. Science 363, 1059 (2019).4. Sharpe, A. L. et al. Emergent ferromagnetism near three-quartersfilling in twisted bilayer Graphene. Science 365, 605 (2019).5. Serlin, M. et al. Intrinsic quantized anomalous Hall Effect in a Moiréheterostructure. Science 367, 900 (2020).6. Lu, X. et al. Superconductors, orbitalmagnets andcorrelated statesin magic-angle bilayer graphene. Nature 574, 653 (2019).7. Liu, X. et al. Tuning electron correlation in magic-angle twistedbilayer graphene using Coulomb screening. Science 371,1261 (2021).8. Arora, H. S. et al. Superconductivity in metallic twisted bilayergraphene stabilized by WSe2. Nature 583, 379 (2020).9. Saito, Y., Ge, J., Watanabe, K., Taniguchi, T. & Young, A. F. Inde-pendent superconductors and correlated insulators in twistedbilayer graphene. Nat. Phys. 16, 926 (2020).10. Stepanov, P. et al. Untying the insulating and superconductingorders in magic-angle graphene. Nature 583, 375 (2020).11. Zondiner, U. et al. Cascade of phase transitions and dirac revivals inmagic-angle graphene. Nature 582, 203 (2020).12. Rozen, A. et al. Entropic evidence for a pomeranchuk effect inmagic-angle graphene. Nature 592, 214 (2021).13. Saito, Y. et al. Isospin Pomeranchuk effect in twisted bilayer gra-phene. Nature 592, 220 (2021).14. Park, J. M., Cao, Y., Watanabe, K., Taniguchi, T. & Jarillo-Herrero, P.FlavourHund’s coupling, chern gaps and charge diffusivity inMoiréGraphene. Nature 592, 43 (2021).15. Tomarken, S. L. et al. Electronic compressibility of Magic-AngleGraphene superlattices. Phys. Rev. Lett. 123, 046601 (2019).16. Yu, J. et al., Correlated Hofstadter spectrum and flavour phasediagram in magic-angle twisted bilayer graphene. Nat. Phys. 18,825–831 (2022)17. Khalaf, E., Chatterjee, S., Bultinck, N., Zaletel, M. P. & Vishwanath, A.Charged Skyrmions and topological origin of superconductivity inmagic-angle graphene. Sci. Adv. 7, eabf5299 (2021).18. Bultinck, N. et al. Ground state and hidden symmetry of magic-angle graphene at even integer filling. Phys. Rev. X 10,031034 (2020).19. Kang, J. & Vafek, O. Strong coupling phases of partially filled twis-ted bilayer graphene narrow bands. Phys. Rev. Lett. 122,246401 (2019).20. Xie, M. & MacDonald, A. H. Nature of the correlated insulatorstates in twisted bilayer graphene. Phys. Rev. Lett. 124, 097601(2020).21. Seo, K., Kotov, V. N. &Uchoa, B. FerromagneticMott state in twistedgraphene bilayers at the magic angle. Phys. Rev. Lett. 122,246402 (2019).22. Lian, B. et al. Twisted Bilayer Graphene. IV. Exact insulator groundstates and phase diagram. Phys. Rev. B 103, 205414 (2021).23. Kwan, Y. H. et al. Kekul\’e spiral order at all nonzero integer fillingsin twisted bilayer graphene. Phys. Rev. X 11, 041063 (2021).24. Wagner, G., Kwan, Y. H., Bultinck, N., Simon, S. H. & Parameswaran,S. A. Global phase diagram of the normal state of twisted bilayergraphene. Phys. Rev. Lett. 128, 156401 (2022).25. Liu, S., Khalaf, E., Lee, J. Y. & Vishwanath, A. Nematic topologicalsemimetal and insulator in magic-angle bilayer graphene at chargeneutrality. Phys. Rev. Res. 3, 013033 (2021).26. Parker, D. E., Soejima, T., Hauschild, J., Zaletel, M. P. & Bultinck, N.Strain-induced quantum phase transitions in magic-angle gra-phene. Phys. Rev. Lett. 127, 027601 (2021).27. Chatterjee, S., Bultinck, N. & Zaletel, M. P. Symmetry breaking andskyrmionic transport in twisted bilayer graphene. Phys. Rev. B 101,165141 (2020).28. Kwan, Y. H., Wagner, G., Bultinck, N., Simon, S. H. & Parameswaran,S. A. Skyrmions in Twisted Bilayer Graphene: Stability, Pairing, andCrystallization, ArXiv211206936 Cond-Mat (2021).29. Khalaf, E. & Vishwanath, A. From Electrons to Baby Skyrmions inChern Ferromagnets: A Topological Mechanism for Spin-PolaronFormation in Twisted Bilayer Graphene, ArXiv211206935 Cond-Mat (2021).30. Pierce, A. T. et al. Unconventional sequence of correlated cherninsulators in magic-angle twisted bilayer graphene. Nat. Phys. 17,11 (2021).31. Grover, S. et al. Imaging Chern Mosaic and Berry-Curvature Mag-netism in Magic-Angle Graphene, arXiv:2201.06901.32. Tschirhart, C. L. et al. Imaging Orbital Ferromagnetism in a MoiréChern Insulator. Science 372, 1323 (2021).33. Kang, J., Bernevig, B. A. & Vafek, O. Cascades between light andheavy fermions in the normal state of magic-angle twisted bilayergraphene. Phys. Rev. Lett. 127, 266402 (2021).34. Bi, Z., Yuan,N. F.Q. & Fu, L.Designingflat bandsby strain.Phys. Rev.B 100, 035448 (2019).Article https://doi.org/10.1038/s41467-023-42275-6Nature Communications |         (2023) 14:6679 635. Hong, J. P., Soejima, T. & Zaletel, M. P.Detecting Symmetry Breakingin Magic Angle Graphene Using Scanning Tunneling Microscopy,ArXiv211014674 Cond-Mat (2021).36. Liu, X. et al. Visualizing broken symmetry and topological defects ina quantum hall ferromagnet. Science 375, 321 (2022).37. Călugăru, D. et al. Spectroscopy of Twisted Bilayer Graphene Cor-related Insulators, arXiv:2110.15300.AcknowledgementsWe thank Erez Berg, Nick Bultinck, and Pablo Jarillo-Herrero forhelpful discussions. Y.H.K. and S.A.P. thank Glenn Wagner, NickBultinck, and Steven H. Simon for collaboration on related theore-tical work. This work was supported by the QSQM, an Energy FrontierResearch Center funded by the US Department of Energy (DOE),Office of Science, Basic Energy Sciences (BES), under award no. DE-SC0021238. B.E.F. acknowledges a Stanford University Terman Fel-lowship and an Alfred P. Sloan Foundation Fellowship. Y.H.K. andS.A.P. acknowledge support from the European Research Councilunder the European Union Horizon 2020 Research and InnovationProgramme, Grant Agreement No. 804213-TMCS. K.W. and T.T.acknowledge support from the JSPS KAKENHI (Grant Numbers20H00354 and 23H02052) and World Premier InternationalResearch Center Initiative (WPI), MEXT, Japan. B.A.F. acknowledges aStanford Graduate Fellowship. M.E.B. acknowledges support fromtheMarvin Chodorow Postdoctoral Fellowship of the Applied PhysicsDepartment, Stanford University. Part of this work was performed atthe Stanford Nano Shared Facilities (SNSF), supported by theNational Science Foundation under award ECCS-2026822.Author contributionsJ.Y., B.A.F. and Y.H.K. contributed equally to this work. J.Y., B.A.F., andB.E.F. designed and conducted the scanning SET experiments. B.A.F.fabricated the sample. Y.H.K. and S.A.P. conducted the theoretical cal-culations. M.E.B and Z.-X.S. designed and conducted the MIM experi-ments. K.W. and T.T. provided hBN crystals. All authors participated indiscussions and in writing of the manuscript.Competing interestsZ.-X.S. is a co-founder of PrimeNano Inc., which licensed the MIMtechnology from Stanford University for commercial instruments. Theremaining authors declare no competing interests.Additional informationSupplementary information The online version containssupplementary material available athttps://doi.org/10.1038/s41467-023-42275-6.Correspondence and requests for materials should be addressed toBenjamin E. Feldman.Peer review information Nature Communications thanks Yonglong Xie,and the other, anonymous, reviewer(s) for their contribution to the peerreview of this work. A peer review file is available.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 licence, and indicate ifchanges were made. The images or other third party material in thisarticle are included in the article’s Creative Commons licence, unlessindicated otherwise in a credit line to the material. If material is notincluded in the article’s Creative Commons licence 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 licence, visit http://creativecommons.org/licenses/by/4.0/.© The Author(s) 2023Article https://doi.org/10.1038/s41467-023-42275-6Nature Communications |         (2023) 14:6679 7https://doi.org/10.1038/s41467-023-42275-6http://www.nature.com/reprintshttp://creativecommons.org/licenses/by/4.0/http://creativecommons.org/licenses/by/4.0/ Spin skyrmion gaps as signatures of strong-coupling insulators in magic-angle twisted bilayer graphene Results Methods Sample fabrication SET measurement Gate capacitance and twist angle determination Extraction of Δ2 Hartree-Fock calculations Microwave impedance microscopy measurements Data availability Code availability References Acknowledgements Author contributions Competing interests Additional information