Spin skyrmion gaps as signatures of strong-coupling insulators in magic-angle twisted bilayer graphene

The flat electronic bands in magic-angle twisted bilayer graphene (MATBG) host a variety of correlated insulating ground states, many of which are predicted to support charged excitations with topologically non-trivial spin and/or valley skyrmion textures. However, it has remained challenging to experimentally address their ground state order and excitations, both because some of the proposed states do not couple directly to experimental probes, and because they are highly sensitive to spatial inhomogeneities in real samples. Here, using a scanning single-electron transistor, we observe thermodynamic gaps at even integer moiré filling factors at low magnetic fields. We find evidence of a field-tuned crossover from charged spin skyrmions to bare particle-like excitations, suggesting that the underlying ground state belongs to the manifold of strong-coupling insulators. From the spatial dependence of these states and the chemical potential variation within the flat bands, we infer a link between the stability of the correlated ground states and local twist angle and strain. Our work advances the microscopic understanding of the correlated insulators in MATBG and their unconventional excitations.


Introduction
Magic-angle twisted bilayer graphene (MATBG) has emerged as a remarkably rich platform to investigate correlated ground states of strongly interacting electrons [1][2][3][4][5][6].The observation of insulating electrical transport when an integer number ν of electrons/holes occupy each moire unit cell [1,3,[6][7][8][9][10] was one of the earliest indications of strong electronic correlations in MATBG and remains among their most salient manifestations.However, the nature of these correlated insulators and the novel excitations they host is still a major puzzle under active investigation.In contrast to the abundance of transport signatures of insulating ground states in samples without hBN alignment, there is little evidence of truly gapped ground states in the few measurements that directly probe the thermodynamic density of states [11][12][13][14][15][16] (the only reported example is [13]).Instead, such thermodynamic measurements consistently reveal a robust cascade of spin/valley flavor phase transitions which are present up to temperatures T ≈ 100 K [13].These transitions occur between gapless broken symmetry states that may be viewed as the parent correlated states from which the flavor-polarized insulating ground states emerge at low temperatures [17].It is therefore imperative to corroborate the existence of true gapped ground states at commensurate filling from a thermodynamic perspective.
The spin, valley, and sublattice degrees of freedom in MATBG conspire with the narrow moiré bands to form a near-degenerate manifold of competing ground states [18][19][20][21][22][23][24].In the absence of extrinsic factors such as strain or substrate alignment, theoretical considerations single out a specific "Kramers intervalley-coherent" (KIVC) insulator as the energetically favored candidate at moiré filling factors ν = 0 and ±2 [18,19,22].However, the KIVC state is part of a larger family of "strong-coupling" insulators, driven by interaction strength that is much larger than bandwidth.These insulators exhibit similar topological structure: they are formed by a coherent superposition of flavor-polarized flat Chern bands, and their energetic hierarchy is sensitive to details of the modeling.Intriguing real-space patterns as well as unconventional collective excitations are predicted for these many-body ground states [17,18].Furthermore, realistic departures from this idealized limit such as heterostrain between the two graphene sheets and local variations in twist angle can qualitatively alter the phase diagram, and have been theoretically shown to stabilize distinct but closely-competing alternatives to the strong-coupling states [23][24][25][26].Establishing the precise nature of the correlated states that emerge from this delicate competition is thus a challenging experimental question, whose answer can vary even within a single sample.
Studying the spectrum of excitations can help identify the ground state, as it can carry an imprint of the underlying broken symmetry insulator.Specifically, it has been predicted [17,[27][28][29] that the combination of strong correlations and nontrivial band topology that stabilizes a given strong-coupling insulator can also trigger a local deformation of the flavor degrees of freedom when charges are added or removed, so that the lowest-energy charged excitations are spin/pseudospin skyrmions rather than single electrons or holes.This could also have important implications for superconductivity in MATBG: it has been proposed that charge-e pseudospin skyrmions of the KIVC state could pair into charge-2e bosons that condense to drive the system into a superconducting state [17].While there is some numerical support for this scenario [28], experimental evidence has remained elusive, in part due to the challenge of directly accessing the pseudospin degrees of freedom.

Results
In this work, we use high-resolution local electronic compressibility measurements conducted with a scanning single-electron transistor (SET) to demonstrate the existence of thermodynamically gapped states at ν = ±2.These states are present at zero magnetic field B and their gap sizes are constant or increasing at low fields, but decrease linearly at moderate fields.Using Hartree-Fock calculations and symmetry analysis, we show that this nontrivial gap dependence can be rationalized in terms of a ground state with spin skyrmions as the lowest-lying charge excitations at small fields, yielding to conventional electrons and holes at larger B. The experimental observation of spin skyrmions at ν = ±2 is most naturally explained by a ground state whose pseudospin order corresponds to that of a strong-coupling insulator, among which the KIVC insulator is found to be the most energetically favorable candidate.However, since our experiment does not directly couple to the flavor degree of freedom, we cannot rule out the possibility of a related strong-coupling state which can also exhibit spin skyrmions.In the same regions where we observe a thermodynamic gap at ν = ±2, we also find evidence for a thermodynamic gap at the charge neutrality point (CNP, ν = 0) as well as a smaller total chemical potential change across the flat bands relative to areas in which the gaps are absent.These observations are consistent with a theoretical picture in which the strong-coupling state is destabilized as strain and/or twist angle variations broaden the flat bands and suppress the effect of interactions.Our work thus provides experimental evidence for spin skyrmions in MATBG, and gives microscopic insight into the role of spatial inhomogeneities in tuning the delicate balance between its competing ground states.
We first discuss the signatures of correlated insulating ground states at zero magnetic field.
Figure 1a shows a line cut of the inverse compressibility dμ/dn as a function of spatial position within the sample.We observe a sawtooth-like pattern, which has been extensively reported in MATBG [11][12][13][14][15][16]30], throughout the measured region.Superimposed on this background, the data also reveal incompressible peaks at densities corresponding to ν = ±2 in certain locations (Fig. 1b).These features represent thermodynamically gapped states, and the location at which each respective state is most robust is distinct.The incompressible peak at ν = 2 is most pronounced at θ ≈ 1.14º, and is observed over 300 nm along the measured trajectory.Its gap size as a function of spatial position is shown in Fig. 1c (see Methods).We also indicate the corresponding θ in Fig. 1c, but note that strain, which is not directly probed in our experiment, may also play a role in the observed spatial dependence.A weaker incompressible state occurs at ν = -2 at angles near θ = 1.18°.The gaps are observed only within a relatively small area and were absent in other previously measured low-disorder areas within the sample [16].This suggests that stabilizing these fragile states requires fine-tuning parameters such as twist angle and strain.
To clarify the nature of the correlated ground states, we investigate their dependence on the perpendicular magnetic field.Figure 2a-c shows dμ/dn as a function of ν and B, measured at locations 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 is not sloped in the ν -B plane, indicating it has Chern number C = 0.The high-field Hofstadter (Chern) insulators in Fig. 2a (see Supplementary Information Section 3) follow the same universal pattern reported previously, independent of the existence of the ν = ±2 insulators.The field dependence of the thermodynamic gaps Δν to the lowest available charged excitations at filling factor ν is plotted in Fig. 2d (see Methods and Supplementary Fig. S2).In the θ = 1.14º location, Δ2 decreases linearly when B ≳ 1.8 T, with a slope corresponding to an effective g-factor g* = 4.85 ± 0.03 (Fig. 2c, the error primarily comes from uncertainty in determining the lower bound of the linear fit).
However, the magnitude of the gap at low fields is smaller than predicted by extrapolating this linear trajectory.Even more strikingly, in the θ = 1.18º location, Δ±2 slightly increases with B up to about 4 -4.5 T, and then decreases linearly with a similar g* as in the 1.14º location (Fig. 2d).Qualitatively similar gap dependence is also observed at ν = ±2, suggesting similar excitations are realized (see Supplementary Information Section 2).The large g* extracted suggests the presence of an electron orbital moment that couples to the Zeeman field, even in the absence of substrate-induced C2Z symmetry breaking and without any zero-field quantum anomalous Hall signatures at odd integer fillings.This orbital moment has its origin in the Berry curvature carried by the underlying bands [31,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 understood within the framework of strongcoupling states.In the absence of strain, strong-coupling theory predicts [18][19][20][21][22][23][24] that the ground state at even integer fillings is part of a family of correlated insulators, whose topological structure at ν = ±2 is reminiscent of a quantum Hall ferromagnet consisting of two filled Landau levels with opposite C (Fig. 2e-f).In the following, we assume spin ferromagnetism appropriate for a ferromagnetic Hund's coupling for simplicity, but our conclusions are similar for the spin-valley-locked state favored by an antiferromagnetic coupling (see Supplementary Information Section 5).We find that the gap at B = 0 is controlled by spin skyrmions (Fig. 2e-g, Supplementary Fig. S3) that dominantly involve a single Chern sector (see Supplementary Information Section 4) and are the lowest-energy charged excitations at low fields.The dominant effect of an applied B field is to impose a large Zeeman penalty and reduce the skyrmion size, which also raises its Coulomb energy cost.Hence the gap initially increases at small fields, until the skyrmion excitations become energetically more costly than a single particle-hole pair.This occurs when the dashed skyrmion line in Fig. 2g crosses the lowest available particle-hole excitation, whose nature depends on quantitative details such as the energy difference δ between filled bands of the renormalized band structure.Above this critical field, the thermodynamic gap is controlled by the lowestenergy particle-hole excitation, which closes linearly because of the coupling to the orbital magnetization of the Chern bands (Fig. 2g, Supplementary Fig. S3).
To confirm the above picture, we have performed Hartree-Fock calculations on the Bistritzer-MacDonald model that demonstrate the stability 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 alternative mechanisms).Note that within our modeling, the ground state at even integer filling is the KIVC state, but we expect similar behavior from other strong-coupling candidates.We find that the lowest-energy hole excitation at ν = 2 is a spatially localized skyrmion in one Chern sector with a topologically non-trivial spin texture (Fig. 3a,b) and a non-topological pseudospin modulation (see Supplementary Information Section 4).Skyrmions do not form on electron-doping due to the larger bandwidth of the bands above the chemical potential [28,33].The initial increase of Δ2 upon including spin Zeeman coupling at low fields is consistent with the discussion above (blue squares in Fig. 3c), and the gap closes at higher fields due to the orbital coupling to the magnetic field (red line).While precise details of the microscopic modeling and fluctuations beyond mean-field will quantitatively affect the results, our calculations provide a proof-of-principle demonstration that spin skyrmions remain the lowest-energy hole excitations of the ν = 2 strong-coupling state even in realistic settings where the flat bands deviate significantly from idealized Landau level-like topological structure.
In the presence of strain, the only gapped state energetically competitive with the strong-coupling insulators is the "incommensurate Kekulé spiral" (IKS) state [23].This is a spin-singlet state at ν = ±2 whose gap strictly decreases with increasing B, and is hence difficult to reconcile 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 IKS state, is further evidence against this possibility for both the θ = 1.14º and 1.18º regions.
Under the low-strain conditions for which strong-coupling order is stabilized at ν = ±2, theory also predicts strong-coupling order at ν = 0 and hence an interaction-driven charge gap is expected at the CNP [18][19][20][22][23][24], even in the absence of substrate alignment.So far, only one transport experiment [6] has reported an activation gap at the CNP in non hBN-aligned devices.This could be due to the sensitive dependence of the low-energy physics on fine-tuning parameters like twist angle, local disorder, and strain, which complicates the interpretation of global measurements.In contrast, the scanning SET can measure locally in widely separated areas, which allows us to probe and compare multiple locations whose microscopic parameter configurations differ substantially.We find that the slope of μ(ν) in the vicinity of CNP is significantly larger in the areas where ν = ±2 gaps are observed (denoted as Area 1) compared to a far-separated region of the same device in which no gaps are present at B = 0 (denoted as Area 2 and characterized in detail in [16]), suggestive of a spontaneous gap at the CNP (Supplementary Information Section 7).In Fig. 4a, we compare μ(ν) at the θ = 1.14º location (Area 1) in Fig. 1b and the θ = 1.06º location (Area 2).The shape of μ(ν) near the CNP in Area 2 is consistent with that expected from a massless Dirac dispersion on the electron side [16], whereas in Area 1, it changes significantly more rapidly (Fig. 4a, inset), despite the significantly smaller total change in chemical potential across the flat bands, Δμtot = μ(ν=4) -μ(ν=-4).
While it is experimentally challenging to distinguish a small thermodynamic gap from a vanishing density of states at the Dirac point at B = 0, the evolution with magnetic field can provide additional evidence to help differentiate them.Figure 4b shows the extracted ν = 0 zLL gap ∆0 as a function of B. In Area 2, ∆0 follows a linear trajectory that extrapolates to 0 at zero field, whereas in Area 1, it saturates at approximately 2 meV below about 2 T, exhibiting close quantitative agreement in the 1.14º and 1.18º locations (Fig. 4b).We note that the magnitude of Δ0 in nonzero B is also much larger for Area 1, whereas the other zLL gaps show a similar magnitude (Supplementary Fig. S9) and are consistent with a phenomenological moiré valley splitting model [16].Taken together, the data strongly suggest that the CNP in Area 1 is gapped (See Supplementary Section 7 for more detail, including microwave impedance microscopy measurements).Previously reported single-particle gaps induced by hBN alignment in MATBG were 7 -10 meV [30,32], much larger than the extrapolated gap size in this work, suggesting that the CNP gap may form spontaneously as a result of interactions rather than substrate alignment.This is consistent with theoretical predictions as noted above.
The mechanism that gives rise to two distinct types of behaviors in the same device is unknown; we speculate that the ground states are highly sensitive to twist angle and strain [23][24][25][26]34], and hence postulate a minimal explanation in terms of local heterostrain variations.Whereas for small strain the strong-coupling insulator is the energetically favored ground state at ν = 0 and ±2, even modest heterostrains are known to degrade the ν = 0 gap in favor of a gapless semimetal [24][25][26], and to strongly suppress the strong-coupling gap at ν = ±2.This interpretation is also consistent with the Δμtot of the low energy bands, which is much smaller in Area 1 than Area 2 (Fig. 4c).Our calculations generally predict a trend of decreasing Δμtot with decreasing strain and with increasing twist angle (Fig. 4d), with quantitative values dependent on microscopic parameters.However, the large magnitude of the difference in Δμtot between Area 1 and Area 2 suggests that twist angle alone cannot explain the difference and strain is lower in Area 1 (Supplementary Information Section 8).Although a gap is expected to eventually open as the IKS state emerges for sufficiently large strains, mean-field theory predicts a steep suppression of the strong-coupling gap before this occurs.Our data are consistent with a scenario in which Area 2 has an intermediate level of strain, such that a gap is absent or below the resolution of our measurement.
In conclusion, our measurements provide evidence of thermodynamically gapped ground states at ν = 0, ±2.The magnetic field dependence of the ν = ±2 gapped state is consistent with a strong-coupling insulator whose low field charged excitations are spin skyrmions.Spatially resolved imaging shows that these states are stabilized when the change in chemical potential across the flat bands is small, indicating relatively low strain.In a broader context, our work establishes measurement of thermodynamic gaps as a powerful diagnostic tool for probing unconventional charged excitations.These results further motivate direct imaging of flavor ordered ground states and their topological excitations [35][36][37], especially efforts to experimentally determine whether certain types of strong-coupling ground states support pseudospin skyrmion excitations, and their relation to superconductivity.

