Theory of correlated insulating behaviour and spin-triplet superconductivity in twisted double bilayer graphene.

Two graphene monolayers twisted by a small magic angle exhibit nearly flat bands, leading to correlated electronic states. Here we study a related but different system with reduced symmetry - twisted double bilayer graphene (TDBG), consisting of two Bernal stacked bilayer graphenes, twisted with respect to one another. Unlike the monolayer case, we show that isolated flat bands only appear on application of a vertical displacement field. We construct a phase diagram as a function of twist angle and displacement field, incorporating interactions via a Hartree-Fock approximation. At half-filling, ferromagnetic insulators are stabilized with valley Chern number [Formula: see text]. Upon doping, ferromagnetic fluctuations are argued to lead to spin-triplet superconductivity from pairing between opposite valleys. We highlight a novel orbital effect arising from in-plane fields plays an important role in interpreting experiments. Combined with recent experimental findings, our results establish TDBG as a tunable platform to realize rare phases in conventional solids.

Introduction

The recent discovery of correlated insulating states and superconductivity in twisted bilayer graphene (TBG)[1-4] has opened a new window to exploring strong correlation effects in systems whose doping can be easily tuned, enabling the exploration of a rich range of interaction-driven phenomena. Although the underlying reason for the correlated physics is understood to arise from a relatively narrow electronic bandwidth induced by the long wavelength Moiré pattern[5,6], several details, including the symmetry breaking within the insulating phase and the nature and mechanism of pairing in the neighboring superconductor, remain under debate[7-19]. One of the difficulties in addressing these questions arises from the complexity of the theoretical treatment of TBG, which involves at least a pair of narrow bands per spin per valley with a symmetry-protected band touching, leading to eight bands in total. On top of that, the limited tunability of the band structure makes it experimentally difficult to explore the dependence of different phases on microscopic parameters.

Motivated by recent experimental report[20], we study a related system—twisted double-bilayer graphene (TDBG)—which consists of a pair of bilayer-graphene sheets, twisted with respect to one another with AB–AB-stacking structure. Due to the absence of rotation symmetry, TDBG has a lower symmetry compared with TBG, which simplifies the problem by removing the band touching at the Dirac points, leading to a low energy effective description involving one rather than two narrow bands per spin and valley. Moreover, the band separation can be controlled by applying a vertical displacement field enabling the exploration of different regimes of band isolation and bandwidth within the same device.

We identify three main ingredients necessary to explain the emergence of insulating and superconducting behavior in TDBG. First, we perform an accurate calculation of the single-particle band structure to identify ranges of displacement field and twist angle for which a single band is isolated and relatively flat. We show that lattice relaxation, known to be important in TBG[21,22], as well as several other effects such as trigonal warping, which are absent in TBG, significantly influence the band structure in TDBG, in excellent agreement with experiments. Moreover, we identify a hitherto-neglected in-plane orbital effect which is used to explain the experimentally observed deviation of the in-plane factor from 2[20], as well as the effect of in-plane field on superconducting .

Second, we address the key question of the nature of the interaction-driven insulating state. The similarity between the phase diagram of TBG to that of cuprates was invoked to argue that Mott physics is the underlying mechanism responsible for the correlated insulator[1,7,12]. On the other hand, a different route to correlated insulators is observed in graphene quantum-Hall systems, for instance, when the spin and valley degeneracy of the Landau levels are spontaneously broken by interactions[23]. This usually leads to ferromagnetic insulators, which are otherwise rare in correlated solids where antiferromagnetic order is the norm. For similar reasons, in the TDBG with nonzero valley Chern number, ferromagnetism may be preferred[24] at integer fillings. The situation here is reminiscent of strained graphene, where a suitably chosen strain profile leads to Landau levels arising from the opposite strain magnetic fields applied on the two valleys[25]. At partial fillings that are integers, ferromagnetic ground states were obtained with repulsive interactions[26], and we show that a similar scenario is likely to occur here in TDBG. Indeed a related ground state with spontaneous quantum-Hall response, although metallic, was observed in the twisted monolayer-monolayer graphene (TBG) with -breaking substrate potentials[13,19,24,27-29].

