## Abstract

Topological solitons, a class of stable nonlinear excitations, appear in diverse domains as in the Skyrme model of nuclear forces. Here, we argue that similar excitations play an important role in a remarkable material obtained on stacking and twisting two sheets of graphene. Close to a magic twist angle, insulating behavior is observed, which gives way to superconductivity on doping. Here, we propose a unifying description of both observations. A symmetry breaking condensate leads to the ordered insulator, while topological solitons in the condensate—skyrmions—are shown to be charge 2*e* bosons. Condensation of skyrmions leads to a superconductor, whose physical properties we calculate. More generally, we show how topological textures can mitigate Coulomb repulsion and provide a previously unexplored route to superconductivity. Our mechanism not only clarifies why several other moiré materials do not show superconductivity but also points to unexplored platforms where robust superconductivity is anticipated.

## INTRODUCTION

Typically, when charge is added to an insulator, electron or hole excitations are produced, leading to charge conduction. Can solitons play the role of charge carriers? This unusual scenario can be realized through electrically charged topological textures, which have been theoretically proposed in various contexts (*1*–*7*), although physical realizations are rare. The only experimentally established instance takes place in quantum Hall ferromagnets, where topological textures of spin in the form of skyrmions acquire a charge due to the Landau level topology and are found to be the lowest energy charge excitations (*8*–*14*). On the other hand, finding situations where Cooper pairing occurs between charged topological textures, rather than between electrons, is even harder to come by. For example, in the aforementioned quantum Hall ferromagnets where charged topological textures have been experimentally established, strong breaking of time-reversal symmetry makes superconductivity highly unlikely. Obtaining robust superconductivity from topological textures typically requires simultaneously satisfying two conditions: (i) an unbroken time-reversal symmetry and (ii) the existence of stable low–energy charge 2*e* topological textures. The latter can be achieved if the fundamental defects have charge 2*e* or by pairing charge *e* defects of the same electric charge. Here, we show that all these criteria are satisfied in a simple model consisting of two time reversal–related quantum Hall (or flat Chern band) ferromagnets with purely repulsive interactions, coupled via tunneling. Such a model, we point out, captures the essential physics of magic-angle graphene, making it a promising candidate for superconductivity arising from the pairing of topological textures, a mechanism that is fundamentally different from the conventional electron-phonon mechanism.

Let us start with a brief review of magic-angle twisted bilayer graphene (MATBG). Two sheets of graphene twisted relative to one another generate a moiré pattern and, correspondingly, a reduced Brillouin zone in reciprocal space. Previous work demonstrated that the reconstruction of the electronic bands by the moiré lattice leads to minibands with extremely narrow bandwidth near charge neutrality at a magic angle ∼1° (*15*–*18*), separated by a bandgap from other bands. Recent experiments (*19*–*22*) revealed marked new physics near the same magic angle. Band gaps arising from the moiré potential are expected and observed at electron filling ν_{T} = ± 4 per moiré unit cell. In addition, insulators at other integer filling including ν_{T} = ± 2 (*19*, *21*, *22*) and in some experiments also at ν_{T} = 0 (*22*, *23*) and at certain odd integer fillings have also been observed and are attributed to the effects of electron-electron interaction. Furthermore, superconductivity has been repeatedly observed in twisted bilayer graphene, although its precise relation to the correlated insulating phase remains to be determined. While early experiments observed the superconductivity in the vicinity of the ν_{T} = − 2 correlated insulator, subsequently, a wider extent of superconductivity has been observed (*23*, *24*).

## RESULTS

### Quantum Hall ferromagnetism model of magic-angle graphene

We begin by noting an important feature of the nearly flat bands near charge neutrality. Electrons in graphene carry both a spin (*s* = ↑, ↓) and a twofold valley degeneracy (τ = *K*, *K*′) which would suggest that filling each moiré lattice site would require four electrons. In reality, it takes twice as many electrons to fill the nearly flat bands, since they consist of two connected bands as shown in Fig. 1. Conventionally, one distinguishes the two bands by their kinetic energy, but in view of the narrow bandwidth, other linear combinations may be preferable. A different choice is the sublattice basis, which separates bands on the basis of their weight on the two sites σ = *A*, *B* of the honeycomb lattice of monolayer graphene. This basis arises naturally in the ideal chiral limit of (*25*), where the sublattice polarization is complete but can also be defined at physical parameter values as explained in (*26*). Even for realistic parameters, the sublattice polarization is substantial, so a reasonable approximation is to ignore sublattice off diagonal terms in the density operator (*26*). This sublattice polarized (SP) approximation will be assumed throughout this paper.

