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

TL;DR

This paper theoretically explores twisted double bilayer graphene, revealing that flat bands emerge under a displacement field, leading to ferromagnetic insulating states and spin-triplet superconductivity, with results aligning with recent experimental observations.

Contribution

The study introduces a phase diagram for TDBLG showing flat band emergence under displacement fields and predicts spin-triplet superconductivity driven by ferromagnetic fluctuations.

Findings

Flat bands appear only with a vertical displacement field D.

Half-filled states are ferromagnetic insulators with valley Chern number C_v=2.

Response to magnetic fields matches recent experimental data.

Abstract

Two monolayers of graphene twisted by a small `magic' angle exhibit nearly flat bands leading to correlated electronic states and superconductivity, whose precise nature including possible broken symmetries, remain under debate. Here we theoretically study a related but different system with reduced symmetry - twisted {\em double} bilayer graphene (TDBLG), consisting of {\em two} Bernal stacked bilayer graphene sheets, 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 $D$. We construct a phase diagram as a function of twist angle and $D$, incorporating interactions via a Hartree-Fock approximation. At half filling, ferromagnetic insulators are stabilized, typically with valley Chern number $C_v=2$. Ferromagnetic fluctuations in the metallic state are argued to lead to spin triplet superconductivity from pairing between electrons in opposite valleys. Response of these states to a magnetic field applied either perpendicular or parallel to the graphene sheets is obtained, and found to compare favorably with a recent experiment. We highlight a novel orbital effect arising from in-plane fields that can exceed the Zeeman effect and plays an important role in interpreting experiments.

Figures

Tables

Table 1. Table 1: Examples of symmetry broken states and corresponding remaining symmetries for all possible translation-symmetric gapped states at half ν = 2 𝜈 2 \nu=2 and quarter ν = 1 𝜈 1 \nu=1 fillings. The similar table with the form of M ​ ( 𝒌 ) 𝑀 𝒌 M({\bm{k}}) and symmetry generators is in Supplementary material 3. Here, 𝒯 𝒯 \mathcal{T} is the spinless time-reversal 𝒯 = τ x ​ 𝒦 𝒯 subscript 𝜏 𝑥 𝒦 \mathcal{T}=\tau_{x}\mathcal{K} squaring to + 1 1 +1 whereas 𝒯 ′ superscript 𝒯 ′ \mathcal{T}^{\prime} is the spinful time-reversal 𝒯 ′ = i ​ σ y ​ 𝒯 superscript 𝒯 ′ 𝑖 subscript 𝜎 𝑦 𝒯 \mathcal{T}^{\prime}=i\sigma_{y}\mathcal{T} squaring to − 1 1 -1 (with 𝒦 𝒦 \mathcal{K} denoting complex conjugation). U ​ ( 1 ) x , y , z θ = e i ​ θ ​ σ x , y , z / 2 𝑈 superscript subscript 1 𝑥 𝑦 𝑧 𝜃 superscript 𝑒 𝑖 𝜃 subscript 𝜎 𝑥 𝑦 𝑧 2 U(1)_{x,y,z}^{\theta}=e^{i\theta\sigma_{x,y,z}/2} denotes spin rotation around the x , y , z 𝑥 𝑦 𝑧 x,y,z axis by an angle θ 𝜃 \theta whereas U ​ ( 1 ) v θ = e i ​ θ ​ τ z / 2 𝑈 superscript subscript 1 𝑣 𝜃 superscript 𝑒 𝑖 𝜃 subscript 𝜏 𝑧 2 U(1)_{v}^{\theta}=e^{i\theta\tau_{z}/2} denotes rotation in the valley x − y 𝑥 𝑦 x-y plane by an angle θ 𝜃 \theta . Finally, ℤ 2 z , v superscript subscript ℤ 2 𝑧 𝑣 \mathbb{Z}_{2}^{z,v} is generated by the combined rotation U ​ ( 1 ) z π ​ U ​ ( 1 ) v π 𝑈 superscript subscript 1 𝑧 𝜋 𝑈 superscript subscript 1 𝑣 𝜋 U(1)_{z}^{\pi}U(1)_{v}^{\pi} .

$\nu =2$	Example of the state	Sym. Gen.
SP	$|\uparrow {𝑲}_{+}⟩\otimes |\uparrow {𝑲}_{-}⟩$	$U{\left(1\right)}_z$, $U{\left(1\right)}_v$, $𝒯$
VP	$|\uparrow {𝑲}_{+}⟩\otimes |\downarrow {𝑲}_{+}⟩$	$SU\left(2\right)$, $𝒯$
SVL	$|\uparrow {𝑲}_{+}⟩\otimes |\downarrow {𝑲}_{-}⟩$	$U{\left(1\right)}_z$, $U{\left(1\right)}_v$, ${𝒯}^{\prime }$
IVC	$\begin{array}{c}(|\uparrow {𝑲}_{+}⟩+e^{i\theta }|\uparrow {𝑲}_{-}⟩)\otimes \\ \left(|\downarrow {𝑲}_{+}⟩+e^{i\theta }|\downarrow {𝑲}_{-}⟩\right)\end{array}$	$SU\left(2\right)$, $𝒯$
SIVCL	$\begin{array}{c}(|\uparrow {𝑲}_{+}⟩+e^{i\theta }|\downarrow {𝑲}_{-}⟩)\otimes \\ \left(|\downarrow {𝑲}_{+}⟩+e^{i\theta }|\uparrow {𝑲}_{-}⟩\right)\end{array}$	$U{\left(1\right)}_z$, ${\mathbb{Z}}_2^{xv}$, $𝒯$
$\nu =1,3$	Example of the state	Symmetry
SVP	$|\uparrow {𝑲}_{+}⟩$	$U{\left(1\right)}_z$, $U{\left(1\right)}_v$
SPIVC	$|\uparrow {𝑲}_{+}⟩+e^{i\theta }|\uparrow {𝑲}_{-}⟩$	$U{\left(1\right)}_z$, $𝒯$
SVLIVC	$|\uparrow {𝑲}_{+}⟩+e^{i\theta }|\downarrow {𝑲}_{-}⟩$	${\mathbb{Z}}_2^{zv}$, ${𝒯}^{\prime }$

Table 2. Table 2: Examples of order parameter M ​ ( 𝒌 ) 𝑀 𝒌 M({\bm{k}}) and corresponding independent generators of preserved symmetries for all possible translation-symmetric gapped states at half ν = 2 𝜈 2 \nu=2 and quarter ν = 1 𝜈 1 \nu=1 fillings. Note that the M ​ ( 𝒌 ) 𝑀 𝒌 M({\bm{k}}) can take a more general form. For example, in IVC or SIVCL, τ x subscript 𝜏 𝑥 \tau_{x} can be replaced by c x ​ τ x + c y ​ τ y subscript 𝑐 𝑥 subscript 𝜏 𝑥 subscript 𝑐 𝑦 subscript 𝜏 𝑦 c_{x}\tau_{x}+c_{y}\tau_{y} with c x 2 + c y 2 = 1 superscript subscript 𝑐 𝑥 2 superscript subscript 𝑐 𝑦 2 1 c_{x}^{2}+c_{y}^{2}=1 . Also, for any spin-polarized state, σ z subscript 𝜎 𝑧 \sigma_{z} can be replaced by any 𝝈 = sin ⁡ θ ​ cos ⁡ ϕ ​ σ x + sin ⁡ θ ​ sin ⁡ ϕ ​ σ y + cos ⁡ θ ​ σ z 𝝈 𝜃 italic-ϕ subscript 𝜎 𝑥 𝜃 italic-ϕ subscript 𝜎 𝑦 𝜃 subscript 𝜎 𝑧 \bm{\sigma}=\sin\theta\cos\phi\sigma_{x}+\sin\theta\sin\phi\sigma_{y}+\cos\theta\sigma_{z} . Here, τ x ​ 𝒦 subscript 𝜏 𝑥 𝒦 \tau_{x}{\cal K} is a spinless time-reversal, where 𝒦 𝒦 \cal K is an anti-unitary symmetry. Caveat: For SVLIVC (which is like SVL+SIVCL) state at ν = 1 , 3 𝜈 1 3 \nu=1,3 , only a certain product structure (in this case spin S z subscript 𝑆 𝑧 S_{z} -locked SVL and spin S x subscript 𝑆 𝑥 S_{x} -locked SIVCL) would be allowed.

$\nu =2$	Example of $M\left(𝒌\right)$	Sym. Gen.
SP	$\left(1+{\sigma }_z{\tau }_0\right)/2$	${\sigma }_z$, ${\tau }_z$, ${\tau }_x𝒦$
VP	$\left(1+{\sigma }_0{\tau }_z\right)/2$	${\sigma }_z$, ${\sigma }_x$, ${\tau }_x𝒦$
SVL	$\left(1+{\sigma }_z{\tau }_z\right)/2$	${\sigma }_z$, ${\tau }_z$, ${\sigma }_x{\tau }_x𝒦$
IVC	$\left(1+{\sigma }_0{\tau }_x\right)/2$	${\sigma }_z$, ${\sigma }_x$, ${\tau }_x𝒦$
SIVCL	$\left(1+{\sigma }_x{\tau }_x\right)/2$	${\sigma }_x$, ${\sigma }_z{\tau }_z$, ${\tau }_x𝒦$
$\nu =1,3$	Example of $M\left(𝒌\right)$	Sym. Gen.
SVP	$\begin{array}{c}\left(1+{\sigma }_z{\tau }_0\right)\left(1+{\sigma }_0{\tau }_z\right)/4\end{array}$	${\sigma }_z$, ${\tau }_z$
SPIVC	$\begin{array}{c}\left(1+{\sigma }_z{\tau }_0\right)\left(1+{\sigma }_0{\tau }_x\right)/4\end{array}$	${\sigma }_z$, ${\tau }_x𝒦$
SVLIVC	$\begin{array}{c}\left(1+{\sigma }_z{\tau }_z\right)\left(1+{\sigma }_x{\tau }_x\right)/4\end{array}$	${\sigma }_z{\tau }_z$, ${\sigma }_x{\tau }_x𝒦$

Equations

$${\mathcal{H}}_{\text{int}}=\frac{1}{2\text{Vol}}\sum _{\mathbf{q}}\hat{\rho }(\mathbf{q})V(\mathbf{q})\hat{\rho }(-\mathbf{q}),$$

$$⟨c_{\sigma ,\tau }^{†}(\mathbf{k})c_{{\sigma }^{\prime },{\tau }^{\prime }}({\mathbf{k}}^{\prime })⟩=M_{\sigma \tau ,{\sigma }^{\prime }{\tau }^{\prime }}(\mathbf{k}){\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }}.$$

$$\mathcal{H}=\sum _{\mathbf{k},\tau ,\sigma }{c}_{\sigma ,\tau ,\mathbf{k}}^{†}{\xi }_{\sigma ,\tau ,\mathbf{k}}{c}_{\sigma ,\tau ,\mathbf{k}}-g\sum _{\mathbf{q}}{\mathbf{S}}_{\mathbf{q}}\cdot {\mathbf{S}}_{-\mathbf{q}},$$

$${\mathcal{H}}_{\mathrm{int}}=-g\sum _{\mathbf{k},\mathbf{q}}\mathrm{tr}(\mathbf{\sigma }{\Delta }_{\mathbf{k}})\cdot ({\mathbf{\sigma }}^T{\Delta }_{\mathbf{k}+\mathbf{q}}^{†})$$

$$\mathbf{\sigma }\cdot ({\Delta }_{t,s}{\mathbf{\sigma }}^T)={\lambda }_{t,s}{\Delta }_{t,s},$$

$$\begin{gathered} F=\kappa [(T-T_c+b({\mu }_B\mathbf{B}{)}^2)\mathbf{d}\cdot {\mathbf{d}}^{\ast }+ia{\mu }_B\mathbf{B}\cdot (\mathbf{d}\times {\mathbf{d}}^{\ast }) \\ +c{\mu }_B^2|\mathbf{B}\cdot \mathbf{d}{|}^2+\alpha (\mathbf{d}\cdot {\mathbf{d}}^{\ast }{)}^4+\eta |\mathbf{d}\cdot \mathbf{d}{|}^4] \end{gathered}$$

$$\mathbf{d}=(\frac{{\Delta }_{\uparrow \uparrow }+{\Delta }_{\downarrow \downarrow }}{2},\frac{{\Delta }_{\uparrow \uparrow }-{\Delta }_{\downarrow \downarrow }}{2i},0)$$

$$F=\frac{\kappa }{2}\sum _{s=\uparrow ,\downarrow }{F}_s$$

$$F_s=|{\Delta }_{ss}{|}^2(T-T_c-{\sigma }_{s}a{\mu }_{B}B+b({\mu }_{B}B{)}^2)+\frac{\alpha }{2}|{\Delta }_{ss}{|}^4$$

$$T_{c,\uparrow /\downarrow }(B)=T_c\pm a{\mu }_{B}B-b({\mu }_{B}B{)}^2.$$

$$a=2\chi T_c\frac{N^{\prime }(0)}{N(0)}\ln \frac{\Lambda }{T_c}$$

$$b=\frac{1}{T_c}{\int }_{\mathrm{FS}}d\mathbf{k}({\mathbf{e}}_{\mathbf{B}}\cdot {\mathbf{g}}_{+,\mathbf{k}}{)}^2$$

$$h_{\bm{k}}=\matrixquantity(U_{1}+\Delta&-\gamma_{0}f({\bm{k}})&\gamma_{4}f^{*}({\bm{k}})&\gamma_{1}\,\,\\ -\gamma_{0}f^{*}({\bm{k}})&U_{1}&\gamma_{3}f({\bm{k}})&\gamma_{4}f^{*}({\bm{k}})\\ \gamma_{4}f({\bm{k}})&\gamma_{3}f^{*}({\bm{k}})&U_{2}&-\gamma_{0}f({\bm{k}})\\ \gamma_{1}&\gamma_{4}f({\bm{k}})&-\gamma_{0}f^{*}({\bm{k}})&U_{2}+\Delta),$$

$$f({\mathbf{K}}_{\pm }+\mathbf{k})=\frac{3}{2}(\mp k_x+i{k}_y)a,$$

$$({\gamma }_0,{\gamma }_1,{\gamma }_3,{\gamma }_4,\Delta )=(2610,361,283,138,15)\text{ meV}.$$

$$H_{+}$$

$$\displaystyle+\sum_{n}\quantity(T_{n}c^{\dagger}_{{\bm{k}}+{\bm{q}}_{n},+,b}c^{~{}}_{{\bm{k}},+,t}+T^{\dagger}_{n}c^{\dagger}_{{\bm{k}},+,t}c^{~{}}_{{\bm{k}}+{\bm{q}}_{n},+,b})\Big{]},$$

$$T_{n}=\matrixquantity(0&1\\ 0&0)_{\textrm{layer}}\otimes\quantity(w_{0}+w_{1}e^{2\pi n\sigma_{3}/3}\sigma_{1}e^{-2\pi n\sigma_{3}/3})_{\textrm{sublattice}},$$

$$R_{2\pi /3}|\mathbf{k},n⟩=e^{2\pi i{L}_{n,\mathbf{k}}/3}|\mathbf{k},n⟩,\mathbf{k}=K_1,K_2,\Gamma$$

$$C_n\equiv L_{n,\Gamma }+L_{n,K_1}+L_{n,K_2}\mathrm{mod}3$$

$$e^{i\frac{q}{\hslash }{\int }_{\mathbf{R}}^{\mathbf{R}+\mathbf{\delta }}d\mathbf{r}\cdot \mathbf{A}(\mathbf{z})}=e^{-i\frac{e}{\hslash }{\mathbf{\delta }}_{xy}\cdot \left[\text{quantity}({\mathbf{R}}_z+\frac{{\mathbf{\delta }}_z}{2})\times \mathbf{B}\right]},$$

$$\sum _{\mathbf{R},\mathbf{\delta }}{e}^{i\frac{q}{\hslash }{\int }_{\mathbf{R}}^{\mathbf{R}+\mathbf{\delta }}d\mathbf{r}\cdot \mathbf{A}(\mathbf{z})}{c}_{\mathbf{R}+\delta }^{†}{c}_{\mathbf{R}}=\sum _{\mathbf{k},\mathbf{\delta }}{e}^{-i\text{quantity}(\mathbf{k}+\mathbf{\alpha })\cdot \mathbf{\delta }}{c}_{\mathbf{k}}^{†}{c}_{\mathbf{k}},$$

$${\mathcal{H}}_{l,m}(\mathbf{k},\mathbf{B})={\mathcal{H}}_{l,m}\left(\mathbf{k}-\frac{e}{\hslash }\frac{(l+m)d}{2}{\mathbf{e}}_z\times \mathbf{B}\right)$$

$${\xi }_{n,\tau }(\mathbf{k},\mathbf{B})={\xi }_{n,\tau }(\mathbf{k})+{\mu }_B{\mathbf{g}}_{n,\tau }^{xy}(\mathbf{k})\cdot \mathbf{B}$$

$${\mathbf{g}}_{n,\tau }^{xy}(\mathbf{k})=\frac{1}{{\mu }_B}⟨{\psi }_{n,\tau }(\mathbf{k})|{\nabla }_{\mathbf{B}}{\mathcal{H}}_{\tau }(\mathbf{k},\mathbf{B}){|}_{\mathbf{B}=0}|{\psi }_{n,\tau }(\mathbf{k})⟩$$

$${\mathbf{g}}_{n,\tau }^{xy}(R_{\pm 2\pi /3}\mathbf{k})=R_{\mp 2\pi /3}{\mathbf{g}}_{n,\tau }^{xy}(\mathbf{k})$$

$$g_{n,\tau }^z(\mathbf{k})=-\frac{4m}{{\hslash }^2}\text{imaginary}\sum _{l\ne n}\frac{⟨n,\tau |{\partial }_{k_x}{\mathcal{H}}_{\tau }|l,\tau ⟩⟨l,\tau ,|{\partial }_{k_y}{\mathcal{H}}_{\tau }|n⟩}{{ϵ}_{n,\tau ,\mathbf{k}}-{ϵ}_{l,\tau ,\mathbf{k}}}$$

$${\xi }_{n,\mathbf{\sigma },\tau }(\mathbf{k},\mathbf{B})={\xi }_{n,\tau }(\mathbf{k})+{\mu }_B\text{quantity}(g\mathbf{\sigma }\cdot \mathbf{B}+{\mathbf{g}}_{n,\tau }(\mathbf{k})\cdot \mathbf{B}),$$

$${\mathbf{g}}_{n,\tau }(\mathbf{k})=({\mathbf{g}}_{n,\tau }^{xy}(\mathbf{k}),g_{n,\tau }^z(\mathbf{k})).$$

$$a_1=a(\sqrt{3},0),a_2=a(-\frac{\sqrt{3}}{2},\frac{3}{2})\iff G_1=\frac{4\pi }{3a}(\frac{\sqrt{3}}{2},\frac{1}{2}),G_2=\frac{4\pi }{3a}(0,1),G_3=G_2-G_1$$

$$G_1^M=\frac{8\pi \sin \frac{\theta }{2}}{3a}\text{quantity}(\frac{1}{2},\frac{\sqrt{3}}{2}),G_2^M=\frac{8\pi \sin \frac{\theta }{2}}{3a}\text{quantity}(-\frac{1}{2},\frac{\sqrt{3}}{2}),a_1^M=\frac{\sqrt{3}a}{2\sin \frac{\theta }{2}}\text{quantity}(\frac{\sqrt{3}}{2},\frac{1}{2}),a_2^M=\frac{\sqrt{3}a}{2\sin \frac{\theta }{2}}\text{quantity}(-s\frac{\sqrt{3}}{2},\frac{1}{2})$$

Full text

I Introduction

The recent discovery of correlated insulating states and superconductivity in twisted bilayer graphene (TBG) PabloMott ; PabloSC ; Dean-Young ; Efetov 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 MacDonald ; Santos , several details, including the symmetry breaking within the insulating phase and the nature and mechanism of pairing in the neighboring superconductor, remain under debate Balents18 ; Po2018 ; IsobeFu ; Thomson18 ; YouAV ; Vafek ; Xie2018 ; Nandkishore ; Kivelson ; PhilipPhillips ; phononMacDonald ; phononLianBernevig ; Zou2018 . 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 8 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 TDBGexp2019 , 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 to 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 Koshino2017 ; Koshino2018 , 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 TDBGexp2019 , 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 PabloMott ; Balents18 ; Vafek . 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 Goerbig . 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 non-zero valley Chern number, ferromagnetism may be preferred Zhang2018 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 Crommie2010 . At partial fillings that are integers, ferromagnetic ground states were obtained with repulsive interactions Ghaemi2012 , 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 Sharpe2019 ; Xie2018 ; Zou2018 ; Zhang2018 ; Zhang2019 ; Bultinck19 .

Third, we investigate the nature of the superconducting phase by highlighting that the valley degree of freedom, which behaves as a pseudo-spin, 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 $T_{c}$ on in-plane field TDBGexp2019 .

II Results

II.1 Single-particle physics

We consider a system consisting of two AB-stacked graphene bilayers twisted relative to ABAB stacking by a small angle $\theta$, illustrated in Fig. 1. For a detailed discussion on the Hamiltonian and model parameters, see Methods. 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})$ Koshino2017 ; Koshino2018 . In the original Bistritzer-Macdonald model, $w_{0}$ and $w_{1}$ are taken to be equal MacDonald2011 . However, in a realistic twisted model, the ratio $r\equiv 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 around 0.75 for the first magic angle Koshino2017 ; Koshino2018 . Here, we similarly include lattice relaxations by taking $r$ to be smaller than 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 $\nu=\pm 4$ filling. In this work, we take $(w_{0},w_{1})=(88,100)\textrm{ 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 Tarnopolsky ; Khalaf2019 . In the quadratic approximation of the bilayer-graphene dispersion, the first conduction and valence bands in TDBG become almost perfectly flat at the angle $\theta\approx 1.05$ Zhang2018 . However, once trigonal warping ($\gamma_{3}$) and particle-hole asymmetry ($\gamma_{4}$) terms are included, the flat-bands acquire a significant dispersion and become overlapped with each other ( Fig. 2 a,b). Theses bands can only be separated by applying a strong enough gate voltage between top and bottom layers (Fig. 2 c). Using numerical simulations, we identify the parameter space of twist angle $\theta$ and applied voltage $U$ where the first conduction band is isolated (Fig. 3 a). On the other hand, we find that there is barely any regime where the first valence band is isolated (Fig. 3 c). Such a particle-hole asymmetry in the band structure is originated from $\gamma_{4}$ and $\Delta$ terms. The results are consistent with the experimental findings TDBGexp2019 , 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 to 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-40 meV ChoiSTM2019 .

Another crucial difference compared to TBG is the absence of two-fold 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$) three-fold rotation symmetry $C_{3}$, ($ii$) time-reversal symmetry $\mathcal{T}$, and, ($iii$) mirror reflection about the $x$-axis $M_{y}$ 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 $U(1)_{\textrm{v}}$, arising from the decoupling of Moiré and atomic lattice scale physics.

In addition, the conduction band within each valley carries a non-zero Chern number. In ordinary condensed matter systems, ${\cal 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 ${\cal T}$-symmetry which exchanges the two valleys. Therefore, spontaneous valley polarization would lead to a Chern band without explicitly breaking ${\cal T}$-symmetry Zhang2018 ; Ghaemi2012 ; Xie2018 ; Bultinck19 . At $U=0$, the reflection symmetry $M_{y}$ enforces $C=0$ for both valleys. At $U\neq 0$, the conduction band develops a non-vanishing Chern number computed numerically in Fig. 3 c which is equal to $\pm 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 $\bm{K}_{+}$ and $\bm{K}_{-}$ valleys, which amounts to a Chern number $C_{v}=\pm 2$ for each valley. BLG_valley1 ; BLG_valley2 ; BLG_valley3 ; BLG_valley_exp .

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\mu_{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 Koshino2011 ; McEuen2017 . It leads to an energy correction of $\mu_{B}g^{z}_{\tau}({\bm{k}})B_{z}$, with a ${\bm{k}}$-dependent $g$-factor $g^{z}_{\tau}({\bm{k}})$ satisfying $g^{z}_{-\tau}(-{\bm{k}})=-g^{z}_{\tau}({\bm{k}})$ due to time-reversal symmetry ($\tau$ is a valley index). As shown in Fig. 3 f, $g^{z}_{\tau}({\bm{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 ${\bm{B}}$, we can choose the gauge ${\bm{A}}({\bm{z}})=-{\bm{z}}\times{\bm{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 $-\frac{e}{\hbar}\frac{(l+m)d}{2}\,{\bm{e}}_{z}\times{\bm{B}}$ to the hopping connecting layers from $l$ to $m$, where $d$ is the interlayer separation (See Methods). This leads to an energy correction of the form $\mu_{B}(g^{x}_{\tau}({\bm{k}})B_{x}+g^{y}_{\tau}({\bm{k}})B_{y})$ to the leading order in ${\bm{B}}$ with $g^{x,y}_{-\tau}(-{\bm{k}})=-g^{x,y}_{\tau}({\bm{k}})$. The orbital effect due to in-plane field amounts to a very small relative momentum shift $\sim\frac{eda}{\hbar}\approx 10^{-5}$. However, it cannot be neglected since it is of the same order of magnitude as the Zeeman effect, $\frac{ev_{F}d}{\mu_{B}}\sim 1$ (see Fig. 3 d,e). In general, the in-plane orbital contribution changes the band dispersion due to its ${\bm{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 in-plane field on the insulating gap and the superconducting temperature (See Methods and Supplementary Material 6).

II.2 Correlated insulating states

In the band isolation regime, the first conduction band carries a non-zero Chern number as shown in Fig. 3 a,b which prevents the existence of exponentially localized Wannier functions Vanderbilt1997 . 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 TDBGexp2019 , we perform a self-consistent Hartree-Fock mean field theory similar to the one employed in Ref. Po2018 ; Zhang2018 . Below, we sketch the derivation from the microscopic theory, relegating most details to Supplementary Material 2 and 3. The interacting Hamiltonian in momentum space is given by

[TABLE]

where $V({\bm{q}})$ is the Fourier-transformed screened Coulomb interaction Young_FCI_2018 ; Screening_2 . Since the screening coming from the distance between the system and the gate is comparable to the Moiré length scale, the screening length can be important for the interaction effects. The density $\hat{\rho}({\bm{q}})$ consists of an intravalley part $\rho^{+}\sim c_{\pm}^{\dagger}c_{\pm}$ and an intervalley part $\rho^{-}\sim c_{\pm}^{\dagger}c_{\mp}$, where $c^{\dagger}_{\pm}$ is the electron creation operator for ${\bm{K}}_{\pm}$-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, $\mathcal{H}_{\rm int}=\mathcal{H}_{0}+\mathcal{H}_{J}$, where $\mathcal{H}_{0}$ contains the coupling between intravalley densities $\rho^{+}\rho^{+}$ whereas $\mathcal{H}_{J}$ contains the coupling between intervalley densities $\rho^{-}\rho^{-}$. Rough estimation for the relative energy scales for $H_{0}$ and $H_{J}$ gives ${V_{0}}\sim 35\textrm{ meV}$ and $J\sim 0.6\textrm{ 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 inter-valley Hund’s coupling as in Zhang2018 .

Within the self-consistent Hartree-Fock mean field theory, we consider the order parameter defined as

[TABLE]

For a gapped phase, matrix $M({\bm{k}})$ must be a projector, i.e. $M({\bm{k}})^{2}=M({\bm{k}})$ satisfying $\mathop{\mathrm{tr}}M({\bm{k}})=\nu$ for all ${\bm{k}}$. Given that there are four flavors of fermions due to spin ($\sigma$) and valley ($\tau$) degeneracies, any possible order parameter $M$ can be expanded in terms of the generators of SU(4) $\sigma_{i}\otimes\tau_{j}$, which can be grouped based on their symmetry breaking into 5 categories: (i) $\{\sigma_{0}\tau_{z}\}$ only breaks $\mathcal{T}$ and corresponds to a valley-polarized (VP) state, (ii) $\{\sigma_{x,y,z}\tau_{0}\}$ breaks spin-rotation symmetry and correspond to a spin-polarized (SP) state. (iii) $\{\sigma_{x,y,z}\tau_{z}\}$ breaks both spin-rotation and time-reversal (but preserve some combination of the two) and corresponds to a spin-valley locked (SVL) state, (iv) $\{\sigma_{0}\tau_{x,y}\}$ breaks $U(1)$ valley-charge conservation and corresponds to an inter-valley coherent (IVC) state, and (v) $\{\sigma_{x,y,z}\tau_{x,y}\}$ breaks both spin-rotation and U(1)v valley-charge conservation, corresponds to spin-IVC locked (SIVCL) state (see Table 2). We note that any of these orders may break or preserve $C_{3}$ symmetry depending on its ${\bm{k}}$ dependence.

The results of the self-consistent Hartree-Fock calculation are summarized in the following (See Supplementary Material 3 for details). Restricting ourselves to translation-symmetric gapped states, we find there are five options: SP, VP, SVL, IVC, and SIVCL at half-filling $\nu=2$ and three options: spin-valley-polarized (SVP), spin-polarized-IVC (SPIVC), and spin-valley-locked-IVC (SVLIVC) at quarter-filling $\nu=1,3$, as in Tab. 2. By solving the Hartree-Fock self-consistency condition, the ground state energy $E$ and the correlation gap $\Delta$ are computed for different states, Fig. 4 a. 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 (see supplemental material). 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. 4 a 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 $\Delta$ for the SP state in the band isolation region ranges between 4 and 8 meV, see Fig. 4 b.

To understand the reason why IVC order is energetically unfavorable, we can employ the argument of Ref. Bultinck19 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 $|C|$ 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 supplemental material.

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 $\Delta_{\rm SP}$. On the other hand, it reduces the energy of some of the empty bands for the VP-state, reducing $\Delta_{\rm VP}$ (see Fig. 4 c,e). The Hund’s coupling term similarly reduces $\Delta_{\rm SVP}$ at quarter filling by lowering the energy of one of the excited states (see Fig. 4 d). 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 to the latter and requires the application of a magnetic field TDBGexp2019 .

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. 4 e), which is in agreement with the experimental data TDBGexp2019 . From the numerical calculation, we confirmed that such a reduction in gap also depends on the in-plane field direction, which exhibits 3-fold periodicity (see Methods). 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\pi/3$ in the angle.

II.3 Superconductivity

When the correlated insulator is doped away from half-filling, a superconducting phase is observed below 3.5 K TDBGexp2019 . Our proposed scenario for the observed superconductivity is illustrated in Fig. 5 $\mathbf{a}$ where pairing takes place between time-reversal partners in opposite valley. Such an intervalley pairing between time-reversal partners has also been proposed SC_TMD_theory1 ; SC_TMD_theory2 ; SC_TMD_theory3 and observed in transition metal dichalcogenides (TMD) SC_TMD_exp . 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 saxena_superconductivity_2000 ; aoki_coexistence_2001 ; huy_superconductivity_2007 , ferromagnetic spin fluctuations can act as a pairing glue responsible for superconductivity TBG_SC_triplet2018 . This motivates the following simplified Hamiltonian,

[TABLE]

where the spin operator ${\bm{S}}^{a}_{\bm{q}}=\sum_{{\bm{k}},\tau,\sigma,\sigma^{\prime}}c^{\dagger}_{\sigma,\tau,{\bm{k}}+{\bm{q}}}{\bm{\sigma}}^{a}_{\sigma,\sigma^{\prime}}c_{\sigma^{\prime},\tau,{\bm{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 ${\bm{q}}=0$ which justifies a ${\bm{k}}$-independent coupling.

Next, we consider the simplest possible intervalley superconducting pairing function $\Delta$ which is ${\bm{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, $\Delta$ is proportional to $\tau_{x}$ or $\tau_{y}$ corresponding to valley triplet or singlet, respectively. The overall antisymmetry of $\Delta$ implies that the former scenario corresponds to a spin-singlet $i\sigma_{y}$ whereas the latter corresponds to a spin-triplet $i\sigma_{y}{\bm{d}}\cdot{\bm{\sigma}}$. Here, ${\bm{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

[TABLE]

We now assume ${\bm{k}}$-independent $\Delta$ and decompose it into spin-singlet/velly triplet $\Delta_{s}$ and spin-triplet/valley singlet $\Delta_{t}$. We now use

[TABLE]

where $\lambda_{t}=1$ and $\lambda_{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 detailed discussion of these pairing channels within the linearized BCS equation is provided in Supplementary Material.

We highlight here that spin-triplet pairing is only known to occur in liquid He3 Ambegaokar73 and a few Uranium compounds saxena_superconductivity_2000 ; aoki_coexistence_2001 ; huy_superconductivity_2007 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 ${\bm{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 $\Delta=\tau_{y}\sigma_{y}{\bm{d}}\cdot{\bm{\sigma}}$ in the presence of a magnetic field ${\bm{B}}$. Restricting ourselves to terms up to quartic order in ${\bm{d}}$ or ${\bm{B}}$, we can write the following free energy functional

[TABLE]

Detailed microscopic derivation of the coefficients $a,b,c,\kappa,\alpha,\eta$ is provided in supplemental material. 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, ${\bm{B}}=B{\bm{e}}_{z}$, we can then write

[TABLE]

Substituting in the free energy (6) and using the fact that $\eta=-\alpha/2$ (see supplemental material) yields

[TABLE]

One important feature is that $\alpha>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

[TABLE]

The most remarkable feature of this result is that, for non-zero $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

[TABLE]

where $\Lambda$ is the bandwidth, $N(0)$ is the density of states at the Fermi energy, and $\chi$ is the dimensionless magnetic susceptibility (see supplemental material). Similar linear field-dependence of $T_{c}$ is known in superfluid He3 Ambegaokar73 , 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 spin-singlet superconductor. One crucial observation here is that $a$ seems to depend on several details and is expected to be very small since $T_{c}\ll\frac{N(0)}{N^{\prime}(0)}\sim\epsilon_{F}$. Surprisingly, the measured value of $a$ is of order 1 TDBGexp2019 which suggests the vicinity of a quantum critical point where the scaling of the susceptibility cancels exactly against the other parameters. Indeed, the scaling $\chi\sim\epsilon_{F}/(T\log T)$ predicted by Herz-Millis theory in the quantum critical regime for an itinerant ferromagnet Herz76 ; Millis93 leads to such cancellation resulting in $a\sim 1$.

The origin of the quadratic term in Eq. 9 can be understood in terms of the in-plane orbital effect discussed in Sec. II.1. First, note that Zeeman splitting cannot break Cooper pairs between aligned spins. Instead, it yields an initial linear increase in $T_{c}({\bm{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 Material 6)

[TABLE]

where ${\bm{e}}_{{\bm{B}}}$ is the direction of the external magnetic field. The average value of $({\bm{e}}_{\bm{B}}\cdot{\bm{g}}_{+}({\bm{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 superconductivity as $\mu_{B}B_{c}\sim\sqrt{T_{c}/b}$ yielding a value about 3 Teslas which compares favorably to the experimental value TDBGexp2019 . Furthermore, if we consider an out-of-plane field instead, $|g_{z}|$ is on average about 1-2 orders of magnitude larger than $|g_{x,y}|$, yielding a critical field of about $\sim 0.1T$ which is very close to the experimentally observed result TDBGexp2019 .

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 Hattori11 . Such effects are neglected within our simplified analysis (S74) which assumes a constant coupling $g$.

III Conclusion

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. TDBGexp2019 .

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 non-zero 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 to the bandwidth, this effect is comparable to Zeeman splitting, leading to a modified $g$-factor which compares favorably to the experimental value TDBGexp2019 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 TDBGexp2019 . 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.

Note: After completing this work we noticed two experimental papers PabloTDBG ; IOP_TDBG which are consistent with Ref. TDBGexp2019 and theoretical discussion contained here.

IV 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:

[TABLE]

which is labelled in the order of $A_{\textrm{1}}$, $B_{\textrm{1}}$, $A_{\textrm{2}}$, $B_{\textrm{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 $\Delta$ for these sites. Here, $f({\bm{k}})\equiv\sum_{l}e^{-i{\bm{k}}\cdot\delta_{l}}$, where $\delta_{1}=a(0,-1)$, $\delta_{2}=a(-\sqrt{3}/2,1/2)$, and $\delta_{3}=a(\sqrt{3}/2,1/2)$ are vectors from $B$-site to $A$-sites. One can expand $f({\bm{k}})$ near ${\bm{K}}_{\pm}=\pm(4\pi/3\sqrt{3}a,0)$ as

[TABLE]

where $a$ is the distance between carbon atoms. Throughout, we will use the phenomenological parameters extracted from Ref. Jeil2014

[TABLE]

where $\gamma_{0,1,3,4}$ and $\Delta$ are the parameters illustrated in Fig. 1. Additionally, 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$, ABAB system’s dielectric constant $\epsilon$ and the thickness of the BLG/BLG system $d$, $U={\epsilon^{-1}}D\cdot 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 ${\bm{K}}_{\pm}$ points. In the continuum model approximation MacDonald2011 , Moire bands from ${\bm{K}}_{\pm}$ valleys decouple; for the Moire band from ${\bm{K}}_{+}$ valley, the Hamiltonian is given by

[TABLE]

where $c^{\dagger}_{{\bm{k}},+,t/b}$ is a 4-components electron creation operator for top/bottom layer with momentum ${\bm{K}}_{+}+{\bm{k}}$. Here, $h_{\theta}({\bm{k}})=h(R_{-\theta}{\bm{k}})$ with $R_{\theta}$ denoting the counter-clockwise rotation matrix by angle $\theta$ relative to the $x$-axis. The momenta ${\bm{q}}_{0,1,2}$ are given by ${\bm{q}}_{0}=R_{\theta/2}K-R_{-\theta/2}K=\frac{8\pi\sin(\theta/2)}{3\sqrt{3}a}(0,-1)$, ${\bm{q}}_{1}=R_{\phi}{\bm{q}}_{0}$, and ${\bm{q}}_{2}=R_{-\phi}{\bm{q}}_{0}$ where $\phi=2\pi/3$. The hopping matrices $T_{n}$, $n=0,1,2$ are given by

[TABLE]

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\in[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. 7 a,b,c expand in the parameter space of $(\theta,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 ${\bm{K}}_{\pm}$-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}\mapsto-k_{y}$ so Berry curvature flips its sign Zhang2018 . 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 $\pm 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. 3 b. 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 $\pm 3$. Then, the band inversion would occur at the center Dirac cone, exchanging Chern number by $\mp 1$. In total, it will still change the Chern number by $\pm 2$. At larger values of the gate voltage $U$, the band inversion happens between first and second conduction band at $\Gamma$ point, and the Chern number then changes by $\mp 1$ (It can change by $\mp 2$ for other parameter setting), decreasing the Chern number.

This can be further checked by inspecting symmetry indicators Bernevig2012 ; Turner2012 ; Watanabe2017 . There are three $C_{3}$-invariant momenta $\Gamma$, $K$, and $K^{\prime}$. For a Bloch state with these momenta, $C_{3}$ rotation symmetry would map the state back to itself with a rotation eigenvalue:

[TABLE]

where $L_{n,{\bm{k}}}$ is an angular momentum associated with the Bloch state $\ket{{\bm{k}},n}$. Then, the Chern number of the $n$-th band can be determined modulo 3 by

[TABLE]

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 ${\bm{K}}_{+}$ valley at $\theta=1.33^{\circ}$. At $U=0$ meV, we start with $(n_{\Gamma},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\mapsto 1$. At $U=90$ meV, Chern number again changes by $-1$, manifested by $n_{\Gamma}:0\mapsto-1$. See Fig. 6 for the detail.

Magnetic Field Effect Under in-plane magnetic field ${\bm{B}}=(B_{x},B_{y},0)$, one can choose the gauge ${\bm{A}}({\bm{z}})=-{\bm{z}}\times{\bm{B}}$. Then, the effect of a magnetic field on hopping terms is evaluated via Peierl’s substitution, where the hopping term from ${\bm{R}}$ to ${\bm{R}}+{\bm{\delta}}$ is multiplied by the phase factor

[TABLE]

such that

[TABLE]

where $\bm{\alpha}=-\frac{q}{\hbar}{\bm{A}}({\bm{R}}_{z}+{\bm{\delta}}_{z}/2)=-\frac{e}{\hbar}\left[\quantity({\bm{R}}_{z}+\frac{{\bm{\delta}}_{z}}{2})\times{\bm{B}}\,\right]$ since ${\bm{A}}({\bm{z}})$ is linear function of ${\bm{z}}$. Hence, the effect of in-plane field can be included by simply replacing all ${\bm{k}}$-dependent matrix elements of Bloch Hamiltonians by ${\bm{k}}+\bm{\alpha}$ as follows (we take $c_{\bm{k}}=\sum_{\bm{R}}e^{-i{\bm{k}}\cdot{\bm{R}}}c_{\bm{R}}$):

[TABLE]

where $\mathcal{H}_{l,m}$ is the matrix element connecting layers $l$ and $m$ ($l,m=0,\dots,3$ from bottom to top) in Eq. IV, $d=3.42\textrm{\AA}$ is the interlayer distance, and ${\bm{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

[TABLE]

where ${\bm{g}}^{xy}_{n,\tau}({\bm{k}})$ is given by

[TABLE]

where $\tau$ is the valley index. Time-reversal symmetry implies that ${\bm{g}}^{xy}_{n,\tau}(-{\bm{k}})=-{\bm{g}}^{xy}_{n,-\tau}({\bm{k}})$. The in-plane orbital $g$-factor transforms under $C_{3}$ rotation as

[TABLE]

provided that the band $n$ is non-degenerate at ${\bm{k}}$. This implies that ${\bm{g}}^{xy}_{n,\tau}({\bm{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 Material 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. Koshino2011 ; McEuen2017 whose $g$-factor is given by

[TABLE]

In summary, the single particle energies has the following dependence on magnetic field

[TABLE]

where ${\bm{\sigma}}$ is the electron spin operator (which is $\pm 1/2$ for up/down spins) and $\tau=\pm$. The valley orbital $g$-factor is defined as

[TABLE]

We have also assumed that the spin-quantization axis is parallel to the field.

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

Acknowledgements.

We thank Shiang Fang, Yahui Zhang, Yizhuang You, Erez Berg, and Bertrand Halperin for helpful discussion. In particular, we thank Mikito Koshino for clarification on his earlier works on BLG parameters. A. Vishwanath, J.Y. Lee and E. Khalaf were supported by a Simons Investigator Fellowship. P. Kim, X. Liu and Z. Hao acknowledge partial support from the Gordon and Betty Moore Foundation’s EPiQS Initiative through Grant GBMF4543 and the DoD Vannevar Bush Faculty Fellowship N00014-18-1-2877.

I Details of Twisted Double Bilayer Graphene Continuum Model

Before getting into the detailed calculation, here we clarify some subtleties in the Moiré continuum approach and tight-binding parameters. As we will see, this is crucial because it can affect the physical band structures of a twisted double bilayer graphene. First, let us fix the lattice convention. Let $a=1.42\textrm{ \AA}$ be the distance between carbon atoms. The original hexagonal lattice is defined by

[TABLE]

Here, ${\bm{K}}_{\pm}$-point is given by $\frac{4\pi}{3\sqrt{3}a}(\pm 1,0)$. Once we twist two layers, Moiré structure is formed with a spatial modulation given by linear combinations of $\{{\cal R}_{\theta/2}G_{1},{\cal R}_{\theta/2}G_{2},{\cal R}_{-\theta/2}G_{1},{\cal R}_{-\theta/2}G_{2}\}$, where ${\cal R}_{\theta}$ rotates a vector by angle $\theta$ counterclockwise. Periodicity of the given system is then governed by the smallest reciprocal lattice vector that can be obtained by linear combination. These are $({\cal R}_{\theta/2}-{\cal R}_{-\theta/2})G_{1}$ and $({\cal R}_{\theta/2}-{\cal R}_{-\theta/2})G_{2}$, which correspond to Moiré reciprocal lattice vectors. Therefore, we obtain $G_{1,2}^{M}=\frac{8\pi\sin\theta/2}{{3}a}(-1/2,\pm\sqrt{3}/2)$. For convenience, use the following Moiré reciprocal lattice vectors from now on:

[TABLE]

Before getting into detail, let us fix the Fourier transform convention by $c^{\dagger}_{\bm{k}}=\sum_{\bm{R}}e^{i{\bm{k}}\cdot{\bm{R}}}c^{\dagger}_{\bm{R}}$. Here, $c^{\dagger}_{\bm{R}}$ creates a Wannier orbital $W({\bm{r}}-{\bm{R}})$ centered at ${\bm{R}}$. This is consistent with the other convention used throughout the paper, where $\ket{\psi_{\bm{k}}}=\sum_{R}e^{i{\bm{k}}\cdot{\bm{R}}}\ket{{\bm{R}}}$. Under this choice, whenever there is a hopping from ${\bm{R}}$ to ${\bm{R}}^{\prime}$ in real space, one obtains the term proportional to $e^{-i{\bm{k}}\cdot({\bm{R}}^{\prime}-{\bm{R}})}$ in the Bloch Hamiltonian. Here, the tight-binding model for bilayer graphene can be fully characterized by the parameters $(\gamma_{0},\gamma_{1},\gamma_{3},\gamma_{4},\Delta)$, where the hopping term for nearest neighbor $\gamma_{0}$ is intentionally taken to have an additional minus sign from the hopping integral so that all $\gamma_{i}$’s are positive. The sign difference between vertical ($V_{pp\pi}$) and horizontal hopping ($V_{pp\sigma}$) overlap integral is originated from the phase structure of $2p_{z}$ orbital. (See how Slater-Koster parameters SlaterKoster1954 are calculated) In Ref. McCann2013 , the sign of the trigonal warping $\gamma_{3}$ is taken to be negative, which is an inaccurate choice of the parameter because different sign convention would flip the shape of trigonal warping. Following the DFT result in Ref. Charlier1991 , positive sign in front of $\gamma_{3}$ should be a proper choice for realistic materials.

Under this convention, where the phase structure for $2p_{z}$ orbitals at every carbon site is taken to be equivalent, one can derive the Moiré hopping term as the following. Once we have a Moiré structure, between momentum points of top and bottom layers, there exists Moiré-hopping term whose momentum transfer is given by linear combinations of $G^{M}_{1}$ and $G^{M}_{2}$. As we are interested in a Moire band structure near charge neutral point, we only consider electron momenta near Dirac points for top and bottom layers, ${\cal R}_{\pm\theta/2}{\bm{K}}_{\pm}$-points. For example, with respect to ${\cal R}_{\pm\theta/2}{\bm{K}}_{+}$, the momentum transfer condition can be written as

[TABLE]

where $q_{1}=({\cal R}_{-\theta/2}-{\cal R}_{\theta/2}){\bm{K}}_{+}=K_{M}(0,-1)$ with ${K_{M}}=\frac{8\pi\sin\frac{\theta}{2}}{3\sqrt{3}a}$. Let us denote ${\bm{K}}_{+}^{t}=\quantity({\cal R}_{\theta/2}{\bm{K}}_{+})$ and ${\bm{K}}_{+}^{b}=\quantity({\cal R}_{-\theta/2}{\bm{K}}_{+})$. The Moiré state with a momentum $p$ is given by a superposition of Bloch states of top and bottom layers with (absolute) momenta $\{p+n_{1}G_{1}^{M}+n_{2}G_{2}^{M}|n_{1},n_{2}\in\mathbb{Z}\}$. Since we want to solve the equation in terms of Dirac Hamiltonians with respects to ${\bm{K}}_{+}^{b,t}$-points, we do the following precedure. Define $\tilde{p}=p-{\bm{K}}_{+}^{t}$. Then, the Moiré state with momenta $p$ is composed of the Bloch states with the following momenta defined with respect to ${\bm{K}}_{+}^{t}$ and ${\bm{K}}_{+}^{b}$:

[TABLE]

where $q_{2}=(G_{1}^{M}+q_{1})=K_{M}(\sqrt{3}/2,1/2)$ and $q_{3}=(G_{2}^{M}+q_{1})=K_{M}(-\sqrt{3}/2,1/2)$. Therefore, the state with a momentum $p=\tilde{p}+{\bm{K}}_{+}^{t}$ can be solved by considering coupling among bloch states from $H({\cal R}_{-\theta/2}[\tilde{p}-K_{\textrm{lat}}^{t}])$ and $H({\cal R}_{\theta/2}[\tilde{p}-K_{\textrm{lat}}^{b}])$. Here, $K_{\textrm{lat}}^{t,b}$ is a vector denoting the location of each lattice point in $\{n_{1}G_{1}^{M}+n_{2}G_{2}^{M}:n_{1,2}\in\mathbb{Z}\}$ for the top layer and $\{q_{1}+n_{1}G_{1}^{M}+n_{2}G_{2}^{M}:n_{1,2}\in\mathbb{Z}\}$ for the bottom layer in the $k$-space.

Let top and bottom layer have relative shift $\delta$. For convenience, take a frame where the bottom layer is fixed and the top layer is rotated by $\theta$. Then, Each original Bloch state is represented as

[TABLE]

where $\alpha,\beta$ are sublattice indices and $\tau_{\alpha}$ is associated displacement. Here, primed coordinate $R^{\prime}={\cal R}_{\theta}(R+\delta)$ and $\tau^{\prime}={\cal R}_{\theta}\tau$ are for the top layer, meaning that it is rotated by $\theta$ along counter-clockwise direction. Here, $R=n_{1}a_{1}+n_{2}a_{2}$. The initial displacement $\delta$ between layers is not important in the end, as we will see. Let $\tau_{0}=a(0,1)$. In the mono-mono case $\delta=0$ for $A-A$ stacking and $\delta=\tau_{0}$ for $A-B$ stacking. $A-B$ stacking means that $A$-site of the top layer is placed at the location of $B$-site of the bottom layer.

By definition, $A$-site is original site spanning the lattice, and $B$-site is displaced by $\tau_{0}$ with respect to the $A$-site. In the bilayer(AB)-bilayer(AB) case, $\delta=-\tau_{0}$ because between layers it is $BA$ stacking. For the following calculation, we take $\delta=-\tau_{0}$ and $d=0$. Now, hopping matrix element from top second layer to bottom first layer $H_{b,t}$ can be evaluated by

[TABLE]

where we used $t(R)=\frac{1}{N}\sum_{q}e^{iq\cdot R}t_{q}\approx\frac{1}{\Omega}\int d^{2}q\,e^{iq\cdot R}t_{q}$, and ${\cal G}$ and ${\cal G}^{\prime}$ are sets of reciprocal lattice vectors for original bottom and top lattices, respectively. In principle, for a different gauge choice $t({\bm{R}})$ depends on $\alpha,\beta$ indices as well, but since we took the gauge choice where all phase structures for $2p_{z}$ orbitals are the same, the dependence will be trivial. In the last line, we just did change of variables. Now we can see that it is nonvanishing only when

[TABLE]

Considering that $t_{q}$ is decreasing fast with $q$, we can only retain most relevant terms where $k+g_{1}$ is minimized. In terms of a momentum relative to the ${\bm{K}}^{t,b}_{+}$ points, we have

[TABLE]

Naively, when we consider $k+g_{1}$, since $k$ does not deviate much from ${\bm{K}}_{+}^{b}$-point, the most relevant $t_{k+g_{1}}$ would be given when $g_{1}=g_{2}=0,G_{3},-G_{1}$ so that $\absolutevalue{{\bm{K}}_{+}^{b}+g_{1}}=\absolutevalue{{\bm{K}}_{+}^{b}}$. For cases with $g_{1}\neq g_{2}$, $(i)$ the energy difference between top and bottom electrons are very large, and $(ii)$ $t_{k+g_{1}}$ would be small, and therefore we ignore such cases. In fact, each of these cases corresponds to when $\tilde{p}-\tilde{k}=q_{1},q_{2},q_{3}$:

[TABLE]

With this understanding, one can write down three hopping matrices as the following:

[TABLE]

where $z=e^{2\pi i/3}$ and since $G_{3}\cdot\tau_{0}=G_{1}\cdot\tau_{0}=2\pi/3$. Due to the fact that $\delta=-\tau_{0}$ instead of $\tau_{0}$, the form is slightly different from the TBG case MacDonald2011 ; Po2018 . By proper phase redefinition of Bloch states represented by momentum lattices (gauge degrees of freedom for Bloch states), we can absorb $z$ and $z^{*}$ factors in front of matrices. Thus, the form of the hopping matrices can be simplified. Moreover, an initial displacement between two layers is not important. In this derivation, it is not difficult to notice that diagonal and off-diagonal entries for Moiré hopping matrices can be different. If there is an additional spatial modulation with a Moiré scale, given differently between AA(BB) and AB(BA) sites in $t(R+\tau_{\alpha}-R^{\prime}-\tau_{\beta})$, one would obtain a different values for $t_{k}$, as explained in Ref. Moon2013 ; Koshino2018 . Finally, to obtain an energy spectrum at Moiré momentum $k$, one needs to diagonalize the following Hamiltonian with a certain cutoff:

[TABLE]

Now, we want to point out some subtlety for the generic Moiré hopping matrix:

[TABLE]

where the Pauli operator $\sigma_{i}$ acts on the sublattice basis. This is the form of the hopping term written in Ref. MacDonald2011 , where the ${\bm{K}}$ and ${\bm{K}}^{\prime}$ Dirac Hamiltonian was written as111$H({\bm{K}}+{\bm{k}})$ in the MacDonald’s paper is a Hamiltonian with respect to the ${\bm{K}}^{\prime}$-point in our convention.

[TABLE]

Now, imagine we choose a different basis choice, for example multiplying $(-)$ sign for the $B$-sublattices. This is equivalent to apply $\sigma_{3}$ transformation to the operators, and as a result, both $H({\bm{K}}+{\bm{k}})$ and $H({\bm{K}}^{\prime}+{\bm{k}}^{\prime})$ would change its sign. This is the basis chosen in Ref. Po2018 . Accordingly, interlayer hopping term $T_{0}$ would change as well, from $w_{0}+w_{1}\sigma_{1}$ to $w_{0}-w_{1}\sigma_{1}$ and similarly for others.

II Interaction Projection and Intervalley Hund’s Coupling

Here, we provide the details for the procedure of projecting the Coulomb interaction on the isolated flat band and how to derive the intervalley Hund’s coupling. The interaction Hamiltonian can be written as

[TABLE]

Here, ${\bm{r}}$ integration is over the whole space not just the unit cell. For the screed Coulomb interaction, $V({\bm{r}})$ is given by

[TABLE]

where $\kappa$ denotes the inverse screening length. The density operator is given by

[TABLE]

where $\sigma$, $\sigma^{\prime}$ sum over spin states $\uparrow$, $\downarrow$ and $\tau$, $\tau^{\prime}$ sum over vallyes $\pm$, and $n$ sums over the relevant set of bands. In the following, we will restrict ourselves to the isolated Moiré band and drop the band index $n$. Expansion in the Bloch basis is done by writing

[TABLE]

where $N$ is the number of momentum point in the first Brillouin zone which equals to the total number of Moiré unit cells in system, and $\psi_{\bm{k}}({\bm{r}})$ are the Bloch states satisfying $\psi_{\bm{k}}({\bm{r}}+{\bm{R}})=e^{i{\bm{k}}\cdot{\bm{R}}}\psi_{\bm{k}}({\bm{r}})$ for a given Moiré lattice translation ${\bm{R}}$. We now split the density into intra- and intervalley components

[TABLE]

Here, we used the fact that different spin states are orthogonal. If valley symmetry is exact, states belonging to different valleys would also be orthogonal leading to a vanishing intervalley density $\rho^{-}_{\sigma,\tau}$. However, valley symmetry is broken on the scale of $~{}|K-K^{\prime}|^{-1}$ leading to a very small intervalley Hund’s coupling term. This term can be usually neglected since it is much smaller than the interaction between intravalley densities. Nevertheless, contributions from this term can lift the degeneracy between different broken symmetry states which are otherwise exactly degenerate, which makes it important to include it in our analysis. We note that the Bloch states are generally vectors with some internal index denoting layer, sublattice, etc which means that the combination $\psi^{\dagger}\psi$ above denotes an inner product in these internal indices.

The Bloch states can be written in terms of the periodic function $u_{\bm{k}}({\bm{r}})$ which can be expanded in a Fourier series in reciprocal lattice vectors ${\bm{G}}$ leading to

[TABLE]

where ${\bm{G}}$ is the Moiré reciprocal lattice vector, and $\Omega$ is the area of the Moiré unit cell. Here, $u_{{\bm{k}}}({\bm{G}})$ are normalized such that $\sum_{\bm{G}}u^{\dagger}_{{\bm{k}}}({\bm{G}})u_{{\bm{k}}}({\bm{G}})=1$. In addition, we can choose the gauge such that the Bloch states satisfy

[TABLE]

If the band has a non-vanishing Chern number, it is impossible to choose a smooth and periodic gauge and there would be an additional phase factor in front of the RHS Kohmoto1985 . In this case, the condition (Eq. S21) implies a discontinuity of the phase of $u_{\bm{k}}$ at the Brillouin zone boundaries.

The interacting Hamiltonian in momentum space is given by

[TABLE]

where $V({\bm{q}})=\int d{\bm{r}}V({\bm{r}})e^{-i{\bm{q}}{\bm{r}}}$ and $\textrm{Vol}=N\Omega$. We note that the Fourier transform of $\rho^{\pm}({\bm{r}})$ is not restricted to momenta inside the Moiré BZ and it should be expressed in terms of a general momentum ${\bm{q}}$. The density $\hat{\rho}({\bm{q}})$ is generally non-periodic in ${\bm{q}}$ under reciprocal Moiré lattice translations since the Bloch states have a non-trivial spatial structure inside the Moiré unit cell. Instead, it decays over some momentum scale comparable to the Moiré Brillouin zone size. On the other hand, the Bloch states has no structure inside the unit cell of the original bilayer graphene where a tight-binding description of the orbitals was employed. Hence, the density $\hat{\rho}({\bm{q}})$ is periodic under any reciprocal lattice translation for the original system. As a result, $\hat{\rho}({\bm{q}})$ consists of several identical narrow peaks centered at reciprocal lattice vectors of the original bilayer graphene $\tilde{\bm{G}}$ for the intravalley density $\rho^{+}$ or at $K-K^{\prime}+\tilde{\bm{G}}$ for the intervalley density $\rho^{-}$. This poses a problem since it implies that the summation over ${\bm{q}}$ in Eq. S22 diverges.

To resolve this issue, we notice that the periodicity of $\hat{\rho}({\bm{q}})$ in reciprocal space for the original lattice is an artifact of the tight-binding approximation, where an atomic orbital is taken to be point-like. If we instead use the actual shape of the Wannier orbital, the density operator $\hat{\rho}({\bm{q}})$ will decay for momenta larger than a certain cutoff $\Lambda$ which is given by the inverse size of the Wannier orbitals. Rather than attempting to precisely determine the value of $\Lambda$ from the graphene Wannier orbitals, we will consider $\Lambda$ as a phenomenological parameter of the same order as the size of the original Brillouin zone. This will have the effect of restricting the sum over momenta in Eq. S22 to the vicinity of ${\bm{q}}=0$ for the intravalley density $\hat{\rho}^{+}({\bm{q}})$ and the vicinity of $K-K^{\prime}$ and $R_{\pm 2\pi/3}(K-K^{\prime})$ for the intervalley density $\hat{\rho}^{-}({\bm{q}})$.

Therefore, we restrict ourselves to the vicinity of 0 for $\rho^{+}({\bm{q}})$ and $K-K^{\prime}$ (and its rotation related points) for $\rho^{-}({\bm{q}})$. In the following, we perform Fourier transform in terms of small deviations around these momenta by defining $\rho_{\sigma,\tau}^{\alpha}({\bm{q}})$ as(note that $c_{{\bm{k}}+{\bm{G}}}=c_{\bm{k}}$):

[TABLE]

Here, we introduced $K_{+}=K$ and $K_{-}=K^{\prime}$ and we used that $\psi_{\tau,{\bm{k}}}({\bm{r}})=\sum_{\bm{G}}e^{i({\bm{k}}+K_{\tau}+{\bm{G}})\cdot{\bm{r}}}u_{\tau,{\bm{k}}}({\bm{G}})$. In addition, we introduced the intra- and intervalley form factors defined by

[TABLE]

The function ${\bm{p}}({\bm{q}})$ and ${\bm{G}}({\bm{q}})$ are defined to give the projection onto the first BZ and the reciprocal lattice vector corresponding to ${\bm{q}}$, respectively, such that ${\bm{q}}={\bm{G}}({\bm{q}})+{\bm{p}}({\bm{q}})$. The last equality is important for the numerical implementation because the summation over ${\bm{q}}$ can go outside the first BZ whereas the numerical calculation is only carried out in the first BZ.

Time-reversal symmetry dictates that

[TABLE]

for some reciprocal lattice vector ${\bm{G}}_{0}$. This relation can be exploited for the evaluation of form factors. In fact, a direct numerical evaluation gives ${\bm{G}}_{0}$ in our setting.

The form factors satisfy the identities

[TABLE]

The first identity follows from the definition of the form factor, the second from time-reversal symmetry (Eq. S25) and the third from our periodic gauge choice (Eq. S21).

Finally, the resulting interaction can be expanded as a sum of four terms: $\rho^{+}\rho^{+}$ containing intravalley densities, $\rho^{-}\rho^{-}$ containing intervalley densities and two cross terms $\rho^{+}\rho^{-}$. The latter ones have to vanish since they necessarily involve densities at large momenta ${\bm{q}}\pm(K_{+}-K_{-})$ (due to the factor $\lambda_{\bm{q}}$ which is assumed to decay with ${\bm{q}}$). The $\rho^{-}\rho^{-}$ terms is only non-vanishing when $\tau=-\tau^{\prime}$. In addition, since ${\bm{q}}$ is much smaller than $|K_{+}-K_{-}|$, we can ignore the ${\bm{q}}$ dependence in the interaction term and replace it by the constant $V(|K_{+}-K_{-}|)$. Thus, the resulting Hamiltonian consists of two parts

[TABLE]

$\mathcal{H}_{0}$ contains the coupling between intravalley densities $\rho^{+}\rho^{+}$ whereas $\mathcal{H}_{J}$ contains the coupling between intervalley densities $\rho^{-}\rho^{-}$. They are given explicitly by

[TABLE]

[TABLE]

where the intravalley and intervalley form factors $\lambda^{\pm}_{\tau,{\bm{q}}}({\bm{k}})$ are defined as

[TABLE]

All momenta in Eq. S28 and Eq. S29 are measured in units of $q_{M}=\frac{4\pi\theta}{3\sqrt{3}a}$ with $v_{\bm{q}}=\absolutevalue{q_{M}}/\sqrt{{\bm{q}}^{2}+\kappa^{2}}$ denoting the dimensionless screened Coulomb interaction with a screening length $1/\kappa$. The main source of screening is from the gate, which has the distance about 30-50 nm from the sample. The distance is comparable to the Moiré length scale, implying that the screening length can be important. In the following calculation, we would use $\kappa=5\times 10^{7}\textrm{ m}^{-1}$. Rough estimations for $V_{0}$ and $J$ provide the scale of the two interaction terms and are given by

[TABLE]

Here, we used $|\Omega|=\frac{3\sqrt{3}a^{2}}{2\theta^{2}}$ and used $\theta^{o}$ to denote the value of $\theta$ in degrees. Using a value of $\epsilon$ of about 5 at twist angles around $1^{o}$ yields $V_{0}\approx 35$ meV and $J=0.6$ meV. We see that the $J$ term is significantly smaller than the $V_{0}$ term. It can be important, however, since it identifies the two separate spin-rotation symmetry for ${\bm{K}}_{\pm}$ valleys $\textrm{SU(2)}_{+}\times\textrm{SU(2)}_{-}$ down to the single spin-rotation SU(2) symmetry, while preserving valley U(1) symmetry. Thus, it can lift the degeneracy between some symmetry breaking states which are degenerate on the level of the $V_{0}$ interaction. The $J$ term generally has the effect of favoring spin alignment and can be written in the form of inter-valley Hund’s coupling as in Zhang2018 .

We notice that the interaction term is invariant under the gauge transformation

[TABLE]

Time-reversal symmetry imposes an additional constraint on the gauge transformation, $\theta_{\tau}({\bm{k}})=-\theta_{\tau}(-{\bm{k}})$.

III Hartree-Fock calculation

Here, we provide the details for the Hartree-Fock calculation. Throughout this section, we neglect the Hund’s coupling term which is discussed in the previous section and drop the superscript $\pm$ from the form factor $\lambda$ such that $\lambda_{\tau,{\bm{q}}}({\bm{k}})=\lambda^{+}_{\tau,{\bm{q}}}({\bm{k}})$ since we only consider the intravalley form factor here.

We now move on to the general setup for the Hartree-Fock mean field theory. Define the expectation value

[TABLE]

which we will assume to be diagonal in ${\bm{k}}$ and ${\bm{k}}^{\prime}$, $M({\bm{k}},{\bm{k}}^{\prime})=\delta_{{\bm{k}},{\bm{k}}^{\prime}}M({\bm{k}})$. In the following, we will introduce the combined index $\alpha=(\sigma,\tau)$ such that $M({\bm{k}})$ is a matrix with components $M_{\alpha,\alpha^{\prime}}({\bm{k}})$. Next, we expand the interaction (Eq. S27) in the difference $c_{\alpha}^{\dagger}c_{\alpha^{\prime}}-M_{\alpha,\alpha^{\prime}}$ and neglect terms beyond linear order.

The resulting mean field Hamiltonian has the form

[TABLE]

Here, $c_{\bm{k}}$ is a column vector in the index $\alpha$, $\xi({\bm{k}})$ is a diagonal matrix containing the single particle energies $\xi_{\uparrow/\downarrow,\pm}({\bm{k}})$ and $h_{0,1}({\bm{k}})$ are 4 $\times$ 4 matrices in $\alpha$ given by

[TABLE]

and

[TABLE]

The matrix $\Lambda^{\pm}_{\bm{q}}({\bm{k}})$ simply contains the form factors defined in Eq. S30

[TABLE]

and $P_{\pm}=\frac{1}{2}(1\pm\tau_{z})$ is the projector on the $\pm$ valley with $\tau_{x,y,z}$ denoting the Pauli matrices in the valley space.

In both Eq. S36 and Eq. S37, the first term is a Hartree term whereas the second is a Fock term. Hartree terms were neglected in some of the previous mean-field studies Po2018 ; Zhang2018 since they are expected to couple only to the density which is determined by the filling in the gapped phase and is independent of the symmetry-breaking order. This is, however, not true in the presence of the form factors which are not the same for the two valleys $\lambda^{\pm}_{+,{\bm{q}}}({\bm{k}})\neq\lambda^{\pm}_{-,{\bm{q}}}({\bm{k}})$. As a result, the Hartree-term also couples to the valley density and it cannot be neglected.

It is important here to point out one major difference between our approach and the one employed recently in a self-consistent Hartree-Fock mean field study in twisted bilayer graphen Xie2018 . In that work, the Hartree-Fock corrections to the flat bands coming from all other ($\sim 150$) bands were taken into account. Here, we will instead make the assumption that the effect of the Hartree-Fock contributions from the other bands is already included at some level in the model parameters which should be either fit to experiments or obtained from ab initio studies at charge neutrality Jeil2014 ; Carr2019 . Thus, we only include the effects arising from filling the isolated band.

To write the self-consistency condition, we diagonalize $h_{0}({\bm{k}})+h_{1}({\bm{k}})$ by introducing the variables $d_{\bm{k}}=U_{\bm{k}}c_{\bm{k}}$ for some unitary $U_{\bm{k}}$. We then impose the constraint $M_{\alpha,\alpha^{\prime}}({\bm{k}})=\langle c_{\alpha}^{\dagger}({\bm{k}})c_{\alpha^{\prime}}({\bm{k}})\rangle$. In the following, we will only consider possible gapped phases at integer fillings $\nu$. In this case, the self-consistency condition has the form

[TABLE]

where $\chi$ is a ${\bm{k}}$-independent matrix containing $\nu$ ones along the diagonal and zeroes everywhere else. This means that $M({\bm{k}})$ is a projection operator satisfying

[TABLE]

Our assumption that the phase is gapped has to be checked self-consistently by computing the mean field band structure

[TABLE]

and ensuring that correlation induced gap for filling $\nu$ defined as

[TABLE]

is positive. Here, we assumed that the mean field bands $\epsilon_{\alpha,{\bm{k}}}$ are sorted in order of increasing energy. ($\alpha=1,2,3,4$)

We notice that $M({\bm{k}})$ is, in general, not gauge invariant. Instead it transforms as

[TABLE]

under the gauge transformation (Eq. S33). In the following mean field analysis, we will choose the gauge such that $\theta_{-}({\bm{k}})=\theta_{+}({\bm{k}})$ which guarantees the gauge independence of $M({\bm{k}})$.

III.1 Half-filling $\nu=2$

To understand the symmetry breaking at $\nu=2$, we notice that the order parameter can be written as

[TABLE]

$Q({\bm{k}})$ can then be expanded in terms of the generators $\sigma_{i}\tau_{j}$ as described in the main text. In the absence of inter-valley Hund’s coupling, the problem possesses an SU(2)$\times$SU(2) symmetry corresponding to independent spin rotations in each valley which are generated by $\sigma_{x,y,z}\tau_{0,z}$ in addition to UV(1) valley charge conservation generated by $\tau_{z}$ and time-reversal symmetry given by $\mathcal{T}=\tau_{x}K$. Inter-valley Hund’s coupling further breaks the SU(2)$\times$SU(2) to SU(2) corresponding to overall rotations. The generators can be grouped into 5 categories according to the symmetries they break as summarized in Table I. We notice that all these terms commute or anticommute with the generators of spin rotation $\sigma_{x,z}$, of UV(1) valley-charge conservation $\tau_{z}$ and with time-reversal symmetry. In fact, when considering possible symmetry broken states in the limit of flat bands and decoupled valleys, we can always restrict ourselves to matrices $Q({\bm{k}})$ which satisfy this requirement (for some choice of the generators of the symmetries). The reason is that such order are always energetically more favorable. To see this, consider a ’mixed’ order given by

[TABLE]

The Fock contribution to the mean-field energy is given by

[TABLE]

Substituting the mixed order (S45), we find that the mixed term containing both $Q_{1}$ and $Q_{2}$ has to vanish since there is some symmetry generator which commutes with $Q_{1}$ and anticommutes with $Q_{2}$ (note that the form factors are invariant under all symmetries). This implies that

[TABLE]

Since the Hartree-Fock solutions has to be extrema of the Hartree-Fock energy functional, we conclude that only pure orders which either commute or anticommute with each symmetry generator are possible self-consistent solution. This justifies restricting ourselves to the list of orders provided in Table I in the main text: SP, VP, SVL, IVC, and SIVCL (such order parameters are in general ${\bm{k}}$-dependent and may have more complicated forms than the ones written in the second column of the figure, but they have to respect the same symmetries).

If we first neglect the intervalley Hund’s coupling, we notice that the mean-field energies of the SP and SVL are equal as well as the IVC and SIVCL since they are related by rotating the spin in one of the valleys. Thus, in the following discussion, we can restrict ourselves to VP, SP, and IVC orders.

III.1.1 Valley polarized (VP) state

A valley polarized states breaks time-reversal but preserves spin rotation and valley charge. Together with the requirement that the order parameter has the form (S44), this yields

[TABLE]

The eigenvalues of $h_{\bm{k}}$ are given by

[TABLE]

We now need to check the correlated gap defined as

[TABLE]

is positive, so that the fully valley polarized state is a proper gapped state.

The total energy of the valley polarized state is obtained by adding the kinetic energy of the filled bands and the potential energy leading to

[TABLE]

III.1.2 Spin polarized (SP) state

Next we assume a spin polarized state along the $z$-direction in both valleys. Such state breaks spin-rotation but preserves time-reversal and valley charge conservation with the order parameter given by

[TABLE]

The energy eigenvalues are

[TABLE]

We also require the gap $\Delta_{\rm SP}=\min_{{\bm{k}},\tau}\epsilon_{\uparrow,\tau,{\bm{k}}}-\max_{{\bm{k}}^{\prime},\tau^{\prime}}\epsilon_{\uparrow,\tau^{\prime},{\bm{k}}^{\prime}}$ to be positive so that the spin polarized state is a proper gapped state.

The total energy of the spin polarized state is given by

[TABLE]

Comparing to the VP state, we find that the two phases have exactly the same ground state energy. This follows from time-reversal symmetry which implies that $\xi_{-}({\bm{k}})=\xi_{+}(-{\bm{k}})$ and $\lambda_{-,{\bm{G}}}({\bm{k}})=\lambda_{+,{\bm{G}}}(-{\bm{k}})$ as well as $|\lambda_{-,{\bm{q}}}({\bm{k}})|=|\lambda_{+,-{\bm{q}}}(-{\bm{k}})|$. Using the relation $\sum_{{\bm{k}}}\lambda_{+,{\bm{G}}}({\bm{k}})=\sum_{{\bm{k}}}\lambda_{-,{\bm{G}}}({\bm{k}})$, one can show that the total energies as well as the gaps are the same for VP and SP states.

III.1.3 Intervalley coherent (IVC) order

The intervalley coherent order parameter is given by

[TABLE]

We note that it is not possible in general to take the fully polarized limit in the $x-y$ plane and at the same time fulfill the self-consistency conditions. Hence, we include a small $z$ valley polarization parametrized by the angle $\theta_{\bm{k}}$. We notice that the state (Eq. S57) will not break time-reversal symmetry provided that $\theta_{-{\bm{k}}}=\pi-\theta_{\bm{k}}$ and $\phi_{-{\bm{k}}}=-\phi_{\bm{k}}$ which implies that the average valley polarization $\sum_{\bm{k}}\cos\theta_{\bm{k}}$ vanishes.

The mean field Hamiltonian has the form

[TABLE]

with $f_{\bm{k}}$, $A_{\bm{k}}$, $B_{\bm{k}}$ given by

[TABLE]

The self-consistency condition reads

[TABLE]

where energy eigenvalues are given by

[TABLE]

with the gap given by $\Delta_{\rm IVC}=\min_{{\bm{k}}}\epsilon_{+,{\bm{k}}}-\max_{{\bm{k}}^{\prime}}\epsilon_{-,{\bm{k}}^{\prime}}$ which should be positive for a proper gapped phase. The results of the energy competition between the VP/SP phase and the IVC state obtained by numerically solving the self-consistency equation are given in the main text.

III.1.4 Effect of intervalley Hund’s coupling

As we have seen above, the three distinct states with spin polarization, valley polarization or spin-valley locking are degenerate in the absence of intervalley Hund’s and their energy is always lower than the energy of the valley off-diagonal orders (IVC and SIVCL). In the following, we want to investigate the effect of intervalley Hund’s coupling on these three states. Since this term $J$ is much smaller than the main part of the interaction $V_{0}$, it suffices to compute it for the three valley-diagonal orders since the valley off-diagonal orders are already energetically unfavorable on the level of $V_{0}$. Substituting in (S37) we find that

[TABLE]

III.2 Quarter-filling $\nu=1$

At quarter filling, $\nu=1$ (similarly for $\nu=3$, with some caveats), we can always write the order parameter as

[TABLE]

This leads to three distinct possibilities: (i) $Q_{1}=\sigma_{x,y,z}\tau_{0}$ and $Q_{2}=\sigma_{0}\tau_{z}$ which corresponds to a spin and valley polarized state, (ii) $Q_{1}=\sigma_{x,y,z}\tau_{0}$ and $Q_{2}=\sigma_{0}\tau_{x,y}$ which corresponds to a spin-polarized IVC, and (iii) $Q_{1}=\sigma_{z}\tau_{z}$ and $Q_{2}=\sigma_{x}\tau_{x},\sigma_{y}\tau_{y}$ which correspond to a spin-valley locked IVC.

The SPIVC and SVLIVC are related by a spin rotation in one of the valleys, thus we can focus only on the competition between SVP and SPIVC. Compared to the VP vs IVC states at half-filling these differ by a factor of 2 in the Fock energy and a factor of 4 in the Hartree energy. Since the former is the main deciding factor in the competition between the phases, the results for the gaps and energy difference between SVP and SPIVC at quarter filling are very similar to those between VP and IVC at half-filling.

IV Perturbative solution and competition between VP/SP and IVC

In this section, we would like to discuss the competition between inter-valley coherent order and valley/spin polarized order in a more general setting that is not too sensitive to the details of the model parameters. To this end, it is useful to derive an approximate solution to the self-consistency equations and compute an analytic expression for the energy difference between the IVC phase and the VP/SP phase.

In order to make progress analytically, we can write the IVC order parameter as

[TABLE]

where $\gamma_{\bm{k}}\sim\beta_{\bm{k}}\sim\delta\ll 1$. This approximation can be justified as follows: the starting symmetry of the isolated band is SU(4) which is broken to SU(2) $\times$ SU(2) due to the asymmetry between the two valley in energies and form factors ($\xi_{+}({\bm{k}})\neq\xi_{-}({\bm{k}})$, $\lambda_{+,{\bm{q}}}({\bm{k}})\neq\lambda_{-,{\bm{q}}}({\bm{k}})$). In the following, we will assume that breaking SU(4) to SU(2) $\times$ SU(2) is not very strong so that the deviation from the situation where the valleys are identical is weak. This condition can be written more explicitly as the requirement that $\frac{|\xi_{+}({\bm{k}})-\xi_{-}({\bm{k}})|}{V_{0}}\sim|\lambda_{+,{\bm{q}}}({\bm{k}})-\lambda_{-,{\bm{q}}}({\bm{k}})|\sim\delta\ll 1$. The first part is guaranteed by the small bandwidth whereas the second one can be checked numerically and shown to hold at least for most values of ${\bm{k}}$ and ${\bm{q}}$. This is equivalent to expanding in time-reversal symmetry breaking terms within each valley.

The variables $\gamma_{\bm{k}}$ and $\beta_{\bm{k}}$ can be obtained by solving a linearized version of the self-consistency equation as follows. We start by expanding $\theta_{\bm{k}}$ and $\phi_{\bm{k}}$ in terms of small deviations $\delta$ from a perfect IVC state in the $\tau_{x}$ as shown in Eq. S66. Substituting in Eq. S62 and expanding to leading order in $\delta$ yields the following set of linear equations given by

[TABLE]

where $a_{\bm{k}}$ is given by

[TABLE]

and $F_{{\bm{k}},{\bm{k}}^{\prime}}$ and $b_{\bm{k}}$ are given by

[TABLE]

We notice that $a_{\bm{k}}$ and $\imaginary b_{\bm{k}}$ are of order $\delta$. Substituting in the expression for the energy, the energy difference between the IVC state and the SP/VP state can be written (up to second order in $\delta$) as

[TABLE]

where $F_{{\bm{k}},{\bm{k}}^{\prime}}$ and $b_{\bm{k}}$ are defined in (S69). The first term in Eq. S70 reproduces the non self-consistent Hartree-Fock energies obtained in Ref. Zhang2018 in which case VP/SP is always favored to IVC.

The second and third terms are corrections coming from solving the self-consistency condition. It is instructive to reproduce the results of Ref. Po2018 which considers the simplified setting where all form factors are taken equal to 1. In addition, $v_{\bm{q}}$ was taken equal to a constant which is cutoff at large momenta ${\bm{q}}\sim\Lambda$ yielding the interaction strength $g=\frac{V_{0}}{2N}\sum_{|{\bm{q}}|<\Lambda}=\frac{V_{0}}{2}\sum_{|{\bm{G}}|<\Lambda}$. In this case, $\gamma_{\bm{k}}=\frac{\xi_{+}({\bm{k}})-\xi_{-}({\bm{k}})}{2g}$ and $\beta_{\bm{k}}=0$ leading to

[TABLE]

which implies that the IVC phase is energetically favored to the VP/SP phase in agreement with the conclusion of Ref. Po2018 222The result differs by a factor of 2 due to the incorrect way the large $g$ limit was implemented Po2018 ..

Our result (Eq. S70) interpolates between these two limits with the first two terms favoring spin or valley polarization and the last term favoring intervalley coherence. The competition between SP/VP and IVC is then settled by the details of the band structure, form factors and interaction. We notice that this expression underestimates the energy of the IVC states when the bands have non-zero Chern number. In this case, it was shown in Ref. Bultinck19 that vortices in the IVC order parameter are unavoidable. The existence of vortices is neglected in the expansion (Eq. S66) which assumes that $\phi_{\bm{k}}$ is small everywhere. This implies that the expression (Eq. S70) underestimates the IVC ground state energy for non-zero Chern number.

In order to gain some insights about what parameters control this competition, let us consider a very simplified setting where the Berry curvature is uniform in momentum space with the form factor assuming the simple form Zhang2018

[TABLE]

Here, $C$ is the Chern number for the $+$ valley and the parameter $\alpha$ determines how quickly the form factor decays with ${\bm{q}}$ which we take equal to ${2\pi}/{A_{\rm BZ}}$ to reproduce the Landau level form factors for $C=1$. These form factors would be obtained in a Landau level if it is folded into a Brillouin zone of the lattice with the flux density one usov1988theory . In addition, we will consider a very simple form of the dispersion corresponding to nearest neighbour tight-binding model on a triangular lattice with hopping amplitude $te^{\pm i\phi}$ for the $\pm$ valleys. For $C=0$, we know that it is possible to write such a tight-binding model. For non-zero $C$, it is generally impossible to write such tight binding model. However, we can still use the same resulting dispersion and assume that the non-zero Chern number only affects the form factors. This will enable us to disentangle the effects of the band dispersion from those related to band topology.

For the form factors given in Eq. S72, the self-consistency equations can be solved by performing Fourier transform to real space. Following a series of straightforward steps, we get

[TABLE]

where $I_{b}(x)$ is the modified Bessel function of the first kind and $\eta=(t/V_{0})\sin\phi$. The proportionality here indicates that we have dropped a constant positive factor given by $\frac{\pi}{|A_{\rm BZ}|}\sqrt{\frac{\pi}{\alpha}}$ which does not influence the competition between the two phases.

The expression (Eq. S73) depends only on two dimensionless parameters: (i) the Chern number $C$ and (ii) $\eta$ which measures the bandwidth relative to the interaction strength multiplied by the strength of time-reversal symmetry breaking within each valley. The first term in Eq. S73 is always positive and favors SP/VP state. It vanishes for zero Chern number and increases as the Chern number increases. This suggests that increasing the Chern number favors valley/spin polarization over intervalley coherent order. The second term, on the other hand, favors IVC and increases with increasing the bandwidth or the time-reversal symmetry breaking within each valley.

The phase diagram for different values of $C$ and $\eta$ is given in Fig. S2. For $C=0$, IVC order always wins. This is an artifact of our simple choice for the form factors which corresponds to uniform Berry curvature. In a more realistic model where the Berry curvature vanishes on average but does not vanish everywhere, we expect some region of VP/SP. This is expected to be particularly manifest in the vicinity of topological phase transitions where the valley Chern number changes leading to a large concentration of the Berry curvature at some momenta. For $C\neq 0$, we find that VP/SP is always favored for relatively small values of the bandwidth whereas IVC is favored for relatively large values. Since our approach underestimates the IVC energy for non-zero Chern number (since it ignores vortices Bultinck19 ), we expect the transition from VP/SP to IVC to happen at even larger values of $\eta$ implying that VP/SP is the most energetically favorable insulator at half-filling whenever the bandwidth is relatively narrow.

V Spin-triplet superconductivity

In the following, we provide some details on the discussion related to spin-triplet superconductivity in the main text. The interaction term

[TABLE]

can be rewritten as

[TABLE]

with $g_{\alpha,\beta;\gamma,\delta}=\sum_{a}\sigma_{\alpha\beta}^{a}\sigma_{\gamma\delta}^{a}$. When performing the BCS decoupling, we restrict ourselves to pairing between time-reversed pairing which corresponds to ${\bm{k}}^{\prime}=-{\bm{k}}$ and $\tau^{\prime}=-\tau$. In this case, we can define the gap function $\Delta_{\alpha,\beta,\tau,\tau^{\prime}{\bm{k}}}=\delta_{\tau,-\tau^{\prime}}\langle c^{\dagger}_{\alpha,\tau,{\bm{k}}}c^{\dagger}_{\beta,\tau^{\prime},-{\bm{k}}}\rangle$, which satisfies the linearized BCS equation

[TABLE]

where $v_{F}({\bm{k}})$ is the Fermi velocity at point ${\bm{k}}$ on the Fermi s surface $v_{F}({\bm{k}})=|\nabla_{\bm{k}}\epsilon_{\bm{k}}|$. Choosing $\Delta_{\bm{k}}$ to be ${\bm{k}}$-independent, we can simplify (Eq. S76)

[TABLE]

where $\tilde{\lambda}$ is related to $\lambda$ by some constant rescaling (coming from the Fermi surface integral), ${\bm{\sigma}}$ is the Pauli matrix vector in spin space and $\Delta$ is a matrix in spin and valley spaces. As discussed in the main text, intervalley pairing is proportional to $\tau_{x}$ or $\tau_{y}$ which corresponds to valley triplet or singlet respectively, which, due to the overall antisymmetry of the gap function, implies the former scenario corresponds to a spin-singlet $i\sigma_{y}$ whereas the latter corresponds to a spin-triplet $i\sigma_{y}{\bm{d}}\cdot{\bm{\sigma}}$. Here, ${\bm{d}}$ is the vector which captures the direction of the spin state.

The symmetry of the superconducting order parameter is obtained by finding the pairing channel for which $\tilde{\lambda}$ is positive and maximum. Substituting the spin-singlet and triplet gap functions in (Eq. S77) yields

[TABLE]

which implies a valley-singlet spin-triplet superconductor.

VI Dependence of $T_{c}$ on magnetic field

In the following, we will write a simple mean field theory to relate the parameters in the Ginzburg-Landau free energy in Eq. 6 of the main text to the microscopic parameters. We start by writing the following imaginary time mean-field action

[TABLE]

Here, $\psi$ is a (grassman-valued) spinor in valley and spin spaces, $\sigma$ and $\tau$ are Pauli matrices for the spin and valley degrees of freedom, respectively. $\chi$ is dimensionless magnetic susceptibility and $\Delta$ is a matrix in the valley and spin spaces. Following the discussion of the main text, we take $\Delta_{\bm{k}}$ to be ${\bm{k}}$-independent, spin-triplet and valley singlet

[TABLE]

The magnetic field enters (Eq. S80) through Zeeman and orbital couplings with the ${\bm{k}}$-dependent $g$-factor arising form the orbital effect (see the main text). (In (Eq. S80), ${\bm{g}}_{\bm{k}}$ is a diagonal matrix in spin and valley spaces given by $\sigma_{0}{\rm diag}({\bm{g}}_{+,{\bm{k}}},{\bm{g}}_{-,{\bm{k}}})_{\tau}$). If the parent state is either a weak ferromagnet or close to a ferromagnetic quantum critical point which we anticipate to be the case, then dimensional suscepbtility $\chi$ can be relatively large and cannot be put to 1.

We can go now to matsubara frequency by writing

[TABLE]

leading to

[TABLE]

Here, we introduced the Green’s function $G_{p}$ as

[TABLE]

where $\xi_{{\bm{k}}}$ depends on the valley index such that $\xi_{+,-{\bm{k}}}=\xi_{-,{\bm{k}}}$. The fermions can be integrated out leading to a Pfafian which can be written in the exponential as the logarithm of the trace of some operator. The resulting free energy can be expanded in powers of ${\bm{B}}$ and $\Delta$.

The term proportional to $\Delta\Delta^{\dagger}$ provides the standard BCS instability which is given by

[TABLE]

where $f(\xi)$ is defined as

[TABLE]

The integral over $\xi$ is cut off by the bandwidth $\Lambda$. However, we can choose to perform the $\xi$ integral before the frequency sum in which case, the integral is automatically cut off by $\omega_{n}$ so that it can be extended to infinity with the cutoff $\Lambda$ moved to the $\omega_{n}$ sum instead. This leads to

[TABLE]

The final result is given by

[TABLE]

The term proportional to $\Delta\Delta^{\dagger}{\bm{B}}\cdot{\bm{\sigma}}$ is given by

[TABLE]

The linear term corresponding to the orbital effect $\Delta\Delta^{\dagger}{\bm{B}}\cdot{\bm{g}}$ vanishes due to time-reversal symmetry which can be seen as follows

[TABLE]

The last equality follows from the fact that ${\bm{g}}_{\tau,{\bm{k}}}$ is odd under time-reversal symmetry ${\bm{g}}_{+,-{\bm{k}}}=-{\bm{g}}_{-,{\bm{k}}}$.

The term proportional to $\Delta\Delta^{\dagger}({\bm{B}}\cdot{\bm{\sigma}})^{2}$ is given by

[TABLE]

The first term can be simplified by noting that

[TABLE]

whereas the second term can be evaluated as

[TABLE]

leading to

[TABLE]

The term proportional to $\Delta\Delta^{\dagger}({\bm{B}}\cdot{\bm{g}})^{2}$ is given by

[TABLE]

Here, we used ${\bm{g}}_{\sigma,{\bm{k}}}=-{\bm{g}}_{-\sigma,-{\bm{k}}}$ to go from the second to the third line and (Eq. S92) and (Eq. S93) to go from the third to the fourth line.

Finally, we evaluate the quartic term $(\Delta^{\dagger}\Delta)^{2}$ as

[TABLE]

The summation over $p$ is given by (S93), whereas the trace can be evaluated as

[TABLE]

leading to

[TABLE]

The Free energy now has the form

[TABLE]

We notice that the second derivative of the density of states can be estimated as $1/\epsilon_{F}^{2}$ which is much smaller that $\beta^{2}$, thus we can throw away all terms containing $N^{\prime\prime}(0)$. Expanding in $T$ close to $T_{c}=\Lambda e^{-\frac{1}{gN(0)}}$, we get

[TABLE]

Comparing with Eq. 6 in the main text, we find that the coefficients $\kappa,a,b,c,\alpha,\eta$ are given by

[TABLE]

where ${\bm{e}}_{{\bm{B}}}$ is the direction of the external magnetic field. We notice that the term $a$ was obtained in the description of the superfluid transition in He3 Ambegaokar73 .