Third, we investigate the nature of the superconducting phase by highlighting that the valley degree of freedom, which behaves as a pseudospin, allows for exotic pairing possibilities which are relatively rare in other materials. In particular, we show that spin triplet with valley-singlet pairing, which is momentum-independent within each valley, is favored. We investigate the consequences of such scenario and show it can be used to explain the measured dependence of on in-plane field[20].

Results

Single-particle physics

We consider a system consisting of two AB-stacked graphene bilayers twisted relative to AB–AB stacking by a small angle , illustrated in Fig. 1. For a detailed discussion on the Hamiltonian and model parameters, see the Methods section. The bottom layer of the top BLG and the top layer of the bottom BLG are coupled via Moiré hopping between and sites, parametrized by [21,22]. In the original Bistritzer–Macdonald model, and are taken to be equal[30]. However, in a realistic twisted model, the ratio is smaller than one due to the lattice relaxation which expands (shrinks) AB (AA) regions. In TBG, is taken to be ~0.75 for the first magic angle[21,22]. Here, we similarly include lattice relaxations by taking to be <1. This is crucial for the existence of a gap between first and second conduction (valence) bands in TDBG which is necessary to explain the band insulator at filling. In this work, we take corresponding to . For different values of , we obtained qualitatively similar features (Methods).

Fig. 1

Twisted double BLG model (AB–AB stacking) with the gating voltage across the system. Throughout this work, we assume the voltage drop across the layers is uniform, .

Unlike TBG, a realistic description of TDBG does not exhibit magic-angle physics whose origin is the vicinity to a chiral symmetric model with perfectly flat bands at specific angles[31,32]. In the quadratic approximation of the bilayer-graphene dispersion, the first conduction and valence bands in TDBG become almost perfectly flat at the angle [24]. However, once trigonal warping () and particle–hole asymmetry () terms are included, the flat bands acquire a significant dispersion and become overlapped with each other (Fig. 2a, b). Theses bands can only be separated by applying a strong enough gate voltage between top and bottom layers (Fig. 2c). Using numerical simulations, we identify the parameter space of twist angle and applied voltage where the first conduction band is isolated (Fig. 3a). On the other hand, we find that there is barely any regime where the first valence band is isolated (Fig. 3c). Such a particle–hole asymmetry in the band structure is originated from and terms. The results are consistent with the experimental findings[20], showing that the system at charge neutrality remains metallic unless a rather large vertical electric field is applied. Furthermore, a correlated insulating phase is only observed on electron-doping side, consistent with the theoretically expected particle–hole asymmetry. Note that the bandwidth is not as flat as that of magic-angle TBG. However, the bandwidth is still small compared with the interaction scale which implies that strongly correlated physics can still arise. Indeed, there is some debate regarding the bandwidth of magic-angle TBG itself, with reported bandwidths ranging from 10 to 40 meV[33].

Fig. 2

Moiré band structures of TDBG. a, b At . Solid (dotted) line represents the band originated from () valley. Red, blue and black represent conduction, valence, and the other bands, respectively. a The band structure for the idealized model with only and being nonzero. The flat band is observed with the bandwidth  meV. b The band structure for the realistic model with overlapping bands. The “magic angle” does not exist in this case. c Moiré band structure at . The first conduction band (red) is isolated and relatively flat.

Fig. 3

Summary of single-particle calculations of TDBG. a Isolation region for the first conduction band (colored) with the bandwith indicated by the color. We observe two seperate isolation regions for smaller or larger than . The former is not very robust, and is sensitive to fine-tuning of parameters whereas the latter is very robust and is associated with a valley Chern number of 2 (See b). b The Chern number of the first conduction band from valley. Note, the Chern number is defined as long as a direct bandgap is present. c A schematic plot for the insulating (black) regions and the first conduction/valence band isolated region (red/blue) in the TDBG at . The red dot is charge neutrality point (CNP). In the shaded region, strongly correlated physics is expected near integer fillings. Asymmetry between electron and hole dopings is predicted from the theory. d–g Color plots for -factor associated with orbital magnetic effects , , , and single-particle dispersion over the Moiré Brillouin zone for the first conduction band at , where the band is isolated. are in the unit of , and is in the unit of meV. Both and vanish at high symmetric points , , and .

