Charged skyrmions and topological origin of superconductivity in magic-angle graphene

See allHide authors and affiliations

Science Advances  05 May 2021:
Vol. 7, no. 19, eabf5299
DOI: 10.1126/sciadv.abf5299


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 2e 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.


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 (17), 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 (814). 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 2e topological textures. The latter can be achieved if the fundamental defects have charge 2e 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° (1518), separated by a bandgap from other bands. Recent experiments (1922) 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).


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.

Fig. 1 Defining pseudospin.

Linear combinations of the nearly flat bands of twisted bilayer graphene (left) give rise to pseudospin bands with opposite Chern numbers. Two states with the same pseudospin that lie in opposite Chern sectors have identical valley wave functions but are distinguished by A-B sublattice polarization (red and blue; third column from left). In this figure, we have omitted spin. Half-filling of this spinless model with net Chern number C = 0, at energies below the interaction scale, corresponds to quantum Hall pseudospin ferromagnets in each Chern band, described by unit vectors n±.

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 C2𝒯. Thus, the flat bands can be conveniently studied by introducing the pseudospin spinors Ψ+ = (cKA, cKB)T and Ψ = (cKB, cKA)T 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.

Table 1 Definition of new variables—Chern sector and pseudospin, from the valley and sublattice degrees of freedom.

Actions of symmetries on the internal indices are shown on the right, in both set of variables. In addition, the independent pseudospin rotations in the Chern sectors are generated by the ηP± where the projectors P±=12(1±γz) single out a specific Chern sector.

View this table:

The Hamiltonian arises from projecting the Coulomb interaction V(r) into the flat bandsH=d2k ΨkhkΨk+12d2r δn(r)V(rr)δn(r)(1)where δn(r) is the deviation of the density n(r)=γ=±Ψγ(r)Ψγ(r) from its average value. The Fourier components of 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)nq=k,γ=±Fq(k)eiγΦq(k)Ψγ,kΨγ,k+q(2)

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 nq, 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 HC=12qδnqVqδnq can be determined exactly for any purely repulsive interaction, Vq > 0 ∀q. In this case, ℋC is a sum of positive definite terms, so any state that is annihilated by each δnq ∣Ω⟩ = 0 is a ground state. This is satisfied by a manifold of 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ΩΨ±ηΨ±Ω=n±(3)

As a brief aside, we note that there is another possible set of ground states of HC 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 C2𝒯 symmetryh(k)=γxhx(k)+γyhy(k)(4)

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 C2𝒯 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 hx, 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 Jh2/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 nx, ny plane and select the K-IVC order.

Fig. 2 Tunneling-induced superexchange J.

(Top) Action of symmetries in the Chern and pseudospin basis. (Middle) Starting from the two pseudospin ferromagnets that are related by time-reversal symmetry, dispersion h acts as tunneling between C2z𝒯-related states (bottom) leading to a superexchange term Jh2/U that acts as an antiferromagnetic coupling between pseudospins in opposite Chern sectors.

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.

Fig. 3 Energetics of charged excitations.

Ratio of the elastic energy of the 2e skyrmion to the particle-hole band gap in the nearly chiral (orange) and realistic (purple) limits, with stiffness and gap values extracted from self-consistent Hartree-Fock. Below the gray dashed line, the skyrmion energy lies within the particle-hole gap. Note that w0 and w1 are the interlayer hoppings in AA- and AB-stacked regions [see (25) for details].

Fig. 4 Sigma model parameters.

Pseudospin stiffness (top), antiferromagnetic coupling J and pseudospin easy-plane anisotropy λ (bottom) as a function of the twist angle θ for different values of the gate screening distance d for dielectric constant ϵ = 9.5. All curves are computed using the analytic expressions provided in the Supplementary Materials using the band renormalization scheme of (68).

Fig. 5 Pairing of charged topological textures.