A remarkable feature of the sublattice basis is that individual bands in this basis carry Chern number. This allows for a mapping of MATBG to a pair of time reversal–related copies of quantum Hall–like systems. The presence of flavors then relates our problem to quantum Hall ferromagnetism, as advocated in (*26*). This mapping to a generalized spin-valley ferromagnet is broadly consistent with the observation of a cascade of polarization transitions on varying electron density (*27*, *28*), and it acquires a particularly simple form in the vicinity of half-filling ν = ± 2 when the spin degree of freedom can be neglected. In this case, the low-energy physics is dominated by four flat bands, each band being labeled by a valley and sublattice index.

The valley/sublattice resolved flat bands are characterized by Chern number *C* = ± 1. Labeling sublattice (valley) by σ(τ), we have *C* = στ, which is opposite for opposite valleys due to time-reversal symmetry 𝒯 and also opposite for opposite sublattices within the same valley due to the combination of twofold rotation and time-reversal *C*_{2}𝒯. Thus, the flat bands can be conveniently studied by introducing the pseudospin spinors Ψ_{+} = (*c _{KA}*,

*c*

_{K′B})

*and Ψ*

^{T}_{−}= (

*c*,

_{KB}*c*

_{K′A})

*for the ± Chern sectors as shown in Fig. 1. A summary of the two different basis (valley/sublattice and Chern sector/pseudospin) and the relationship between them is provided in Table 1, together with the implementation of each symmetry in both basis.*

^{T}The Hamiltonian arises from projecting the Coulomb interaction *V*(**r**) into the flat bands*n*(**r**) is the deviation of the density *n*(**r**) depend on the detailed structure of the Bloch wave functions *u*_{γη;}** _{k}**(

**r**), which can be conveniently encapsulated in the “form” factor expression for the density, which can be written within the SP approximation as (

*26*)

Note that the form factor is diagonal in the Chern sectors with equal amplitude *F* and opposite phase Φ for opposite Chern sectors. It follows that *n*** _{q}**, and hence the interaction, is symmetric under independent pseudospin rotations

*U*

_{+}(2) ×

*U*

_{−}(2) in each

*C*= ± sector.

At half-filling, the ground state of the interaction *V*** _{q}** > 0 ∀

**q**. In this case, ℋ

*is a sum of positive definite terms, so any state that is annihilated by each δ*

_{C}*n*

**∣Ω⟩ = 0 is a ground state. This is satisfied by a manifold of**

_{q}*C*= 0 states, ∣Ω⟩, obtained by choosing a direction

**n**

_{±}in pseudospin space for each of the ± Chern sectors and “filling” the pseudospin bands independent of

**k**, yielding a pseudospin ferromagnet in each Chern sector. Thus

As a brief aside, we note that there is another possible set of ground states of H* _{C}* with Chern number ∣

*C*∣ = 2, obtained by filling only one of the two Chern sectors. These states, however, do not admit nontrivial topological pseudospin textures and strongly break time-reversal symmetry, making them incompatible with superconductivity and, as a result, not relevant for our discussion.

Now let us reintroduce the single-particle dispersion *h* whose form is constrained by the pseudospin and *C*_{2}𝒯 symmetry

It is important to stress here that the form of *h*(**k**) above is the most general form of the dispersion allowed by the symmetries of the chiral model, and we do not assume the flat band limit. Deviation from the chiral limit introduces an extra term proportional to η* _{z}*γ

_{0}(other terms are prohibited either by

*C*

_{2}𝒯 or particle-hole symmetry), which was shown in (

*26*) to be relatively small. This means that the main effect of deviation from the chiral limit is altering the values of

*h*

_{x, y}(

**k**) rather than introducing new terms in Eq. 4. Thus, Eq. 4 takes realistic dispersion fully into account.

