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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{{\rm{v}}}=\pm 2$$\end{document}Cv=±2. 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.

T he recent discovery of correlated insulating states and superconductivity in twisted bilayer graphene (TBG) [1][2][3][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][8][9][10][11][12][13][14][15][16][17][18][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 C 2 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 g factor from 2 20 , as well as the effect of in-plane field on superconducting T c .
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 C 2 -breaking substrate potentials 13,19,24,[27][28][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 momentumindependent within each valley, is favored. We investigate the consequences of such scenario and show it can be used to explain the measured dependence of T c 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 AA and AB sites, parametrized by ðw 0 ; w 1 Þ 21,22 . In the original Bistritzer-Macdonald model, w 0 and w 1 are taken to be equal 30 . However, in a realistic twisted model, the ratio r w 0 =w 1 is smaller than one due to the lattice relaxation which expands (shrinks) AB (AA) regions. In TBG, r is taken to be~0.75 for the first magic angle 21,22 . Here, we similarly include lattice relaxations by taking r 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 ν ¼ ± 4 filling. In this work, we take ðw 0 ; w 1 Þ ¼ ð88; 100Þ meV corresponding to r ¼ 0:88. For different values of ðw 0 ; w 1 Þ, we obtained qualitatively similar features (Methods).
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 θ % 1:05 24 . However, once trigonal warping (γ 3 ) and particle-hole asymmetry (γ 4 ) 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 U 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 γ 4 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 .
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 C 3 , (ii) time-reversal symmetry T , and, (iii) mirror reflection about the x-axis M y , which only exists in the    We observe two seperate isolation regions for θ smaller or larger than 1:1 . 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 K þ 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 θ ¼ 1:33 . 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 g-factor associated with orbital magnetic effects g x þ ðkÞ, g y þ ðkÞ, g z þ ðkÞ, and single-particle dispersion ξ þ ðkÞ over the Moiré Brillouin zone for the first conduction band at ðθ; UÞ ¼ ð1:33 ; 60 meVÞ, where the band is isolated. g x;y;z ðkÞ are in the unit of μ B , and ξðkÞ is in the unit of meV. Both g x and g y vanish at high symmetric points Γ, K 1 , and K 2 .
In addition, the conduction band within each valley carries a nonzero Chern number. In ordinary condensed matter systems, T -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 T -symmetry, which exchanges the two valleys. Therefore, spontaneous valley polarization would lead to a Chern band without explicitly breaking T -symmetry 13,24,26,29 . At U ¼ 0, the reflection symmetry M y enforces C ¼ 0 for both valleys. At U ≠ 0, the conduction band develops a nonvanishing Chern number computed numerically in Fig. 3c which is equal to ± 2 for the parameter region corresponding to band isolation. The evolution of Chern number as a function of U 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 K þ and K À valleys, which amounts to a Chern number C v ¼ ± 2 for each valley. [34][35][36][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 2μ B B. 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 μ B g z τ ðkÞB z , with a k-dependent g-factor g z τ ðkÞ satisfying g z Àτ ðÀkÞ ¼ Àg z τ ðkÞ due to time-reversal symmetry (τ is a valley index). As shown in Fig. 3f, g z τ ðkÞ 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 B, we can choose the gauge AðzÞ ¼ Àz B which does not depend on x or y, 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 À e _ ðl þ mÞd 2 e z B to the hopping connecting layers from l to m, where d is the interlayer separation (see the Methods section). This leads to an energy correction of the form μ B ðg x τ ðkÞB x þ g y τ ðkÞB y Þ to the leading order in B with g x;y Àτ ðÀkÞ ¼ Àg x;y τ ðkÞ. The orbital effect due to in-plane field amounts to a very small relative momentum shift $ eda _ % 10 À5 . However, it cannot be neglected since it is of the same order of magnitude as the Zeeman effect, ev F d Fig. 3d, e). In general, the in-plane orbital contribution changes the band dispersion due to its k-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 inplane 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 realspace 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 selfconsistent 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 by where VðqÞ 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ρðqÞ consists of an intravalley part ρ þ $ c y ± c ± and an intervalley part ρ À $ c y ± c Ç , where c y ± is the electron creation operator for K ± 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 H 0 contains the coupling between intravalley densities ρ þ ρ þ , whereas H J contains the coupling between intervalley densities ρ À ρ À . Rough estimation for the relative energy scales for H 0 and H J gives V 0 $ 35 meV and J $ 0:6 meV for the experimentally relevant regime. Although H J is significantly smaller than H 0 , 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 H 0 . Indeed, we found that H J 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 as For a gapped phase, matrix MðkÞ must be a projector, i.e., MðkÞ 2 ¼ MðkÞ satisfying tr MðkÞ ¼ ν for all k. Given that there are four flavors of fermions due to spin (σ) and valley (τ) degeneracies, any possible order parameter M can be expanded in terms of the generators of SU(4) σ i τ j , which can be grouped based on their symmetry breaking into five categories: (i) fσ 0 τ z g only breaks T and corresponds to a valley-polarized (VP) state, (ii) fσ x;y;z τ 0 g breaks spin-rotation symmetry and correspond to a spin-polarized (SP) state. (iii) fσ x;y;z τ z g breaks both spin rotation and time-reversal (but preserve some combination of the two) and corresponds to a spin-valley locked (SVL) state, (iv) fσ 0 τ x;y g breaks Uð1Þ valley-charge conservation and corresponds to an intervalley coherent (IVC) state, and (v) fσ x;y;z τ x;y g breaks both spin rotation and U(1) v valley-charge conservation, corresponds to spin-IVC locked (SIVCL) state (see Table 1). We note that any of these orders may break or preserve C 3 symmetry depending on its k dependence.
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 ν ¼ 2 and three options: spin-valley-polarized (SVP), spinpolarized-IVC (SPIVC), and spin-valley-locked-IVC (SVLIVC) at quarter-filling ν ¼ 1; 3, as in Table 1. By solving the Hartree-Fock self-consistency condition, the ground-state energy E 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).
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 C 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 jCj 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 Δ SP . On the other hand, it reduces the energy of some of the empty bands for the VP-state, reducing Δ VP (see Fig. 4c, e). The Hund's coupling term similarly reduces Δ SVP 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 g ¼ 2 factor. However, the orbital effect discussed earlier leads to a reduction in the effective g-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 inplane 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 g-factor with period 2π=3 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][44][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 spinpolarized insulator is expected to give rise to a ferromagnetic metal with spin-split Fermi surface. Similar to other ferromagnetic metals [47][48][49] , ferromagnetic spin fluctuations can act as a pairing glue responsible for superconductivity 50 . This motivates the following simplified Hamiltonian, H ¼ X k;τ;σ c y σ;τ;k ξ σ;τ;k c σ;τ;k À g where the spin operator S a q ¼ P k;τ;σ;σ 0 c y σ;τ;kþq σ a σ;σ 0 c σ 0 ;τ;k . 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 q ¼ 0, which justifies a k-independent coupling.
Next, we consider the simplest possible intervalley superconducting pairing function Δ, which is k-independent (s-wave) within each valley. Note, however, that the overall orbital symmetry incorporating both momentum and valley may still be anti-symmetric, e.g., p-wave. For the proposed pairing, Δ is proportional to τ x or τ y corresponding to valley triplet or singlet, respectively. The overall antisymmetry of Δ implies that the former scenario corresponds to a spin-singlet iσ y , whereas the latter corresponds to a spin triplet iσ y d Á σ. Here, d 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 as We now assume k-independent Δ and decompose it into spinsinglet/velly triplet Δ s and spin triplet/valley-singlet Δ t . We now use where λ t ¼ 1 and λ s ¼ À3. 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 Table 1 Symmetry broken states and the remaining symmetries for all possible translation-symmetric gapped states at ν ¼ 1; 2; 3.
The similar table with the form of MðkÞ and symmetry generators is in the Supplementary Table 1. Here, T is the spinless time-reversal T ¼ τ x K squaring to þ1 whereas T 0 is the spinful time-reversal T 0 ¼ iσ y T squaring to À1 (with K denoting complex conjugation). Uð1Þ θ x; y; z ¼ e iθσ x; y; z =2 denotes spin rotation around the x; y; z axis by an angle θ whereas Uð1Þ θ v ¼ e iθτ z =2 denotes rotation in the valley x À y plane by an angle θ. Finally, Z z;v 2 is generated by the combined rotation We highlight here that spin-triplet pairing is only known to occur in liquid He 3 51 and a few Uranium compounds [47][48][49] , as it requires pairing that varies over the Fermi surface (eg. p-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 k-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 Δ ¼ τ y σ y d Á σ in the presence of a magnetic field B. Restricting ourselves to terms up to quartic order in d or B, we can write the following free energy functional Detailed microscopic derivation of the coefficients a; b; c; κ; α; η 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 z-axis, B ¼ Be z , we can then write Substituting in the free energy (6) and using the fact that η ¼ Àα=2 yields One important feature is that α > 0 which implies the stability of the phase considered.
The free energy (8) leads to the following dependence of the superconducting T c on the applied field The most remarkable feature of this result is that, for nonzero a, T c 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 T c for the majority spin with the coefficient where Λ is the bandwidth, Nð0Þ is the density of states at the Fermi energy, and χ is the dimensionless magnetic susceptibility (Supplementary Note 6). Similar linear field-dependence of T c is known in superfluid He 3 51 , indicating independent pairing for each spin species. This behavior is in stark contrast to the monotonic decrease of T c under increasing B-field in a spinsinglet superconductor. One crucial observation here is that a seems to depend on several details and is expected to be very small since T c ( Nð0Þ N 0 ð0Þ $ ϵ F . Surprisingly, the measured value of a 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 χ $ ϵ F =ðTlogTÞ predicted by Herz-Millis theory in the quantum critical regime for an itinerant ferromagnet 52,53 leads to such cancellation resulting in a $ 1. 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 T c ðBÞ 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 T c with the applied field whose coefficient is given by (see Supplementary Note 6 where e B is the direction of the external magnetic field. The average value of ðe B Á g þ ðkÞÞ 2 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 superconductiv- yielding a value about 3 Teslas, which compares favorably to the experimental value 20 . Furthermore, if we consider an out-of-plane field instead, jg z j is on average~1-2 orders of magnitude larger than jg x;y j, yielding a critical field of $ 0:1T which is very close to the experimentally observed result 20 .
It is worth noting that the reduction of T c 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 g.

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 C 2 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 g-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 T c 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 T c 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 T c 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 A 1 , B 1 , A 2 , B 2 . Here, we consider a realistic model of BLG illustrated in Fig. 1. AB stacking means that the A-site of the first layer (A 1 ) sits on top of the B-site of the second layer (B 2 ). This gives a small on-site energy Δ for these sites. Here, f ðkÞ P l e ÀikÁδ l , where δ 1 ¼ að0; À1Þ, δ 2 ¼ aðÀ ffiffi ffi 3 p =2; 1=2Þ, and δ 3 ¼ að ffiffi ffi 3 p =2; 1=2Þ are vectors from B-site to A-sites. One can expand f ðkÞ near K ± ¼ ± ð4π=3 ffiffi ffi 3 p a; 0Þ as where a is the distance between carbon atoms. Throughout, we will use the phenomenological parameters extracted from ref. 57 ðγ 0 ; γ 1 ; γ 3 ; γ 4 ; ΔÞ ¼ ð2610; 361; 283; 138; 15Þ meV; ð14Þ where γ 0;1;3;4 and Δ are the parameters illustrated in Fig. 1. In addition, the potential difference between the top and bottom graphene layer, U is an important parameter in the experiment, which is controlled by the gate voltage difference. For a displacement field strength D, AB-AB system's dielectric constant ϵ and the thickness of the BLG/BLG system d, U ¼ ϵ À1 D Á d.
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 K ± points. In the continuum model approximation 30 , Moire bands from K ± valleys decouple; for the Moire band from K þ valley, the Hamiltonian is given by where c y k;þ;t=b is a 4-components electron creation operator for top/bottom layer with momentum K þ þ k. Here, h θ ðkÞ ¼ hðR Àθ kÞ with R θ denoting the counterclockwise rotation matrix by angle θ relative to the x-axis. The momenta q 0;1;2 are given by q 0 ¼ R θ=2 K À R Àθ=2 K ¼ 8π sinðθ=2Þ 3 ffiffi 3 p a ð0; À1Þ, q 1 ¼ R ϕ q 0 , and q 2 ¼ R Àϕ q 0 where ϕ ¼ 2π=3. The hopping matrices T n , n ¼ 0; 1; 2 are given by where w 0 ; w 1 are Moiré hopping parameters. One crucial parameter tunable in experiments is displacement field U. In Fig. 6, we demonstrated how the band structure evolves with increasing U. One can see that the first conduction band becomes isolated in the range of U 2 ½40; 80. 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 r, gapped regimes in Fig. 7a-c expand in the parameter space of ðθ; UÞ, giving arise to a wider band isolation regime (data available upon request).
Chern number. In the main text, we presented Chern number carried by Moire first conduction bands from K ± -valleys. Here, we carefully examine the evolution of Chern. First, at U ¼ 0, the reflection symmetry M y enforces C ¼ 0 for both valleys as M y maps the system back to itself without exchanging valleys, but k y 7 ! À k y so Berry curvature flips its sign 24 . In the quadratic band approximation limit of BLG, as we increase U, the band inversion between conduction and valence bands occurs at the Moiré K 2 -point (K 1 for negative U) with a quadratic touching. Thus, Chern number of ± 2 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 ± 3. Then, the band inversion would occur at the center Dirac cone, exchanging Chern number by Ç1. In total, it will still change the Chern number by ± 2. At  larger values of the gate voltage U, the band inversion happens between first and second conduction band at Γ point, and the Chern number then changes by Ç1 (It can change by Ç2 for other parameter setting), decreasing the Chern number. This can be further checked by inspecting symmetry indicators [58][59][60] . There are three C 3 -invariant momenta Γ, K, and K 0 . For a Bloch state with these momenta, C 3 rotation symmetry would map the state back to itself with a rotation eigenvalue: where L n;k is an angular momentum associated with the Bloch state k; n j i. Then, the Chern number of the n-th band can be determined modulo 3 by C n L n;Γ þ L n;K 1 þ L n;K 2 mod 3 ð18Þ Thus, by tracking how C 3 eigenvalues of the three momenta change with the gating voltage U, 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 K þ valley at θ ¼ 1:33 . At U ¼ 0 meV, we start with ðn Γ ; n K 1 ; n K 2 Þ ¼ ð0; 1; À1Þ. At U ¼ 14 meV, Chern number changes by þ3 but it can be only captured by Berry curvature not by symmetry indicator. At U ¼ 30 meV, Chern number changes by À1, manifested by n K 2 : À1 7 ! 1. At U ¼ 90 meV, Chern number again changes by À1, manifested by n Γ : 0 7 ! À1. See Fig. 6 for the detail.
since AðzÞ is linear function of z. Hence, the effect of in-plane field can be included by simply replacing all k-dependent matrix elements of Bloch Hamiltonians by k þ α as follows (we take c k ¼ P R e ÀikÁR c R ): where H l;m is the matrix element connecting layers l and m (l; m ¼ 0; ; 3 from bottom to top) in Eq. (15), d ¼ 3:42 A is the interlayer distance, and e z is the unit vector in the z 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 energies ξ n;τ ðk; BÞ ¼ ξ n;τ ðkÞ þ μ B g xy n;τ ðkÞ Á B where g xy n;τ ðkÞ is given by g xy n;τ ðkÞ ¼ 1 μ B hψ n;τ ðkÞj∇ B H τ ðk; BÞj B¼0 jψ n;τ ðkÞi; ð23Þ where τ is the valley index. Time-reversal symmetry implies that g xy n;τ ðÀkÞ ¼ Àg xy n;Àτ ðkÞ. The in-plane orbital g-factor transforms under C 3 rotation as g xy n;τ ðR ± 2π=3 kÞ ¼ R Ç2π=3 g xy n;τ ðkÞ ð24Þ provided that the band n is non-degenerate at k. This implies that g xy n;τ ðkÞ vanishes at any C 3 -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 T c (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-ofplane 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 g-factor is given by g z n;τ ðkÞ ¼ À 4m _ 2 Im X l≠n hn; τj∂ k x H τ jl; τihl; τ; j∂ k y H τ jni ϵ n;τ;k À ϵ l;τ;k : ð25Þ In summary, the single-particle energies has the following dependence on magnetic field ξ n;σ;τ ðk; BÞ ¼ ξ n;τ ðkÞ þ μ B ðgσ Á B þ g n;τ ðkÞ Á BÞ; ð26Þ where σ is the electron spin operator (which is ±1=2 for up/down spins) and τ ¼ ± . The valley orbital g-factor is defined as g n;τ ðkÞ ¼ ðg xy n;τ ðkÞ; g z n;τ ðkÞÞ: ð27Þ We have also assumed that the spin-quantization axis is parallel to the field.

Data availability
All relevant data are available from the authors upon reasonable request.

Code availability
All relevant codes are available from the authors upon reasonable request.