(A) Single-charge e skyrmion in one of the Chern sector. Pseudospins in n+ and n not antiferromagnetically aligned in the skyrmion core. (B) Skyrmion-antiskyrmion pair in n+n binding together due to the antiferromagnetic coupling J, which favors local pseudospin antialignment forming a charge 2e object.

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 q±(r)=14πn±·(xn±×yn±) is the topological density in each Chern sector, which integrates to unity for a single skyrmion. Let us neglect the antiferromagnetic coupling 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 ESk = 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πρpsph 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)L[n+,n]=γ=±(12AMA[nγ]·τnγ+ρps2[nγ]2)+Jn+·nμeδρ+12d2rδρ(r)V(rr)δρ(r)(5)where AM is the area of the moiré unit cell, given by 3a22θ2 (with a=3aCC). The first term is the Berry phase term, the second represents the elastic energy of the nonuniform pseudospin configurations with ρ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−, xyn+, 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) = 2e, 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 ∼JR2. 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 2e 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+, zn−, z but not n+,z2). As a result, this term reduces the binding energy of the charge 2e 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 Tc, the condensing 2e 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 Beff=2πħeAM. This leads to a Lorentz force that tends to pull them apart when they move together. Given that a skyrmion-antiskyrmion pair is bound together by 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 RL=eBeff2(R+ × R+R × R)k2(R+R)2,k=4πJ(6)and eliminating the relative coordinate R+R in favor of the center-of-mass coordinates Rs=R++R2 yields L=(eBeff)22kR·s2 from which we can read off the effective massMpair=(eBeff)2k=πħ2JAM2(7)and the transition temperature (30), related to the effective mass and stiffness throughkBTc=νπħ22AMMpair=νJAM2(8)where ν is the skyrmion filling fraction. The effective mass sets the condensation scale of the composite objects. For MATBG, JAM ∼ 1 meV, leading to the scale Tc ∼ 1 to 5 K. We will verify these estimates using a field theoretical calculation of the phase diagram with doping.

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 2e. 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 CP1 representation by writing n = zηz, where we introduced the bosonic spinon field z = (z1, z2)T, which satisfies the constraint zz = 1. (3133). 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, δρ=xayyaxπ=22πfxy (3133). 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 CP1 model. Thus, the dynamics of this model resembles that of the noncompact CP1 theory (34), which explicitly disallows monopoles.

The CP1 action takes the formS[z]=d3r{ΛgDμz2+2eμ2πcfxy+12cd2rfxy(r)πV(rr)fxy(r)π}(9)where Dμ = (∂μiaμ). Note that we rescaled the (imaginary) time as τ=1crz and introduced the cutoff scale Λ=1/AM, the coupling g=JAM2ρps, and velocity c=2Λ2JAMρps.

In the following, we compute the phase diagram of the CP1 model at finite doping and interpret the results for MATBG. We use the large N approximation by extending to N complex fields (z1, z2, …zN). 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 < gc = 4π characterized by a finite expectation value of z with aμ Higgsed and (ii) a disordered phase for g > gc = 4π where the z variables are gapped and aμ is free. The gap is given by ∣Δ∣ = Λ (1 − gc/g), and the Lagrangian for the gauge field aμ has the standard Maxwell form 1Δ(μaννaμ)2 with an additional Chern-Simons coupling to the background U(1) electromagnetic field ieπϵμνλaμμAλ. The disordered phase actually describes a paired superfluid, which can be seen by integrating out the field aμ by introducing a dual-phase variable ϕ, which can be identified with the phase of the superfluid order parameterS=d2rdτ{ρSC2(iϕ2eAi)2+χSC2(τϕ2eAτ)2}ρSC=6ΔcπN,  χSC=6ΔπcN(10)