Sample fabrication
The MATBG device is fabricated using standard dry-transfer techniques, followed by standard lithographic patterning, as detailed in [16].During the device fabrication process, we have avoided hBN substrate alignment to the graphene layers, which has been confirmed with measurements presented above and also in [16] at an independent location in the same device.

SET measurement
An SET sensor with a diameter of 50 -80 nm was brought to <50 nm above the MATBG sample surface.The resulting spatial resolution is about 100 nm.We modulate the MATBG and gate voltages, V2d,ac and Vg,ac, respectively, and measure inverse compressibility dμ/dn ∝ Ig, ac/I2d,ac, where the corresponding currents Ig, ac, I2d,ac are demodulated from the current through the SET probe at their respective modulation frequencies using standard lock-in techniques.A d.c.offset voltage V2d,dc is further applied to the sample to maintain maximum sensitivity of the SET and minimize tip-induced doping.All data presented are taken at 330 mK in an Unisoku USM1300 system with a customized microscope head.

Gate capacitance and twist angle determination
The capacitance between gate and sample (i.e., the conversion between density and applied gate voltage) was calibrated by measuring the slope of the Landau levels emanating from charge neutrality, which yielded a value consistent with geometric considerations.The twist angle is then calculated using where  = 0.246 nm is graphene's lattice constant, and ns(r) is the carrier density at full filling.