Another crucial difference compared with TBG is the absence of twofold rotational symmetry, which protects the Dirac points in TBG. As a result, the physics of TDBG is controlled by a single narrow band (per spin per valley) rather than two as in TBG. The TDBG Hamiltonian has the following symmetries (i) threefold rotation symmetry , (ii) time-reversal symmetry , and, (iii) mirror reflection about the -axis , which only exists in the absence of vertical electric field, and (iv) SU(2) spin-rotation symmetry. Finally, we assume that in the small angle limit, there is valley-charge-conservation symmetry , arising from the decoupling of Moiré and atomic lattice-scale physics.

In addition, the conduction band within each valley carries a nonzero Chern number. In ordinary condensed matter systems, -symmetry forbids the existence of Chern bands. However, in Moiré systems, Chern bands carrying opposite Chern numbers for opposite valleys can arise due to the valley decoupling. The overall system still satisfies -symmetry, which exchanges the two valleys. Therefore, spontaneous valley polarization would lead to a Chern band without explicitly breaking -symmetry[13,24,26,29]. At , the reflection symmetry enforces for both valleys. At , the conduction band develops a nonvanishing Chern number computed numerically in Fig. 3c which is equal to for the parameter region corresponding to band isolation. The evolution of Chern number as a function of is further confirmed using symmetry indicator (Methods). This can be also understood from the well-known behavior of a AB-stacked bilayer graphene under an electric field. Under the electric field, the bilayer graphene becomes gapped and accumulates opposite Berry curvatures at and valleys, which amounts to a Chern number for each valley.[34-37].

Finally, we discuss the effect of applied magnetic field which influences the single-particle physics in two distinct ways. First, it couples to the electron spin via Zeeman effect leading to the splitting of bands with opposite spin by . Second, it couples to the electron orbital motion leading to modifications in the band structure. For out-of-plane field, the orbital effect arises from the magnetic field coupling to the planar motion of the electron[38,39]. It leads to an energy correction of , with a -dependent -factor satisfying due to time-reversal symmetry ( is a valley index). As shown in Fig. 3f, can be much larger than the Zeeman effect. For in-plane field, the orbital effect arises from coupling to the interlayer motion of electrons. For an in-plane field , we can choose the gauge which does not depend on or , thus preserving the Moiré translation symmetry. The resulting change in the hopping parameters is obtained by the Peierl’s substitution, effectively providing an additional momentum shift of to the hopping connecting layers from to , where is the interlayer separation (see the Methods section). This leads to an energy correction of the form to the leading order in with . The orbital effect due to in-plane field amounts to a very small relative momentum shift . However, it cannot be neglected since it is of the same order of magnitude as the Zeeman effect, (see Fig. 3d, e). In general, the in-plane orbital contribution changes the band dispersion due to its -dependence, whereas the Zeeman effect shifts the entire band uniformly. Moreover, it acts oppositely for different valleys. These properties can be crucial in understanding the effect of in-plane field on the insulating gap and the superconducting temperature (see the Methods section and Supplementary Note 6).

Correlated insulating states

In the band isolation regime, the first conduction band carries a nonzero Chern number as shown in Fig. 3a, b which prevents the existence of exponentially localized Wannier functions[40]. As a result, one cannot construct a Hubbard model for the band unless valley symmetry is broken, or the model is enlarged to include more bands so that the net Chern number is zero. Instead of seeking a complicated real-space description, we discuss the interaction effects in the momentum space, as in the case of quantum-Hall ferromagnetism. One major consequence of the absence of localized Wannier orbitals is the inadequacy of the Mott picture, where the insulating phase is driven by strong repulsion between localized orbitals. Thus, we will use the terminology, correlated insulator to refer to the interaction-driven insulating phase for the following physics.