The dispersion *h*(**k**) acts as tunneling between states with the same pseudospin and momentum in opposite Chern sectors. Its main effect in the limit of strong interaction is to introduce a “superexchange” term *J* ∼ *h*^{2}/*U*, where *U* is the typical interaction scale, which antiferromagnetically couples the pseudospins in opposite Chern sectors as shown in Fig. 2. This reduces the SU(2) × SU(2) symmetry of the interaction Hamiltonian, down to a single SU(2). The generation of such superexchange term is quite generic for any pair of opposite Chern number tunneling-coupled ferromagnets with SU(2) spin rotation symmetry, as shown in the Supplementary Materials, and it plays a crucial role in the skyrmion superconductivity discussed below. It also plays an important role in selecting the insulating ground state among the quantum Hall ferromagnetic states, favoring those for which the pseudospins in the two Chern sectors are antialigned, **n**_{+} = − **n**_{−} = **n**. In the language of (*26*), the resulting state corresponds to the Kramers intervalley-coherent state (K-IVC) state if *n* lies in the *XY* plane. In this phase, the conservation of valley charge is spontaneously broken (hence, intervalley coherence), which corresponds to a tripling of the unit cell on the graphene lattice scale. On the other hand, if *n* points along the *Z* direction, then the valley charge is conserved, and this order corresponds to the valley Hall state. Anisotropies neglected here that arise on going beyond the SP approximation prefer **n** ordering in the *n _{x}*,

*n*plane and select the K-IVC order.

_{y}In the following sections we show the deviation from the perfect sublattice polarization (Fig. 3), pairing of charged topological features (Fig. 4) and the energetics of charged excitations (Fig. 5), and discuss these in detail.

### Charged topological textures and effective sigma model

Because of the band topology of the two Chern sectors, pseudospin skyrmions in *n*_{±} carry electric charge (*8*, *9*) of ±*e*. The charge density δρ(**r**) = *e*(**q**_{+}(**r**) − **q**_{−}(**r**)), where *J* between the two Chern sectors for now, postponing discussion of its effect to the next section. Then, we find that the size of a charge “*e*” skyrmion in a single Chern sector is determined only by the Coulomb repulsion, which prefers to spread it out over the entire system. In this case, its energy is given only by the elastic contribution *E*_{Sk} = 4πρ_{ps}, where ρ_{ps} is the pseudospin stiffness associated with the *n*_{±} vector fields (*8*, *10*).

To determine whether skyrmions play the role of charge carriers, we must compare their energy to that of other charge carriers in the system, in particular the particle-hole excitations. This question generally depends on system details—in quantum Hall systems at ν = 1, the energy of a skyrmion pair, 8πρ_{ps}, is smaller than the particle-hole gap, Δ_{ph}, by a factor of 2. For MATBG, the ratio 8πρ_{ps}/Δ_{ph} can be computed numerically from the self-consistent Hartree-Fock calculation, leading to the results shown in Fig. 3. For model parameters close to the ideal chiral limit, the ratio is close to the quantum Hall value of 0.5, whereas for realistic parameters, it is about 0.8 to 0.9 < 1 so that, ignoring anisotropies, skyrmions are the lowest energy charge excitations. While anisotropies are absent in the SP approximation, we will leave a full calculation of skyrmion energetics in MATBG, including the most general kinds of anisotropy to future work. It is sufficient to note that they remain low-energy excitations, and we will proceed based on the premise that skyrmion excitations are important to the charge physics.

To derive the effective field theory for skyrmions, we integrate out the electrons while retaining the charge degree of freedom in the topological textures. We thus derive an effective description of MATBG by taking the pseudospin variables in the two Chern sectors *n*_{±} to be slowly varying fields leading to (see the Supplementary Materials for details)*A*_{M} is the area of the moiré unit cell, given by _{ps} ∼ 1 meV, and the third term includes the effects of antiferromagnetic coupling via the superexchange *J* ∼ 0.5 to 1 meV. The chemical potential couples to the charge deviations from the background, which is given by the skyrmion (antiskyrmion) topological charge in the + (−) Chern sector.

We note here that the deviation from the perfect sublattice polarization assumed so far leads to an additional term λ(*n*_{+, xy} · *n*_{−, xy} − *n*_{+, z} · *n*_{−, z}), with λ≈ 0.5 meV for realistic MATBG parameters (cf. Fig. 4). This term favors out-of-plane ferromagnetic coupling and in-plane antiferromagnetic coupling, thus acting as an easy-plane anisotropy, which selects the *XY* pseudospin antiferromagnetic order [this is the K-IVC order discussed in (*26*)]. The dependence of the sigma model parameter *J* and λ on the angle and interaction parameters is shown in Fig. 4.

### Superconductivity from skyrmion pairing

The antiferromagnetic coupling *J* between opposite Chern sectors leads to the binding of a skyrmion-antiskyrmion pair in the opposite Chern sectors. The net charge of this combined object is *Q* = *e*(**q**_{+} − **q**_{−}) = 2*e*, i.e., the exchange *J* has effectively resulted in an extended Cooper pair. Note the crucial interplay of antiferromagnetic coupling between opposite Chern bands; if we had had *J* < 0, then skyrmions would bind with skyrmions, leading to a charge-neutral object. Despite the long-range Coulomb interaction, such a bound state will form no matter how small *J* is, in the absence of other anisotropies. Roughly speaking, an isolated charge *e* skyrmion in a single Chern sector pays a “Zeeman” energy due to coupling to the uniform ferromagnetic order in the opposite sector via the exchange term *J*. This energy cost scales with the size *R* of the skyrmions as ∼*JR*^{2}. As in the case of quantum Hall skyrmions (*8*), the competition with the Coulomb repulsion ∼*U*/*R* will lead to a finite size for such skyrmions and yields an extra energy penalty on top of the elastic contribution. On the other hand, a pair of antiferromagnetically locked charge 2*e* skyrmions does not pay any exchange energy, which enables it to evade Coulomb repulsion by becoming very large. Hence, the extended nature of the skyrmion allows for a pairing mechanism, which evades the Coulomb repulsion while benefiting locally from the antiferromagnetic coupling.