Extraction of Δ2
To extract the thermodynamic gaps, we identify either the points surrounding an incompressible peak at which dμ/dn crosses zero, i.e. the local extremum of μ(n), or the local minimum of dμ/dn if it never crosses zero.We numerically integrate dμ/dn between the filling factors corresponding to these two densities, and the resulting step in chemical potential Δμ is taken to be the gap size.However, if dμ/dn exhibits a minimum both sides of the incompressible peak, which arises if there is a large negative background from the sawtooth, then prior to integration, we shift the entire dμ/dn curve so that the less negative local minimum in dμ/dn is zero.This procedure will slightly change (< 30 μeV) the bare value of the extracted gap but has negligible effect on the extraction of g.An example demonstrating the gap extraction procedure is shown in Supplementary Fig. S2.

Hartree-Fock calculations
The properties and energetics of skyrmions of the strong coupling insulator at ν = 2 were theoretically investigated using translation symmetry-breaking Hartree-Fock calculations of the Bistritzer-MacDonald model on a periodic system of size 13×13.In agreement with prior theoretical work, the numerical ground state at ν = 2 is identified with the KIVC insulator.Upon hole-doping, the spin/pseudospin structure of the numerically obtained skyrmions is consistent with a non-linear sigma model analysis that captures the main symmetry-breaking terms of the Hamiltonian.The dependence of the thermodynamic gap ∆2 on the external Zeeman field was estimated by calculating the minimum energy Hartree-Fock solutions (including skyrmions) for total electron numbers Nν=2, Nν=2 ± 1, where Nν=2 corresponds to filling factor ν = 2.The chemical potential range Δμtot = μ(ν=4) -μ(ν=-4) was computed by evaluating the total energy for a band insulator with electron numbers Nν=-4 and Nν=4.

Microwave impedance microscopy measurements
Microwave impedance microscopy (MIM) measurements were performed in a 3He cryostat with a custom scanner incorporating Attocube nanopositioner.An etched tungsten wire was used as the MIM probe and was attached to a quartz tuning fork for topographic sensing.During measurement the tip was held approximately 20 nm above the sample's surface.The measurements reported here were carried out at 6.8 GHz.At GHz frequencies the tip and sample are strongly capacitively coupled, enabling subsurface sensing without requiring additional electrical contacts on the sample.MIM operates in the nearfield limit and the spatial resolution is dictated by the tip diameter, ~100 nm, rather than the microwave wavelength.This leads to a crossover from the skyrmion gap at small fields to the particle-hole gap at a critical field B ≈ 5 T. The Hartree-Fock calculations here are expected to overestimate the gap.wAA = 30 meV, εr = 6.
Inset, spin projection of the skyrmion along the Zeeman axis in the C = 1 sector.In-plane projection is denoted with arrows.showing that Δ rises with decreasing θ and increasing strain.