In order to uncover the nature of the possible correlated insulating states at half and quarter-filling[20], we perform a self-consistent Hartree–Fock mean-field theory similar to the one employed in ref. [8,24]. Below, we sketch the derivation from the microscopic theory, relegating most details to Supplementary Notes 2 and 3. The interacting Hamiltonian in momentum space is given bywhere is the Fourier-transformed screened Coulomb interaction[41,42]. Since the screening coming from the distance between the system and the gate is comparable with the Moiré length scale, the screening length can be important for the interaction effects. The density consists of an intravalley part and an intervalley part , where is the electron creation operator for valley. The latter contribution arises from the small coupling between opposite valleys and gives rise to an intervalley Hund’s coupling term.

The resulting Hamiltonian consists of two parts, , where contains the coupling between intravalley densities , whereas contains the coupling between intervalley densities . Rough estimation for the relative energy scales for and gives and for the experimentally relevant regime. Although is significantly smaller than , it breaks the symmetry of the model down from two independent SU(2) spin-rotation symmetries for each valley to a single SU(2). Thus, it can lift the degeneracy between some symmetry breaking states which are degenerate on the level of the . Indeed, we found that favors the spin alignment between opposite valleys and can be written in the form of intervalley Hund’s coupling as in ref. [24].

Within the self-consistent Hartree–Fock mean-field theory, we consider the order parameter defined asFor a gapped phase, matrix must be a projector, i.e., satisfying for all . Given that there are four flavors of fermions due to spin () and valley () degeneracies, any possible order parameter can be expanded in terms of the generators of SU(4) , which can be grouped based on their symmetry breaking into five categories: (i) only breaks and corresponds to a valley-polarized (VP) state, (ii) breaks spin-rotation symmetry and correspond to a spin-polarized (SP) state. (iii) breaks both spin rotation and time-reversal (but preserve some combination of the two) and corresponds to a spin-valley locked (SVL) state, (iv) breaks valley-charge conservation and corresponds to an intervalley coherent (IVC) state, and (v) breaks both spin rotation and U(1) valley-charge conservation, corresponds to spin-IVC locked (SIVCL) state (see Table 1). We note that any of these orders may break or preserve symmetry depending on its dependence.

Table 1

Symmetry broken states and the remaining symmetries for all possible translation-symmetric gapped states at .

ν=2	Example of the state	Symmetry
SP	↑K_+⊗↑K_−	U(1)_z, U(1)_v, T
VP	↑K_+⊗↓K_+	SU(2), T
SVL	↑K_+⊗↓K_−	U(1)_z, U(1)_v, T^′
IVC	↑K_++e^iθ↑K_−⊗↓K_++e^iθ↓K_−	SU(2), T
SIVCL	↑K_++e^iθ↓K_−⊗↓K_++e^iθ↑K_−	U(1)_z, Z2xv, T

The similar table with the form of and symmetry generators is in the Supplementary Table 1. Here, is the spinless time-reversal squaring to whereas is the spinful time-reversal squaring to (with denoting complex conjugation). denotes spin rotation around the axis by an angle whereas denotes rotation in the valley plane by an angle . Finally, is generated by the combined rotation

The results of the self-consistent Hartree–Fock calculation are summarized in the following (Supplementary Note 3). Restricting ourselves to translation-symmetric gapped states, we find there are five options: SP, VP, SVL, IVC, and SIVCL at half-filling and three options: spin-valley-polarized (SVP), spin-polarized-IVC (SPIVC), and spin-valley-locked-IVC (SVLIVC) at quarter-filling , as in Table 1. By solving the Hartree–Fock self-consistency condition, the ground-state energy and the correlation gap are computed for different states (Fig. 4a). Let us first consider what happens in the absence of intervalley Hund’s coupling. In this case, we find that the SP and SVL states at half-filling and similarly the SPIVC and SVLIVC states at quarter-filling are exactly degenerate since they are related by a spin rotation in one of the valleys. Similarly, due to the enlarged symmetry of the mean-field Hamiltonian, the SP and VP states and the IVC and SIVCL states have the same energy. Thus, we only need to numerically investigate the competition between SP and IVC at half-filling and SVP and SPIVC at quarter-filling. The result of such numerical investigation is shown in Fig. 4a, where we clearly see that SP has a lower energy than that of the IVC in most of the parameter regime. Similar results apply for the competition between SVP and SPIVC at quarter-filling. The correlation-induced gap for the SP state in the band isolation region ranges between 4 and 8 meV (see Fig. 4b).