The inclusion of the easy-plane anisotropy λ due to imperfect sublattice polarization affects the above scenario as follows. First, it leads to the deformation of skyrmions to a topologically equivalent configuration of a meron-antimeron pair to evade the penalty of out-of-plane pseudospin alignment. Second, it reduces the binding energy of antiferromagnetically locked skyrmion-antiskyrmion pair. This can be seen by noting that out-of-plane pseudospin only incurs an energy penalty for such skyrmion pairs but not for individual skyrmions (since the energy only has the term *n*_{+, z}*n*_{−, z} but not *e* skyrmion-antiskyrmion pairs and eventually leads to their unbinding when it is sufficiently large. Evidence from a recent numerical study (*29*) suggests that the binding energy remains finite for the physically relevant parameter regime, which we will assume in what follows.

Superconductivity does not follow from pairing alone. To establish a nonzero superfluid stiffness and transition temperature *T _{c}*, the condensing 2

*e*bound state must have a finite effective mass despite the flat-band dispersion of charge

*e*skyrmions. This effective mass can be generated entirely by the Coulomb repulsion through the exchange scale

*J*. This can be understood by noting that the skyrmion and antiskyrmion in opposite Chern sectors feel opposite effective magnetic fields

*J*as shown in Fig. 5B, the restoring “spring” force balances the Lorentz force leading to an effective mass. More precisely, writing the Lagrangian for a skyrmion and an antiskyrmion at

**R**

_{+}and

**R**

_{−}

*R*

_{+}−

*R*

_{−}in favor of the center-of-mass coordinates

*30*), related to the effective mass and stiffness through

*JA*

_{M}∼ 1 meV, leading to the scale

*T*∼ 1 to 5 K. We will verify these estimates using a field theoretical calculation of the phase diagram with doping.

_{c}### Field theory description of the skyrmion superconductor at finite chemical potential

The antiferromagnetic coupling implies that the **n**_{+} and **n**_{−} pseudospins are antiferromagnetically locked in the ground state. In addition, the lowest energy charge excitation are also bound states of skyrmions in **n**_{+} and of antiskyrmions in **n**_{−} where **n**_{+} = − **n**_{−} = **n**, which can be understood as *n*-skyrmions, which carry charge 2*e*. Thus, we can integrate out the massive ferromagnetic fluctuations, **n**_{+} + **n**_{−}, whose mass is proportional to *J*. The resulting field theory is then written solely in terms of the SO(3) vector *n*. Furthermore, we can rewrite this in the well-known CP^{1} representation by writing **n** = *z*^{†}**η***z*, where we introduced the bosonic spinon field *z* = (*z*_{1}, *z*_{2})* ^{T}*, which satisfies the constraint

*z*

^{†}

*z*= 1. (

*31*–

*33*). The overall phase of

*z*is redundant, which leads to the gauge field

*a*

_{μ}= −

*iz*

^{†}∂

_{μ}

*z*. The topological density is given by the flux of

*a*

_{μ}and is tied to the charge density. Thus,

*31*–

*33*). Note that since the flux of

*a*is tied to a conserved charge, monopole fluctuations, which change the flux in units of 2π, are disallowed on symmetry grounds, unlike in the usual CP

^{1}model. Thus, the dynamics of this model resembles that of the noncompact CP

^{1}theory (

*34*), which explicitly disallows monopoles.

The CP^{1} action takes the form*D*_{μ} = (∂_{μ} − *ia*_{μ}). Note that we rescaled the (imaginary) time as

In the following, we compute the phase diagram of the CP^{1} model at finite doping and interpret the results for MATBG. We use the large *N* approximation by extending to *N* complex fields (*z*_{1}, *z*_{2}, …*z _{N}*). We start by reviewing the solution in the absence of doping and Coulomb interaction (λ = μ =

*V*= 0), which was discussed in the pioneering works of (

*31*,

*32*,

*35*). In this limit, the model has two phases: (i) an ordered phase for

*g*<

*g*= 4π characterized by a finite expectation value of

_{c}*z*with

*a*

_{μ}Higgsed and (ii) a disordered phase for

*g*>

*g*= 4π where the

_{c}*z*variables are gapped and

*a*

_{μ}is free. The gap is given by ∣Δ∣ = Λ (1 −

*g*/

_{c}*g*), and the Lagrangian for the gauge field

*a*

_{μ}has the standard Maxwell form

*U*(1) electromagnetic field

*a*