The phase stiffness ρSC yields the temperature scale for superconductivity using the standard formula Tc=πρSC2kB (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 μc=πcΛg(1g/gc) is needed to introduce skyrmions in the ordered phase g < gc. This value reduces to the elastic skyrmion energy 4πρps for ggc, which serves as a check on the CPN 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 ggc, Δ can be expressed in terms of the doping ν leading toρSC3Λcνg2πN=3JAMνπN,χSC3Λgν2πNc=1N3ν8πAM(11)

Fig. 6 Large-N phase diagram for the doped CP1 model.

(A) Plot of the charge 2e skyrmion filling fraction ν, which is directly related to the flux fxy of the gauge field aμ, and (B) the gap for the spinon field ∆, which is directly proportional to the superfluid stiffness (cf. Eq. 10) in the superconducting phase plotted as a function of the chemical potential μ and the dimensionless coupling g. The energies μ and Δ are measured in units of the energy cutoff Λc. The black line represents the phase boundary between the insulator and the superconductor. In (C), we plot an estimate for superconducting Tc in the CPN − 1 model, computed using the relation: Tc=3cΔNkB while setting N = 2 and ρps = 1 meV as a function of doping ν for different values of J.

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 TcJν in agreement with Eq. 8. The phase diagram of the chemical potential tuned CP1 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 Δμ(k)=gμ(k)ckTηyγμck, which gives four possible pairings. Reverting to the picture of skyrmion-antiskyrmion pairing, we note that pairs reside in opposite Chern sectors, which eliminates two of the four pairing channels, leaving us with Δx, y ∝ ηyγx, y. Of these, the first is antisymmetric in internal indices leading to even pairing gx(−k) = gx(k), whereas the second is symmetric in internal indices leading to odd pairing gy(−k) = − gy(k). Further evaluating the rotation quantum number for these pairings (see the Supplementary Materials for details) (3739), one sees that the natural zero angular momentum skyrmion corresponds to Δx ∝ ηyγx = τy, while the other Δy ∝ ηyγy = τxσz option corresponds to a nonzero (odd) angular momentum.

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 Δy, which commutes with 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 CP1 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 KM and KM where we expect a gap to be induced via a Dirac mass term. To deal with the insulating and superconducting states on equal footing, we introduce the Nambu basis defined asχkT=(ψKM,kT,ψKM,k)(12)

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 formHD=kxγxρz+kyγy+(13)

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=15niΓi, where the SO(5) order parameter n̂ is defined as (n,RΔSC,IΔSC). The corresponding orders and the matrices Γi are shown in the Table 2.

Table 2 SO(5) order parameters and Dirac mass terms.

View this table:

We note that the massive Dirac Hamiltonian (Eq. 13) is invariant under the particle-hole symmetry P = γzρyK (i.e., {HD, P} = 0), where K is the complex conjugation. As a result, it belongs to symmetry class C of the Altland-Zirnbauer classification (41). The mass term ℳ parameterizes the symplectic Grassmanian manifold Sp(4n)Sp(2n) × Sp(2n) [our convention here is that Sp(2n) constitutes 2n × 2n symplectic matrices]. For our case, n = 1 and Sp(4)Sp(2) × Sp(2) is isomorphic to the four-sphere parameterized by the five-dimensional unit vector n̂. The topology of the symplectic Grassmanian π4(Sp(4n)Sp(2n) × Sp(2n))=Z is what allows for the existence of a WZW term in the action (4244).

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 theoryS=SWZW+dτd2rL[n̂](14)where L[n̂] is given byL[n̂]=ρ2(n̂)2+χ2(τn̂)2+(ggc+2χμ2)(n2Δ2)+λn32+2χμΔτΔμ2πn·(xn × yn)(15)

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.

SWZW is the well-known Wess-Zumino-Witten term, which can be written by introducing an auxillary integration variable asSWZW=3i4π01dud2rdτϵabcdenaunbτncxndyne(16)

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 n̂ is always either fully in the superconducting (12) phase or the insulating (345) phase, i.e., there is no coexistence regime. At g = gc, λ = 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 (2022). This is also supported by the cascade picture in (27, 28).

Fig. 7 MATBG close to half-filling ν = 2.

Two possible scenarios for the spin-polarized insulating state at half-filling: (A) The opposite spin band lies inside the gap of the pseudospin antiferromagnet (taken to be the K-IVC order here) or (B) the opposite spin band lies outside this gap. In both scenarios, doped electrons go into the upper K-IVC band. On the other hand, doped holes go into the opposite spin band for (A) but into the same spin K-IVC band for (B). We note here that for simplicity, we considered the case where flavor polarization corresponds to full spin polarization. Another possible flavor-polarized state would be the spin-valley locked state for which the same considerations apply.

The pairing symmetry favored by skyrmion condensation (in the spinless model) is an intervalley singlet, with same sublattice pairing, i.e., Δ̂=ψσττττyψστ. To embed the spinless model within spinful TBG at ν = 2, we need to polarize a flavor as explained above. If the polarized flavor is just spin, then we are left with pairing between equal spins. On the other hand, if the polarized flavor is more complex, for example, a spin-valley locked combination, the corresponding spin structure of the Cooper pair is also more involved. A general discussion of the spinful model is contained in the Supplementary Materials.


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 2e 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 2e 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 2e 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 Jh2/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 C2𝒯 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 2e 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 C2 symmetry such as twisted double bilayer graphene (5658) 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 C2 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 C2𝒯 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 C2, would play an antagonistic role. However, it is important to note that, in contrast to the global breaking of C2, the symmetry is not broken, on average, by disorder. For our mechanism, we only need C2 to hold on distances comparable to the moiré scale. Disorder that preserves C2, 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 (6265) 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 Jt2/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.


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πρpsPH 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 Tx, but which preserves the modified translation symmetry Txϕ which is defined as a conventional translation, followed by a global valley rotationTxϕ=eiϕ2Lxrc(r)τzc(r)Tx(17)where Lx is the length of the system in the x direction. Because the single-particle Hamiltonian commutes with Txϕ, we can label the Bloch states with the modified momentak+ϕ=(kx+ϕ2Lx,ky)(18)in one valley and modified momentakϕ=(kxϕ2Lx,ky)(19)in the other valley. In practice, this means that we simply have to shift the momentum grids in the 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 Txϕ symmetryΔϕ(k)=ψ+,k+ϕψ,kϕ(20)

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 asE0[ϕ]=E0[0]+(2ρIVC)2ALx2ϕ2+O(ϕ4)(21)where 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 material for this article is available at

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.


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.

Stay Connected to Science Advances

Navigate This Article