Fig. 4

The results of the Hartree–Fock calculation. a Color plot (meV) for per electron. b Color plot of self-consistency gap for the SP/VP-state in the band isolated region. (No -term included) c, d Effect of the intervalley Hund’s coupling (-term) on the gap for spin- and valley-polarized phases at half and quater fillings, respectively. At half-filling, -term increases (decreases) (). At quarter-filling, -term reduces the gap to the next-excited state, making the quarter-filled insulator (SP + VP) less stable than the half-filled (SP) one. e The correlated gap for half-filling insulators (SP, VP) as a function of in-plane -field. . Solid lines for SP state and dotted lines for VP-state. Zeeman effect would increase (decrease) for the SP (VP) state with increasing . The valley orbital effect leads to a linear decrease in the gap with field, thus effectively decreasing (increasing) the -factor for the SPS (VP) state.

To understand the reason why IVC order is energetically unfavorable, we can employ the argument of ref. [29] as follows. IVC order between two valleys with opposite Chern number is equivalent after a particle–hole transformation in one of the valleys to superconducting pairing between bands with the same Chern number i.e., a superconductor in a background magnetic field. This means that the order parameter necessarily includes vortices within the Brillouin zone leading to increased energy. A more detailed analytic treatment of the energy competition between SP and IVC is provided in the Supplementary Note 4.

The inclusion of the effect of intervalley Hund’s coupling alters the competition between the phases as follows. First, since the term is ferromagnetic, it lowers the energy of the SP state, favoring the SP state over the VP-state, which is in turn favored over the SVL-state. Second, it lowers the energy of the filled bands for the SP state at half-filling, thus increasing . On the other hand, it reduces the energy of some of the empty bands for the VP-state, reducing (see Fig. 4c, e). The Hund’s coupling term similarly reduces at quarter-filling by lowering the energy of one of the excited states (see Fig. 4d). We note here that the reduction of the correlated gap at quarter-filling relative to that at half-filling may explain why the former is more difficult to observe experimentally compared with the latter and requires the application of a magnetic field[20].

In the presence of an in-plane field, the gap of the SP-phase at half-filling is expected to grow with a slope consistent with the Zeeman factor. However, the orbital effect discussed earlier leads to a reduction in the effective -factor by 20–50% depending on the band structure details (Fig. 4e), which is in agreement with the experimental data[20]. From the numerical calculation, we confirmed that such a reduction in gap also depends on the in-plane field direction, which exhibits threefold periodicity (see the Methods section). Therefore, the orbital effect can be directly verified in a rotating in-plane field setup, where we predict the modulation of the -factor with period in the angle.

Superconductivity

When the correlated insulator is doped away from half-filling, a superconducting phase is observed below 3.5 K[20]. Our proposed scenario for the observed superconductivity is illustrated in Fig. 5a, where pairing takes place between time-reversal partners in opposite valley. Such an intervalley pairing between time-reversal partners has also been proposed[43-45] and observed in transition metal dichalcogenides (TMD)[46]. However, unlike in TMD, where strong spin–orbit coupling implies a locking between spin and valley, here the proposed pairing takes place between the electrons with the same spin. To understand this, we first note that doping a spin-polarized insulator is expected to give rise to a ferromagnetic metal with spin-split Fermi surface. Similar to other ferromagnetic metals[47-49], ferromagnetic spin fluctuations can act as a pairing glue responsible for superconductivity[50]. This motivates the following simplified Hamiltonian,where the spin operator . This Hamiltonian can be obtained within an RPA treatment by identifying the ferromagnetic order as the leading instability in the doped itinerant phase. The ferromagnetic susceptibility is peaked at , which justifies a -independent coupling.

Fig. 5

Spin-triplet superconductivity. a Triplet paring between opposite valleys, and with exact energy match. b Schematic plot for the as a function of -field.