_{μ}by introducing a dual-phase variable ϕ, which can be identified with the phase of the superfluid order parameter

The phase stiffness ρ* _{SC}* yields the temperature scale for superconductivity using the standard formula

*30*). Let us now turn to producing the superconductor by varying the chemical potential, rather than the coupling

*g*.

For nonzero chemical potential μ, the phase diagram can be obtained numerically by solving the mean-field self-consistency equation (see the Supplementary Materials) as a function of *g* and μ as shown in Fig. 6. Note that a finite value *g* < *g _{c}*. This value reduces to the elastic skyrmion energy 4πρ

_{ps}for

*g*≪

*g*, which serves as a check on the CP

_{c}*theory. For larger μ, the skyrmion density ν is obtained numerically. For any finite skyrmion density ν, the effective field theory takes the form (Eq. 10) describing a superconductor whose phase stiffness can be computed numerically as shown in Fig. 6. In the dilute limit ν ≪ 1 and for*

^{N}*g*≪

*g*, Δ can be expressed in terms of the doping ν leading to

_{c}The superfluid stiffness (charge compressibility) of the superconductor is inversely proportional to the pseudospin compressibility (psuedospin stiffness) of the zero doping insulator with the proportionality constant being, up to numerical prefactors, just the filling ν. This leads to *T _{c}* ∼

*J*ν in agreement with Eq. 8. The phase diagram of the chemical potential tuned CP

^{1}theory with realistic parameters and screened Coulomb interaction should also be accessible in future studies using the Monte Carlo technique (

*36*).

To fix the pairing symmetry of the superconductor, consider the symmetries of the operator Δ that creates flux, ∇ × *a* = 2π, which corresponds to the Cooper pair creation operator. Since flux is left invariant by pseudospin rotation, Δ is expected to be a pseudospin singlet. Hence, the Cooper pair operator must be one of _{x, y} ∝ η* _{y}*γ

_{x, y}. Of these, the first is antisymmetric in internal indices leading to even pairing

*g*(−

_{x}**k**) =

*g*(

_{x}**k**), whereas the second is symmetric in internal indices leading to odd pairing

*g*(−

_{y}**k**) = −

*g*(

_{y}**k**). Further evaluating the rotation quantum number for these pairings (see the Supplementary Materials for details) (

*37*–

*39*), one sees that the natural zero angular momentum skyrmion corresponds to Δ

*∝ η*

_{x}*γ*

_{y}*= τ*

_{x}*, while the other Δ*

_{y}*∝ η*

_{y}*γ*

_{y}*= τ*

_{y}*σ*

_{x}*option corresponds to a nonzero (odd) angular momentum.*

_{z}Separately, it is worth noting that the pairing channels Δ_{x, y} are also the ones that correspond to the maximal gap in the presence of the insulating pseudospin antiferromagnetic background. This is readily seen by checking that the corresponding matrices anticommute in the Nambu basis (see the Supplementary Materials); thus, the insulating and superconducting gaps add in quadrature when they are simultaneously present. Furthermore, taking into account the kinetic part of the Hamiltonian *h*, in Eq. 4, it is found to anticommute with Δ* _{x}*, leading to a bigger gap compared to Δ

*, which commutes with*

_{y}*h*.

### SO(5) sigma model and WZW term

In the previous sections, we have mainly focused on the transition into the superconductor starting deep inside the K-IVC state by tuning the chemical potential, treated within the CP_{1} theory. However, following (*4*, *7*, *40*), it is instructive to derive an equivalent effective field theory, which deals with the insulator and superconductor on equal footing. Such a field theoretic description is known to include a topological WZW term. In the following, we will outline this derivation in the context of magic-angle graphene including the effect of the chemical potential.

To derive universal aspects of the field theory, such as the presence of a topological term, it is sufficient to adopt a convenient starting point where we take a dispersion with Dirac points for the nearly flat bands with a spontaneously induced single-particle gap (Dirac mass) that is much smaller than the dispersion. The dispersion can then be linearized close to the moiré Dirac points *K*_{M} and

Let us now introduce the Pauli matrices ρ_{x, y, z}, which act within the two-dimensional Nambu space in χ** _{k}**. The Hamiltonian now takes the form

The mass ℳ is a matrix in γ, η, and ρ spaces containing both the K-IVC and the superconducting order parameters, as well as the valley Hall order, which is a part of the antiferromagnetic manifold. All these correspond to anticommuting mass terms for the Dirac equation and hence can be written as * _{i}* are shown in the Table 2.

We note that the massive Dirac Hamiltonian (Eq. 13) is invariant under the particle-hole symmetry P = γ* _{z}*ρ

*K (i.e., {H*

_{y}*, P} = 0), where K is the complex conjugation. As a result, it belongs to symmetry class C of the Altland-Zirnbauer classification (*

_{D}*41*). The mass term ℳ parameterizes the symplectic Grassmanian manifold