Next, we consider the simplest possible intervalley superconducting pairing function , which is -independent (-wave) within each valley. Note, however, that the overall orbital symmetry incorporating both momentum and valley may still be anti-symmetric, e.g., -wave. For the proposed pairing, is proportional to or corresponding to valley triplet or singlet, respectively. The overall antisymmetry of implies that the former scenario corresponds to a spin-singlet , whereas the latter corresponds to a spin triplet . Here, is the vector which captures the direction of the spin state. To see which of these is the dominant pairing channel, it is useful to decouple the interaction in the pairing channel asWe now assume -independent and decompose it into spin-singlet/velly triplet and spin triplet/valley-singlet . We now usewhere and . This means that the interaction is repuslive in the singlet channel, and attractive in the triplet channel making the latter the dominant pairing channel. A more detailed discussion of these pairing channels within the linearized BCS equation is provided in the Supplementary Note 5.

We highlight here that spin-triplet pairing is only known to occur in liquid He[51] and a few Uranium compounds[47-49], as it requires pairing that varies over the Fermi surface (eg. -wave) which is likely to be energetically unfavorable in typical solids. The existence of the valley degree of freedom here enables us to evade this difficulty and obtain a spin-triplet valley-singlet order parameter even for a -independent interaction.

The experimental consequences of the proposed spin-triplet valley-singlet superconductivity can be investigated by writing the Ginzburg–Landau free energy for the order parameter in the presence of a magnetic field . Restricting ourselves to terms up to quartic order in or , we can write the following free energy functional

Detailed microscopic derivation of the coefficients is provided in the Supplementary Note 6. In the absence of spin–orbit coupling, the order parameter’s spin is expected to align with the magnetic field. Assuming the magnetic field is parallel to the -axis, , we can then writeSubstituting in the free energy (6) and using the fact that yields

One important feature is that which implies the stability of the phase considered.

The free energy (8) leads to the following dependence of the superconducting on the applied fieldThe most remarkable feature of this result is that, for nonzero , initially increases upon the application of magnetic field. This can be understood as follows: for a ferromagnetic metal with weakly spin-split Fermi surfaces, the application of the Zeeman field increases (decreases) the density of states for the majority (minority) spin Fermi surface, leading to a linear increase in for the majority spin with the coefficientwhere is the bandwidth, is the density of states at the Fermi energy, and is the dimensionless magnetic susceptibility (Supplementary Note 6). Similar linear field-dependence of is known in superfluid He[51], indicating independent pairing for each spin species. This behavior is in stark contrast to the monotonic decrease of under increasing -field in a spin-singlet superconductor. One crucial observation here is that seems to depend on several details and is expected to be very small since . Surprisingly, the measured value of is of order 1[20], which suggests the vicinity of a quantum critical point where the scaling of the susceptibility cancels exactly against the other parameters. Indeed, the scaling predicted by Herz-Millis theory in the quantum critical regime for an itinerant ferromagnet[52,53] leads to such cancellation resulting in .

The origin of the quadratic term in Eq. (9) can be understood in terms of the in-plane orbital effect discussed in Sec. IA. First, note that Zeeman splitting cannot break Cooper pairs between aligned spins. Instead, it yields an initial linear increase in followed by saturation at large fields when all the spins are aligned. On the other hand, the in-plane orbital effect can induce pair breaking by mismatching the energies of time-reversal partner states in opposite valleys, resulting in a quadratic decrease in with the applied field whose coefficient is given by (see Supplementary Note 6)where is the direction of the external magnetic field. The average value of over the Fermi surface depends strongly on the filling and the field direction with typical value around 1 (cf. Fig. 3d–f). Using this value, we can make a rough estimate for the in-plane field needed to destroy superconductivity as yielding a value about 3 Teslas, which compares favorably to the experimental value[20]. Furthermore, if we consider an out-of-plane field instead, is on average ~1–2 orders of magnitude larger than , yielding a critical field of which is very close to the experimentally observed result[20].

It is worth noting that the reduction of at large field can also arise from the suppression of ferromagnetic fluctuations responsible for the pairing, as has been observed in the ferromagnetic superconductor UCoGe[54]. Such effects are neglected within our simplified analysis 3, which assumes a constant coupling .

Discussion

In this work, we theoretically investigated the physics of twisted double-bilayer graphene (TDBG), addressing the experimental observations of correlated insulating phases at integer fillings and the neighboring superconductor reported in ref. [20].

First, let us summarize a few important features of the band structure. Due to the absence of a symmetry in TDBG, isolated conduction and valence bands with nonzero valley Chern numbers can exist. Moreover, trigonal warping and particle–hole asymmetry in each bilayer graphene lead to (i) a significant broadening of each band so that they overlap in the absence of a displacement field, and (ii) asymmetry between electron- and hole-doped systems. As a result, the parameter space that can host strongly correlated physics is significantly constrained, and the tunability from displacement field at a particular filling becomes essential to realizing correlated states.

Second, we identified an important role played by the coupling of in-plane field to the orbital motion of the electron in TDBG. Despite being small compared with the bandwidth, this effect is comparable with Zeeman splitting, leading to a modified -factor which compares favorably to the experimental value[20] extracted from the slope of the half-filling gap as a function of in-plane field. Moreover, in our theory, this effect is responsible for the reduction of under an in-plane field by providing the main pair-breaking mechanism when pairing takes place between aligned spins in opposite valleys. The resulting decrease in the superconducting with in-plane field agrees qualitatively with the experimental results.

Furthermore, we have performed a self-consistent Hartree–Fock mean-field calculation to identify the possible symmetry broken correlated insulating states at integer fillings. Our prediction of a spin-polarized ferromagnet at half-filling is consistent with the observed increase in the gap with in-plane field.

Finally, here we have proposed a pairing mechanism based on ferromagnetic fluctuations, which is motivated by the evidence for a ferromagnetic parent insulator. Such a mechanism leads naturally to the spin-triplet pairing suggested by experiments. In addition, we showed that the experimentally observed dependence of on in-plane field suggests that the superconductor emerges in the vicinity to a quantum critical point.

In conclusion, our theoretically established phase diagram for twisted double-bilayer graphene, captures all significant observations of the experiments reported in ref. [20]. This includes single-particle features such as the parameter range for band isolation as well as correlation-induced features including a ferromagnetic insulator at half-filling which leads to a spin-triplet superconductor upon doping. In addition to deepening our understanding of correlated Moiré materials, our results highlight how phases which are rare in conventional solids can be readily realized in this novel and tunable platform.

After completing this work, we noticed two experimental papers[55,56] which are consistent with ref. [20] and theoretical discussion contained here.

Methods

Numerical simulations for single particle

Here, we summarize the numerical methods used to calculate the single-particle physics. First, each bilayer-graphene (BLG) layer is modeled by the following bloch Hamiltonian:which is labeled in the order of , , , . Here, we consider a realistic model of BLG illustrated in Fig. 1. AB stacking means that the -site of the first layer () sits on top of the -site of the second layer (). This gives a small on-site energy for these sites. Here, , where , , and are vectors from -site to -sites. One can expand near aswhere is the distance between carbon atoms. Throughout, we will use the phenomenological parameters extracted from ref. [57]where and are the parameters illustrated in Fig. 1. In addition, the potential difference between the top and bottom graphene layer, is an important parameter in the experiment, which is controlled by the gate voltage difference. For a displacement field strength , AB–AB system’s dielectric constant and the thickness of the BLG/BLG system , .

Next, we couple two layers of AB-stacked bilayer graphenes by Moire hoping terms. As we are interested in the physics near charge neutrality point, we focus on band structures mostly originated near points. In the continuum model approximation[30], Moire bands from valleys decouple; for the Moire band from valley, the Hamiltonian is given bywhere is a 4-components electron creation operator for top/bottom layer with momentum . Here, with denoting the counter-clockwise rotation matrix by angle relative to the -axis. The momenta are given by , , and where . The hopping matrices , are given bywhere are Moiré hopping parameters. One crucial parameter tunable in experiments is displacement field . In Fig. 6, we demonstrated how the band structure evolves with increasing . One can see that the first conduction band becomes isolated in the range of . Furthermore, to illustrate the how the band isolation arises, we plot the energy gap between different bands in Fig. 7. For a smaller value of , gapped regimes in Fig. 7a–c expand in the parameter space of , giving arise to a wider band isolation regime (data available upon request).

Fig. 6