*Sp*(

*2n*) constitutes 2

*n*× 2

*n*symplectic matrices]. For our case,

*n*= 1 and

*42*–

*44*).

Following the standard procedure by integrating out the fermions and performing the gradient expansion (see the Supplementary Materials for details), we derive the following effective theory

We see that the chemical potential enters the action in three different places. First, it couples to the topological density of the vector *n* representing the insulator. Note that since *n* now does not have a fixed length, this term is not quantized. Second, the chemical potential couples linearly to the first derivative of the superconducting gap and quadratically to its magnitude. The latter coupling affects the effective potential, favoring superconductivity at finite doping. We note that both couplings were considered before in the context of SO(5) theories for cuprates (*45*, *46*), although the topological term, which we turn to next, was absent in that context.

*S*_{WZW} is the well-known Wess-Zumino-Witten term, which can be written by introducing an auxillary integration variable as

At μ = 0, this theory has the same form as the theory describing the transition between an antiferromagnet and a valence bond solid (*37*, *40*, *47*, *48*). In this case, the vector *g* = *g _{c}*, λ = 0, the model has been conjectured to have an emergent SO(5) symmetry (

*49*,

*50*), reduced to an SO(4) symmetry in the easy-plane limit. The existence of a WZW term implies that one cannot simultaneously disorder both the superconductor and the K-IVC/Valley Hall order without inducing either a gapless critical point (deconfined criticality) or topological order.

### Role of spin quantum number

The spinless model discussed here is directly relevant to the vicinity of the half-filling insulator [at fillings ν = ± (2 + ϵ)], where charge neutrality corresponds to ν = 0. Achieving an insulator at half-filling requires polarizing a flavor (e.g., spin) and opening a gap at the Dirac point, presumably through developing a pseudospin quantum Hall ferromagnet (with **n**_{+} = − **n**_{−}). Spin polarization means that for each valley, *K* and *K*′, only one spin species is filled, which may correspond to a simple spin ferromagnet or a spin-valley locked state with spin anti-aligned in opposite valleys. The main assumption leading to the reduction to a spinless problem is that spin polarization takes place at a higher-energy scale than the scale relevant for superconductivity. This is strongly suggested by the cascade picture of (*27*, *28*), which identified a relatively high temperature scale where flavor polarization takes place, before the development of insulating or superconducting behavior at lower temperatures. As a result, one spin flavor is completely frozen at the scales relevant for superconductivity, allowing us to focus on the spinless problem. We note that most elements of our discussion remain valid even without the assumption of a frozen spin flavor—pseudospin skyrmions still carry charge and pair due to antiferromagnetic coupling—but the ground-state manifold has a more complicated structure and allows for other competing charged excitations.

It is worth noting that while the spinless model is particle-hole symmetric, we do not expect particle-hole symmetry in the real system in the vicinity of half-filling. To see this, note that there are two scenarios for the ν = ± insulator shown schematically in Fig. 7 where the opposite spin band lies outside or inside the Dirac gap. In the first scenario (i), electron (hole) doping at ν = 2(−2) resembles doping a spinless version of the charge neutrality state, whereas hole (electron) doping resembles doping a spinless version of the full (empty) band structure; whereas in the second scenario (ii), there is no difference between particle and hole doping at half-filling and both reduce to a spinless version of the doped state near charge neutrality. The first scenario is strongly supported by the quantum oscillation measurements where the Landau fans are only observed at ν = ± (2 + ϵ) but not at ν = ± (2 − ϵ) with a degeneracy that is half that observed at charge neutrality (*20*–*22*). This is also supported by the cascade picture in (*27*, *28*).

The pairing symmetry favored by skyrmion condensation (in the spinless model) is an intervalley singlet, with same sublattice pairing, i.e.,

## DISCUSSION