The band structure of the model at and . At , Chern number is exchanged by between the conduction and valence band at three momenta which are located not along the symmetric cut. However, at and , Chern number changes by which can be seen by the gap closing between bands at and points.

Fig. 7

The bandgap (meV) for the range of . Uncolored region implies bands being overlapped. a Gap between the first conduction and valence bands. b Gap between the first and second conduction bands. c Gap between the first and second valence bands.

Chern number

In the main text, we presented Chern number carried by Moire first conduction bands from -valleys. Here, we carefully examine the evolution of Chern. First, at , the reflection symmetry enforces for both valleys as maps the system back to itself without exchanging valleys, but so Berry curvature flips its sign[24]. In the quadratic band approximation limit of BLG, as we increase , the band inversion between conduction and valence bands occurs at the Moiré -point ( for negative ) with a quadratic touching. Thus, Chern number of is exchanged.

Next, let us understand the Chern number evolution in the realistic Hamiltonian with parameters of Eq. (14) along the dotted line in Fig. 3b. With a trigonal warping term, the quadratic band touching point splits into four Dirac cones, three with positive and the other with negative chirality. These three Dirac cones are located at generic momenta, thus would not be observed in the band plot along the high symmetric line. Under the presence of particle–hole asymmetry terms, the degeneracy between four Dirac cones split, and the band inversion would happen first at three Dirac cones, exchanging Chern number by . Then, the band inversion would occur at the center Dirac cone, exchanging Chern number by . In total, it will still change the Chern number by . At larger values of the gate voltage , the band inversion happens between first and second conduction band at point, and the Chern number then changes by (It can change by for other parameter setting), decreasing the Chern number.

This can be further checked by inspecting symmetry indicators[58-60]. There are three -invariant momenta , , and . For a Bloch state with these momenta, rotation symmetry would map the state back to itself with a rotation eigenvalue:where is an angular momentum associated with the Bloch state . Then, the Chern number of the -th band can be determined modulo 3 byThus, by tracking how eigenvalues of the three momenta change with the gating voltage , we can understand how Chern number transition happens in the system. Indeed, the aforementioned scenario can be confirmed. For example, consider a Moiré first conduction band for valley at . At  meV, we start with . At  meV, Chern number changes by but it can be only captured by Berry curvature not by symmetry indicator. At  meV, Chern number changes by , manifested by . At  meV, Chern number again changes by , manifested by . See Fig. 6 for the detail.

Magnetic field effect

Under in-plane magnetic field , one can choose the gauge . Then, the effect of a magnetic field on hopping terms is evaluated via Peierl’s substitution, where the hopping term from to is multiplied by the phase factorsuch thatwhere since is linear function of . Hence, the effect of in-plane field can be included by simply replacing all -dependent matrix elements of Bloch Hamiltonians by as follows (we take ):where is the matrix element connecting layers and ( from bottom to top) in Eq. (15), is the interlayer distance, and is the unit vector in the direction.

Due to its small magnitude relative to the energy gap, it suffices to consider the in-plane orbital effect to first order in pertrubation theory. This amounts to adding the following in-plane orbital term to the single-particle energieswhere is given bywhere is the valley index. Time-reversal symmetry implies that . The in-plane orbital -factor transforms under rotation asprovided that the band is non-degenerate at . This implies that vanishes at any -invariant point. As pointed out in the Results, in general, the in-plane orbital contributions affects the bands very differently from the Zeeman effect. For example, it can distort the Fermi surface when the bands are partially filled in an opposite way in the two valleys which can influence the physical properties, e.g., superconducting (see Supplementary Note 6).

The effect of out-of-plane field on the energy bands is generally more complicated since any gauge choice breaks translation symmetry. As a result, the band picture breaks down for large enough out-of-plane fields where Landau level physics form instead. In the following, we will consider the limit of weak out-of-plane fields which can be treated perturbatively. In this case, the out-of-plane field induces an orbital valley Zeeman effect as pointed out in ref. [38,39] whose -factor is given byIn summary, the single-particle energies has the following dependence on magnetic fieldwhere is the electron spin operator (which is for up/down spins) and . The valley orbital -factor is defined asWe have also assumed that the spin-quantization axis is parallel to the field.