Let us compare and contrast the present work with previous theoretical discussions of skyrmion superconductivity. Starting with Dirac fermions, (*3*, *4*) pointed out that charge 2*e* skyrmions can arise when interactions spontaneously create topological insulators from symmetry breaking. Further, if these skyrmion textures are the lowest energy charge excitations, (*4*, *51*, *52*) argued that doping could lead to superconductivity. Here, we have examined a microscopic model with predominantly repulsive interactions and shown that it exhibits a related topological ground state with low–energy charge 2*e* skyrmions. Note that restricting to repulsion is necessary to establishing an electronic mechanism of superconductivity. To the best of our knowledge, this collection of properties has not been established before (*53*) in a microscopic model with repulsive interactions. A distinction from the previous Dirac models is that the 2*e* skyrmion here is actually a composite of a skyrmion-antiskyrmion pair. Furthermore, our model is motivated by magic-angle graphene, which seems to naturally incorporate the requisite physics. We also have calculated the energetics of skyrmions and shown that in a range of parameters, they are the lowest energy charge excitations. On doping skyrmions, we obtain a superconductor, whose transition temperature we calculate. In our theory, an essential role is played by *J*, which induces pairing between opposite skyrmion textures and arises from electron tunneling between Chern sectors. In particular, we provide a scenario where pairing can be achieved in the absence of any extrinsic attractive pairing mechanism, e.g., phonons. In addition, our mechanism explains how the large Coulomb repulsion can be evaded without screening or retardation, allowing for the smaller pseudospin antiferromagnetic coupling *J* ∼ *h*^{2}/*U* ∼ 1 meV to overcome the larger Coulomb repulsion as summarized in Fig. 5.

The role of the coupling *J* as the source of pairing can account for the absence of superconductivity in MATBG samples, which are aligned with the hexagonal boron nitride (hBN) substrate. The latter generates a sublattice potential, which breaks *C*_{2}𝒯 symmetry by inducing an energy gap of about 15 to 20 meV between the bands polarized on sublattices A and B (*13*, *54*, *55*). This acts as an opposite Zeeman field in the opposite Chern sectors (a Zeeman field for the vector *n*), which (i) shrinks the charge 2*e* skyrmions and makes them energetically less favorable and (ii) introduces an energy detuning between previously resonant bands, which suppresses the tunneling term responsible for the coupling *J*. This mechanism is also absent in other moiré systems lacking *C*_{2} symmetry such as twisted double bilayer graphene (*56*–*58*) and ABC graphene on hBN (*59*), which can also be viewed as consisting of two quantum Hall systems with opposite Chern numbers. In addition to the suppression of *J* due to detuning from the broken *C*_{2} symmetry, the direct tunneling *h* is also strongly suppressed in these systems since the two opposite Chern sectors reside in opposite valleys and are thus very weakly coupled due to valley-charge conservation. Apart from pristine MATBG samples, the other platform that retains the relevant *C*_{2}𝒯 symmetry is twisted trilayer (and multilayer) graphene with alternating twist as described in (*60*). In addition to having the same symmetries as MATBG, it was shown (*60*) that the wave functions of these systems can be mapped exactly to those of MATBG, making them very attractive candidates for pseudospin antiferromagnetic order and skyrmion-mediated superconductivity and are promising platforms for future theoretical and experimental studies.

A natural question that arises is whether disorder, which is expected to locally break *C*_{2}, would play an antagonistic role. However, it is important to note that, in contrast to the global breaking of *C*_{2}, the symmetry is not broken, on average, by disorder. For our mechanism, we only need *C*_{2} to hold on distances comparable to the moiré scale. Disorder that preserves *C*_{2}, on average, and is sufficiently smooth on the moiré scale is therefore not expected to substantially affect our conclusions.

It is worth emphasizing that the superconducting mechanism we discussed here is applicable to a much wider class of systems than just magic-angle graphene. The key requirements are a pair of quantum Hall (or Chern) ferromagnets, related by time-reversal symmetry, and coupled to one another by tunneling. This leads, via a superexchange process, to an antiferromagnetic coupling *J*, which is maximized when tunneling connects states at the same energy.

Let us mention a few other promising systems that incorporate all these ingredients. First, we have already briefly mentioned the multilayer graphene platform as described in (*60*), where each successive layer is twisted by an alternating angle (i.e., θ, − θ, θ…), and where alternating layers, which are at the same angle, are in registry. These systems have shifted magic angles but otherwise retain all essential symmetries and are promising platforms for realizing skyrmion superconductivity. For the simplest alternating angle trilayer case, an additional ingredient is a coexisting Dirac node, which intersects the flat bands, but which is not expected to significantly alter the physics due to its rapid dispersion (*60*). A second route involves strain-induced Landau levels in graphene, which have been observed (*61*), with opposite Chern number in opposite valleys. Here, spin plays the role of our pseudospin, although translation symmetry blocks tunneling between opposite valleys. Activating tunneling with an appropriate translation symmetry breaking order (e.g., a substrate with charge density wave order) can help realize the skyrmion superconductor on doping. Other potential realizations include any spinful valley-Chern insulator, e.g., twisted double bilayer graphene, where again the crucial additional ingredient absent in current platforms is the presence of a source of translation symmetry breaking that opens up intervalley tunneling. More generally, two-dimensional materials with electronic bands harboring the same fragile topology as magic angle graphene (*62*–*65*) will be interesting future candidates if they can be also brought into the correlated regime. Alternately, consider a pair of Haldane models with opposite Chern numbers, with spin degeneracy. Here, spins play the role of our pseudospin. We note that such a model, dubbed as the “shift insulator,” has been recently studied in the context of crystalline topology (*38*). Further, adding repulsive interactions and tunneling between the opposite Chern sectors should realize the physics studied here when placed near quarter filling. Pursuing materials realization of these models, such as bilayers of anomalous quantum Hall ferromagnets, will be an interesting direction for future work.

Last, let us briefly discuss our proposal in light of recent experiments on MATBG superconductivity. First, although we have chosen to discuss our scenario in terms of doping a parent pseudospin-ordered insulator, the presence of the insulating gap is not essential. While an insulating state serves as proxy for order, it is by no means necessary, since an ordered metallic phase can also obtain at integer filling, due to low-energy electron-hole pockets. The essential features of the skyrmion pairing mechanism, however, remain intact, although a theoretical description must now include fermions in addition to the bosonic **n** degrees of freedom. Hence, our results are consistent with recent experimental reports of superconductivity in the absence of correlated insulating behavior (*23*, *24*)—an interesting question there is whether pseudospin order is present in the metallic state near integer filling. We note that in the calculation of The large N phase diagram, we find a superconducting phase with doping in the absence of long-range pseudospin order. Thus strictly, the only real requirement is short-ranged pseudospin order in which skyrmions can be defined. Second, an important energy scale in our theory is *J* ∼ *t*^{2}/*U*, which controls the strength of superconductivity. This may provide a conceptual framework to help design better superconductors. For example, it is helpful to remember that in the interacting problem, not just *U* but also *t* increases with increasing the Coulomb scale, since the flat band dispersion is mostly generated by the interaction as discussed in (*26*).

The discovery of superconductivity in MATBG came as a bolt from the blue, and it would be unexpected if an explanation only invoked conventional ingredients. Instead, we believe that this observation calls for an unusual mechanism like the skyrmion pairing one described here, which relies on characteristics unique to MATBG. More generally, we have identified a new mechanism for superconductivity from repulsion, and it will be interesting to explore in the future a variety of other settings that realize the essential ingredients of Chern ferromagnets that are tunnel coupled to their time-reversed conjugates.

## MATERIALS AND METHODS

In this section, we describe the procedure used to calculate the K-IVC stiffness ρ_{IVC} numerically. This pseudospin stiffness, referred to as ρ_{ps} previously, is used to estimate the energy ratio 8πρ_{ps}/Δ_{PH} in Fig. 3.

We start by threading valley flux ϕ in, say, the *x* direction, and we choose a gauge that manifestly breaks translation in the *x*-direction *T _{x}*, but which preserves the modified translation symmetry

*L*is the length of the system in the

_{x}*x*direction. Because the single-particle Hamiltonian commutes with

*x*direction in opposite ways in the two valleys.

On the shifted momentum grids, we numerically solve the Hartree-Fock self-consistency conditions allowing for an intervalley coherence order parameter of the most general form consistent with the modified

In this expression, sublattice and spin indices are implicit.

Following this procedure, we obtain the ground-state energy of the K-IVC state as a function of ϕ. On general grounds [see, e.g., (*66*)], one knows that the ground-state energy depends on the flux as*A* is the total area of the system. Note that we define ρ_{IVC} to be half of the actual stiffness to get rid of a spin degeneracy factor. We obtain the stiffness by fitting to this quadratic function at small ϕ < 0.5. The simulations used to obtain the results shown in the main text were done on a 18 × 18 momentum grid, keeping only the two flat bands per spin and valley. For the Coulomb interaction, a dual-gate screened potential was used with a gate distance of 20 nm and a dielectric constant ϵ = 9.5. Note that after this work was released, (*29*, *67*) appeared, on related topics.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/19/eabf5299/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank S. Liu for a related previous collaboration and S. Sachdev for the helpful feedback. S.C. thanks V. Lohani for the helpful discussions.

**Funding:**A.V. was supported by a Simons Investigator award and by the Simons Collaboration on Ultra-Quantum Matter, which is a grant from the Simons Foundation (651440 to A.V.). E.K. was supported by a Simons Investigator Fellowship, by NSF-DMR 1411343, and by the German National Academy of Sciences Leopoldina through grant LPDS 2018-02 Leopoldina fellowship. S.C. acknowledges support from the ERC synergy grant UQUAM via E. Altman. M.P.Z. was supported by the Director, Office of Science, Office of Basic Energy Sciences, Materials Sciences, and Engineering Division of the U.S. Department of Energy under contract no. DE-AC02-05-CH11231 (van der Waals heterostructures program; KCWF16).

**Author contributions:**All authors contributed to generating ideas and carrying out calculations and took part in the writing of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Supplementary materials includes Supplementary Text detailing derivations and calculation methods for results in the main text and additional references. Additional data related to this paper may be requested from the authors.

- Copyright © 2021 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).