Ferromagnetism and superconductivity in twisted double bilayer graphene

TL;DR

This paper develops a theoretical framework for understanding the coexistence and competition of ferromagnetism and superconductivity in twisted double bilayer graphene, explaining experimental observations and predicting phase behaviors.

Contribution

It introduces a model linking Coulomb-driven ferromagnetism and phonon-mediated superconductivity in TDBG, with calculations of transition temperatures and phase stability.

Findings

Superconducting domes appear on both electron and hole sides of ferromagnetic insulator.

Ferromagnetic insulating gap exhibits a dome shape dependence on layer potential difference.

Half-filled ferromagnetic insulator is stable against spin and valley magnons.

Abstract

We present a theory of competing ferromagnetic and superconducting orders in twisted double bilayer graphene (TDBG). In our theory, ferromagnetism is induced by Coulomb repulsion, while superconductivity with intervalley equal-spin pairing can be mediated by electron-acoustic phonon interactions. We calculate the transition temperatures for ferromagnetism and superconductivity as a function of moir\'e band filling factor, and find that superconducting domes can appear on both the electron and hole sides of the ferromagnetic insulator at half filling. We show that the ferromagnetic insulating gap has a dome shape dependence on the layer potential difference, which provides an explanation to the experimental observation that the ferromagnetic insulator only develops over a finite range of external displacement field. We also verify the stability of the half-filled ferromagnetic insulator against two types of collective excitations, i.e., spin magnons and valley magnons.

Equations

$${\mathcal{H}}_0(\mathbf{k})=\left(\begin{array}{cccc}{\Delta }_0& \hslash v_1{k}_{-}& \hslash v_2{k}_{+}& {\gamma }_1\\ \hslash v_1{k}_{+}& 0& \hslash v_3{k}_{-}& \hslash v_2{k}_{+}\\ \hslash v_2{k}_{-}& \hslash v_3{k}_{+}& 0& \hslash v_1{k}_{-}\\ {\gamma }_1& \hslash v_2{k}_{-}& \hslash v_1{k}_{+}& {\Delta }_0\end{array}\right),$$

$$(v_1,v_2,v_3)=(0.844,-0.045,-0.091)\times 10^6\text{m/s},$$

$${\gamma }_1=361\text{meV},{\Delta }_0=15\text{meV},$$

$${\mathcal{H}}_{+K}=\left(\begin{array}{cc}{h}_{\mathfrak{b}}(\mathbf{k})& \tilde{T}(\mathbf{r})\\ {\tilde{T}}^{†}(\mathbf{r})& h_{\mathfrak{t}}(\mathbf{k})\end{array}\right),$$

$$\tilde{T}(\mathbf{r})$$

$$T(\mathbf{r})$$

$$T_j$$

$$H_0=\sum _{\mathbf{k},\tau ,s}{\epsilon }_{\mathbf{k},\tau }{c}_{\mathbf{k},\tau ,s}^{†}{c}_{\mathbf{k},\tau ,s},$$

$$H_1=\frac{1}{2\mathcal{A}}\sum V_{{\mathbf{k}}_1{\mathbf{k}}_2{\mathbf{k}}_3{\mathbf{k}}_4}^{(\tau {\tau }^{\prime })}{c}_{{\mathbf{k}}_1,\tau ,s}^{†}{c}_{{\mathbf{k}}_2,{\tau }^{\prime },s^{\prime }}^{†}{c}_{{\mathbf{k}}_3,{\tau }^{\prime },s^{\prime }}{c}_{{\mathbf{k}}_4,\tau ,s},$$

$$V_{{\mathbf{k}}_1{\mathbf{k}}_2{\mathbf{k}}_3{\mathbf{k}}_4}^{(\tau {\tau }^{\prime })}=\sum _{\mathbf{q}}V(\mathbf{q})O_{{\mathbf{k}}_1{\mathbf{k}}_4}^{(\tau )}(\mathbf{q})O_{{\mathbf{k}}_2{\mathbf{k}}_3}^{({\tau }^{\prime })}(-\mathbf{q}),$$

$$O_{{\mathbf{k}}_1{\mathbf{k}}_4}^{(\tau )}(\mathbf{q})=\sum _{\sigma ,\ell }\int d\mathbf{r}{e}^{i\mathbf{q}\cdot \mathbf{r}}{\Phi }_{\tau ,{\mathbf{k}}_1,\sigma ,\ell }^{\ast }(\mathbf{r}){\Phi }_{\tau ,{\mathbf{k}}_4,\sigma ,\ell }(\mathbf{r}),$$

$$H_{\text{MF}}=$$

$$E_{\mathbf{k},\tau ,s}=$$

$${\Sigma }_{\mathbf{k},\tau ,s}=$$

$$-\frac{1}{A}\sum _{{\mathbf{k}}^{\prime }}{V}_{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{k}{\mathbf{k}}^{\prime }}^{(\tau \tau )}{n}_F(E_{{\mathbf{k}}^{\prime },\tau ,s})$$

$${\Sigma }_{\mathbf{k},\tau }^{(-)}$$

$$M_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\tau )}$$

$$H_{\text{att}}=-g_0\sum _{\sigma ,{\sigma }^{\prime },\ell ,s}\int d\mathbf{r}{\hat{\psi }}_{+\sigma \ell s}^{†}{\hat{\psi }}_{-{\sigma }^{\prime }\ell s^{\prime }}^{†}{\hat{\psi }}_{-{\sigma }^{\prime }\ell s^{\prime }}{\hat{\psi }}_{+\sigma \ell s},$$

$$H_{\text{p}}$$

$$g_{\mathbf{k}{\mathbf{k}}^{\prime }}$$

$$H_{\text{BCS}}$$

$${\Delta }_{\mathbf{k},s}$$

$${\chi }_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(s)}$$

$$H_2$$

$$=\frac{U_0}{N}\sum _{{\mathbf{k}}_i}^{{}^{\prime }}\sum _{\sigma }{b}_{{\mathbf{k}}_1\sigma \uparrow }^{†}{b}_{{\mathbf{k}}_2\sigma \downarrow }^{†}{b}_{{\mathbf{k}}_3\sigma \downarrow }{b}_{{\mathbf{k}}_4\sigma \uparrow },$$

$$H_2\approx \frac{U_0}{N}\sum _{{\mathbf{k}}_i}^{{}^{\prime }}\sum _{{\tau }_i}^{{}^{\prime }}\sum _{\sigma }{b}_{{\mathbf{k}}_1{\tau }_1\sigma \uparrow }^{†}{b}_{{\mathbf{k}}_2{\tau }_2\sigma \downarrow }^{†}{b}_{{\mathbf{k}}_3{\tau }_3\sigma \downarrow }{b}_{{\mathbf{k}}_4{\tau }_4\sigma \uparrow },$$

$$b_{\mathbf{k}\tau \sigma s}^{†}=\frac{1}{\sqrt{\mathcal{A}}}\int d\mathbf{r}{e}^{i\mathbf{k}\cdot \mathbf{r}}{\hat{\psi }}_{\tau \sigma s}^{†}(\mathbf{r}),$$

$$H_2\approx u_0\sum _{{\tau }_i}^{{}^{\prime }}\sum _{\sigma }\int d\mathbf{r}{\hat{\psi }}_{{\tau }_1\sigma \uparrow }^{†}{\hat{\psi }}_{{\tau }_2\sigma \downarrow }^{†}{\hat{\psi }}_{{\tau }_3\sigma \downarrow }{\hat{\psi }}_{{\tau }_4\sigma \uparrow },$$

$$\displaystyle\tilde{H}_{2}=\sum_{\tau,\boldsymbol{k}}\big{\{}$$

$$+$$

$$u_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\tau {\tau }^{\prime })}=u_0\mathcal{A}\sum _{\sigma \ell }\int d\mathbf{r}|{\Phi }_{\tau ,\mathbf{k},\sigma ,\ell }(\mathbf{r}){|}^2|{\Phi }_{{\tau }^{\prime },{\mathbf{k}}^{\prime },\sigma ,\ell }(\mathbf{r}){|}^2.$$

$$|\mathbf{Q}{⟩}_S=\sum _{\mathbf{k}}{z}_{\mathbf{k},\mathbf{Q}}{c}_{\mathbf{k}+\mathbf{Q},+,\downarrow }^{†}{c}_{\mathbf{k},+,\uparrow }|\text{FM}⟩$$

$${\mathcal{E}}_S(\mathbf{Q})z_{\mathbf{k},\mathbf{Q}}=\sum _{{\mathbf{k}}^{\prime }}{\mathcal{H}}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\mathbf{Q})}{z}_{{\mathbf{k}}^{\prime },\mathbf{Q}},$$

$${\mathcal{H}}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\mathbf{Q})}=(E_{\mathbf{k}+\mathbf{Q},+,\downarrow }-E_{\mathbf{k},+,\uparrow }){\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }}-\frac{1}{A}{V}_{{\mathbf{k}}^{\prime }(\mathbf{k}+\mathbf{Q})({\mathbf{k}}^{\prime }+\mathbf{Q})\mathbf{k}}^{(++)},$$

$$|\mathbf{Q}{⟩}_V=\sum _{\mathbf{k}}{z}_{\mathbf{k},\mathbf{Q}}{c}_{\mathbf{k}+\mathbf{Q},+,\downarrow }^{†}{c}_{\mathbf{k},-,\uparrow }|\text{FM}⟩.$$

$${\mathcal{E}}_V(\mathbf{Q})z_{\mathbf{k},\mathbf{Q}}=\sum _{{\mathbf{k}}^{\prime }}{\mathcal{W}}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\mathbf{Q})}{z}_{{\mathbf{k}}^{\prime },\mathbf{Q}},$$

$${\mathcal{W}}_{\mathbf{k}{\mathbf{k}}^{\prime }}^{(\mathbf{Q})}=(E_{\mathbf{k}+\mathbf{Q},+,\downarrow }-E_{\mathbf{k},-,\uparrow }){\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }}-\frac{1}{A}{V}_{{\mathbf{k}}^{\prime }(\mathbf{k}+\mathbf{Q})({\mathbf{k}}^{\prime }+\mathbf{Q})\mathbf{k}}^{(-+)},$$

Full text

I introduction

Moiré superlattices form in van der Waals bilayers with a small orientation misalignment and/or lattice constant mismatch. Recently moiré bilayers have emerged as a platform to study fundamental physics of strongly interacting systems, in view of the discovery of correlated insulating and superconducting states in twisted bilayer grapheneCao et al. (2018a, b). Moiré superlattices often generate spatial confinement for low-energy electrons, suppress electron kinetic energy, and therefore effectively enhance interaction effects. Evidences of correlated insulating and superconducting states have so far been reported in three graphene-based moiré systems, including twisted bilayer graphene Cao et al. (2018a, b); Yankowitz et al. (2019); Kerelsky et al. ; Choi et al. ; Sharpe et al. ; Cao et al. (a); Polshyn et al. ; Codecido et al. ; Lu et al. ; Tomarken et al. , twisted double bilayer grapheneShen et al. ; Liu et al. (a); Cao et al. (b); Burg et al. (2019); He et al. , and ABC trialyer graphene on hexagonal boron nitrideChen et al. (2019, a, b) .

Twisted bilayer graphene (TBG) is a subject under intense theoretical studyXu and Balents (2018); Po et al. (2018); Koshino et al. (2018); Kang and Vafek (2018); Dodaro et al. (2018); Padhi et al. (2018); Guo et al. (2018); Fidrysiak et al. (2018); Kennes et al. (2018); Liu et al. (2018); Isobe et al. (2018); You and Vishwanath (2019); Tang et al. (2019); Rademaker and Mellado (2018); Xu et al. (2018); Guinea and Walet (2018); Carr et al. (2018); Thomson et al. (2018); González and Stauber (2019); Su and Lin (2018); Ramires and Lado (2018); Tarnopolsky et al. (2019); Ahn et al. (2019); Song et al. (2018); Lian et al. ; Hejazi et al. (2019); Liu et al. (2019); Sherkunov and Betouras (2018); Venderbos and Fernandes (2018); Kozii et al. (2019); Kang and Vafek ; Xie and MacDonald ; Lin and Nandkishore ; Bultinck et al. ; Alidoust et al. (2019); Angeli et al. ; Peltonen et al. (2018); Lian et al. (2018); Choi and Choi (2018); Wu et al. (2018, 2019); Wu (2019); Wu and Das Sarma , but the exact nature of the correlated insulating (CI) and superconducting (SC) states in TBG remains unsettled. The half-filled correlated insulator in TBG crosses over to a metallic state by a strong perpendicular or parallel magnetic field Cao et al. (2018a, b), which possibly rules out spin-polarized ferromagnetic states, but leaves a large number of possible non-FM states as candidates, e.g., valley polarized state, and charge/spin/valley density wave states to name a few.

By contrast, there appears to be good experimental evidence that the half-filled CI in twisted double bilayer graphene (TDBG) with a twist angle $\theta$ around $1.3^{\circ}$ is ferromagnetic, because the correlation driven insulating gap has been found to be enhanced by an in-plane magnetic field Shen et al. ; Liu et al. (a); Cao et al. (b); Burg et al. (2019); He et al. . Possible signature of SC domes in adjacent to the CIs has also been reported in TDBG Shen et al. ; Liu et al. (a). These experimental discoveries identify TDBG as another important moiré system with strong interaction effects. Moreover, TDBG represents a simpler as well as a more tunable system compared to TBG, because moiré bands of TDBG can be controlled by an out-of-plane electric displacement field and its first moiré conduction band can be energetically isolated from neighboring bands, whereas the first moiré valence and conduction bands in TBG are typically connected via Dirac points.

In this paper, we theoretically study TDBG FM and SC orders in its first moiré conduction band. In our theory, ferromagnetism is driven by Coulomb repulsion as in Stoner model, but superconductivity is mediated by electron-phonon interactions. FMCI can occur at half filling when spin majority and minority bands are separated in energy by Coulomb exchange interaction. We find that the FM insulating gap is tunable by a layer potential difference $U$ that is generated by an external out-of-plane displacement field. This tunability originates from the strong dependence of the moiré bands on $U$, and agrees with experimental observations Shen et al. ; Liu et al. (a); Cao et al. (b); Burg et al. (2019); He et al. . We also calculate spin and valley magnon spectrum and verify the stability of the FMCI.

Away from half filling, the state is generally metallic, which can be susceptible to superconducting instability at low temperature. Because electron-acoustic phonon interactions in graphene mediate both spin singlet and spin triplet intervalley Cooper pairingWu et al. (2019), superconductivity can take place even in the presence of ferromagnetism, as long as the spinless time-reversal symmetry is preserved. We estimate the transition temperatures of ferromagnetism and superconductivity as a function of filling factor, and find that superconducting domes can appear on both sides of the half-filled ferromagnetic insulator.

Our paper is organized as follows. We describe the moiré Hamiltonian and band structure of TDBG in Section  II. We find that the moiré bands are tunable by $U$, and the van Hove singularity in the non-interacting density of states can be tuned from below to above half filling of the first moiré conduction bands. This feature allows the controlling of many-body physics using the out-of-plane displacement field. We present theory for ferromagnetism and superconductivity, respectively, in Sections III and IV, and make a brief conclusion in Section  V. Some technical details of the theory are given in Appendices A and B.

II Moiré Bands

We study TDBG with a small twist angle $\theta$ relative to the AB-AB stacking configuration, and calculate the moiré band structure using a continuum Hamiltonian generalized from TBGBistritzer and MacDonald (2011) to TDBG Zhang et al. (2019); Chebrolu et al. ; Koshino (2019); Lee et al. ; Liu et al. (b); Haddadi et al. . Within the continuum approximation, $\pm K$ valleys are treated separately. For each AB bilayer graphene, we use the following $k.p$ Hamiltonian in $+K$ valley

[TABLE]

which is in the basis of $A_{1}$, $B_{1}$, $A_{2}$ and $B_{2}$ sites [Fig. 1(a)] from one AB bilayer graphene. $k_{\pm}$ stands for $k_{x}\pm ik_{y}$. Parameter values are taken as

[TABLE]

which are extracted from ab initio results of Ref. Jung and MacDonald, 2014. The moiré Hamiltonian in $+K$ valley is given by

[TABLE]

where $h_{\mathfrak{b}}(\boldsymbol{k})$ and $h_{\mathfrak{t}}(\boldsymbol{k})$ are $k.p$ Hamiltonians for bottom and top bilayer graphene, and are equal to $\mathcal{H}_{0}[\hat{R}(+\theta/2)(\boldsymbol{k}-\boldsymbol{\kappa}_{+})]$ and $\mathcal{H}_{0}[\hat{R}(-\theta/2)(\boldsymbol{k}-\boldsymbol{\kappa}_{-})]$, respectively. Here $\hat{R}(\pm\theta/2)$ are rotation matrices and $\boldsymbol{\kappa}_{\pm}=[4\pi/(3a_{M})](-\sqrt{3}/2,\mp 1/2)$. The moiré period $a_{M}$ is approximately $a_{0}/\theta$, where $a_{0}$ is the monolayer graphene lattice constant. $\tilde{T}(\boldsymbol{r})$ is the tunneling between bottom and top bilayer graphene, which varies spatially with the moiré period as specified by

[TABLE]

where we only keep tunneling terms between adjacent layers, and $\boldsymbol{b}_{\pm}$ are moiré reciprocal lattice vectors given by $[4\pi/(\sqrt{3}a_{M})](\pm 1/2,\sqrt{3}/2)$. $w_{0}$ and $w_{1}$ are two tunneling parameters, which in general have different numerical values due to layer corrugation in the moiré pattern Koshino et al. (2018); Choi and Choi . We take $w_{0}=88$ meV and and $w_{1}=100$ meV. An out-of-plane electric displacement field generates a layer dependent potential, which can be parametrized using a single parameter $U$ as illustrated in Fig.1(b). The point group symmetry of TDBG is $D_{3}$ in the absence of the displacement field ($U=0$), and is broken down to $C_{3}$ when $U$ is finite.

A representative moiré band structure is shown in Fig. 2 for $\theta=1.24^{\circ}$ and $U=45$ meV. The first conduction band in Fig. 2(a) is isolated in energy from other bands, narrow in bandwidth($\sim$ 13 meV), and topologically nontrivial with a Chern number of $+2$ in $+K$ valley. Because of time-reversal symmetry, the corresponding moiré band in $-K$ valley has the opposite Chern number.

We find that the band dispersion can be drastically controlled by the potential $U$, as demonstrated by the energy contour plots of the first moiré conduction band shown in Fig. 3. In particular, the van Hove saddle points can be effectively moved in the momentum space by tuning $U$. At a critical $U\approx 44$ meV, three van Hove saddle points merge to the corner of the moiré Brillouin zone (MBZ), forming a high-order saddle point Yuan et al. (2019). Correspondingly, the density of states (DOS) for the non-interacting band has a strong dependence on $U$, and the van Hove singularity in the DOS can be tuned from below to above half filling by varying $U$ [Fig. 2(b)]. This strong dependence of the moiré bands on $U$ has implications on interaction physics, as we explain in the following.

III Ferromagnetism

III.1 Ferromagnetic Ground State

Many-body interactions are effectively enhanced for electrons in the nearly flat moiré bands. Here we study flatband ferromagnetism driven by Coulomb repulsion using a momentum-space formalism Zhang et al. (2019); Lee et al. ; Chen et al. (b). We only retain the first conduction band for simplicity, and the single-particle Hamiltonian projected onto this band is

[TABLE]

where $\boldsymbol{k}$ is momentum measured relative to the center of the MBZ, $\tau=\pm$ is the valley index, $s$ represents spin ($\uparrow$, $\downarrow$) and $c^{\dagger}_{\boldsymbol{k},\tau,s}$ is the fermion creation operator. $\varepsilon_{\boldsymbol{k},\tau}$ is the spin independent moiré band energy; its valley dependence is determined by time reversal symmetry, and $\varepsilon_{\boldsymbol{k},\tau}=\varepsilon_{-\boldsymbol{k},-\tau}$.

We project Coulomb interaction onto the first conduction band, and the interacting Hamiltonian has the form

[TABLE]

where $\mathcal{A}$ is the system area and $\Phi_{\tau,\boldsymbol{k}}(\boldsymbol{r})$ is the Bloch wave function for the first conduction band in valley $\tau K$ and at momentum $\boldsymbol{k}$. The indices $\sigma$ and $\ell$ respectively label sublattices and layers. By time-reversal symmetry, $\Phi_{\tau,\boldsymbol{k}}(\boldsymbol{r})=\Phi^{*}_{-\tau,-\boldsymbol{k}}(\boldsymbol{r})$. In the plane wave matrix element $O_{\boldsymbol{k}_{1}\boldsymbol{k}_{4}}^{(\tau)}(\boldsymbol{q})$, the momentum $\boldsymbol{q}$ can differ from $\boldsymbol{k}_{1}-\boldsymbol{k}_{4}$ by moiré reciprocal lattice vectors. Hamiltonian $H_{1}$ represents density-density interaction, and preserves spin SU(2) and valley U(1) symmetry. In fact, $H_{1}$ has an enlarged SU(2)$\times$SU(2) symmetry, which stands for an independent spin rotational symmetry within each valley. Short-range interactions (e.g., atomic scale on-site Hubbard repulsion), which we do not study explicitly, breaks the SU(2)$\times$SU(2) symmetry down to spin SU(2) symmetry.

The Coulomb interaction $V(\boldsymbol{q})$ can be screened by dielectric environment and nearby metallic gates. We assume that TDBG is encapsulated by an insulator (typically boron nitride), and is in the middle to two metallic gates, which generate an infinite series of equally spaced image charges with alternating signs. Under this image charge approximation, the screened Coulomb potential in momentum space is $V(\boldsymbol{q})=2\pi e^{2}\tanh(qd)/(\epsilon q)$, where $\epsilon$ is the dielectric constant of the encapsulating insulator, and $d$ is the vertical distance between the top (bottom) metallic gate and TDBG. We take $d$ to be 50 nm for all calculations presented in the following. The Coulomb interaction energy scale is set by $E_{C}=e^{2}/(\epsilon a_{M})$. At $\theta=1.24^{\circ}$, $a_{M}\approx 11.4$nm, and $E_{C}\approx$12 meV for $\epsilon=10$. Since the typical Coulomb interaction energy scale is comparable to the bandwidth ($\sim$ 10 meV), there is a strong tendency towards symmetry-breaking phases driven by interactions. The system is characterized by almost-flat narrow noninteracting bands with large Coulomb energy, a classic situation for the manifestation of strong correlation physics.

We use Hartree-Fock (HF) approximation and assume that both moiré periodicity and valley U(1) symmetry are preserved, but allow spin polarization, motivated by the experimental evidence of ferromagnetism Shen et al. ; Liu et al. (a); Cao et al. (b) in TDBG. This leads to the following mean-field HF Hamiltonian

[TABLE]

where the quasiparticle energy $E_{\boldsymbol{k},\tau,s}$ includes moiré band energy $\varepsilon_{\boldsymbol{k},\tau}$ and Hartree-Fock self energy $\Sigma_{\boldsymbol{k},\tau,s}$, and $n_{F}$ is the Fermi-Dirac occupation number. By projecting the interaction onto the first conduction band, we neglected self-energy induced by interaction with all other occupied bands. This is expected to be a reasonable approximation for TDBG because of the energetic separation of neighboring bands from the first conduction band. The neighboring band effects could, in principle, be included in the theory, if necessary, either by developing a brute-force multiband mean field theory or perturbatively, with the higher bands contributing paramterically smaller terms because of the large energy denominators associated with band separations.

We denote the electron filling factor as $\nu=n/n_{s}$, where $n$ is the electron density, and $n_{s}$ the density for 4 electrons per moiré cell. At half filling $\nu=1/2$ that corresponds to 2 electrons per moiré cell, the mean-field theory in Eq. (7) leads to two distinct symmetry-broken states, namely, valley-polarized state and valley-unpolarized state, which are degenerate at this particular filling. Short-range interactions can break this symmetry, as explained in Appendix A. Because the moiré conduction bands carry a valley-contrast Chern number, the valley polarized state supports quantum anomalous Hall effect (QAHE). To our knowledge, QAHE has not yet been observed in TDBG. Therefore, we leave the valley polarized state at $\nu=1/2$ in TDBG to future study, and focus on valley unpolarized states in the following.

Because the Hamiltonian $H=H_{0}+H_{1}$ has the enlarged SU(2)$\times$SU(2) symmetry, a valley unpolarized state can have independent spin polarization in the two valleys. For example, the two valleys can have either parallel or antiparallel spin polarization. However, atomic-scale on-site Hubbard interaction explicitly breaks SU(2)$\times$SU(2) down to SU(2) symmetry, and selects the ferromagnetic state in which spins of the two valleys are polarized along the same direction (see Appendix A). In the following, we only consider valley unpolarized state with identical spin polarization in the two valleys. Therefore, we make the ansatz that $\Sigma_{\boldsymbol{k},\tau,s}=\Sigma_{-\boldsymbol{k},-\tau,s}$. This mean-field ansatz preserves spinless time-reversal symmetry.

With the above ansatz, we solve the mean-field theory in Eq. (7) self consistently. One characteristic quantity is the zero-temperature ($T=0$) ferromagnetic gap $\Delta_{\text{FM}}$ that separates the occupied spin majority states from the unoccupied spin minority states at $\nu=1/2$. We plot $\Delta_{\text{FM}}$ as a function of the layer dependent potential $U$ for different dielectric constant $\epsilon$ in Fig. 4, which shows that $\Delta_{\text{FM}}$ has a dome shape and is positive only over a finite range of $U$ for large $\epsilon$. The dome shape in $\Delta_{\text{FM}}$ correlates with the non-interacting DOS [Fig. 2(b)], which peaks at $\nu=1/2$ at $U\approx 41$meV. A larger non-interacting DOS at $\nu=1/2$ implies a stronger instability towards symmetry breaking, and therefore, a larger interaction driven energy gap. However, we note that this argument is only qualitative, as $\Delta_{\text{FM}}$ does not exactly follow the non-interacting DOS. An important conclusion we can draw from Fig. 4 is that the FM insulating gap is tunable by an external displacement field, which agrees with the experimental observation that the FMCI at $\nu=1/2$ only develops over a finite range of displacement field Shen et al. ; Liu et al. (a); Cao et al. (b); Burg et al. (2019); He et al. .

III.2 Magnon Spectrum

A positive FM insulating gap $\Delta_{\text{FM}}$ indicates that the FM state is a good ansatz at the mean-field level. To examine whether the FM state is stable beyond mean-field theory, we calculate the energy spectrum for one-magnon collective excitations. There are two types of magnons for the FM insulator at half filling, namely, spin magnons and valley magnons. The spin magnons involve collective particle-hole transitions from the occupied spin majority band to the unoccupied spin minority band within the same valley, as illustrated in Fig. 5(a); the valley magnons are collective particle-hole transitions that flip the valley index, as shown in Fig. 5(b). We calculate the spin and valley magnon spectrum separately by solving their corresponding Bethe-Salpeter equations, following the theory developed in Ref. Wu and Das Sarma, 2020. Details of the theory can also be found in Appendix B.

Representative spectra for spin and valley magnons are shown in Figs. 5(a) and 5(b), respectively. The spin magnon spectrum has gapless spin wave modes, consistent with the Goldstone’s theorem, as the continuous SU(2)$\times$SU(2) symmetry is spontaneously broken in the ferromagnetic state. In fact, the SU(2) symmetry associated with each valley is broken, so there are two spin wave modes, one for each valley. The overall spin excitation spectrum is nonnegative in Fig. 5(a), showing the stability of the FM insulator against spin-magnon excitations.

By contrast, the valley magnon spectrum shown in Fig. 5(b) is gapped. This is consistent with the fact that there is only U(1) symmetry in the valley space, and the FM insulator does not break this valley U(1) symmetry. The positive valley magnon spectrum indicates the stability of the FM insulator against valley-magnon excitations, and also implies that the FM insulator is energetically more favorable than inter-valley density wave state.

Based on the spectrum shown in Fig. 5, we conclude that the half-filled FMCI in TDBG can be stable against one-magnon collective excitations.

III.3 Mean-Field Transition Temperature

We now turn to finite temperature physics and calculate the mean-field transition temperature $T_{\text{FM}}$ for the FM phase. To determine $T_{\text{FM}}$, we define $\Sigma_{\boldsymbol{k},\tau}^{(\pm)}=(\Sigma_{\boldsymbol{k},\tau,\uparrow}\pm\Sigma_{\boldsymbol{k},\tau,\downarrow})/2$. At $T_{\text{FM}}$, $\Sigma_{\boldsymbol{k},\tau}^{(-)}$ is infinitesimally small, and the self-consistent equation (7) can be linearized as follows

[TABLE]

from which $T_{\text{FM}}$ can be obtained by requiring the largest eigenvalue of the matrix $M^{(\tau)}$ to be 1.

A representative plot of $T_{\text{FM}}$ as a function of filling factor is shown in Fig. 6, where ferromagnetism develops over a large range of filling factors with $T_{\text{FM}}$ up to few tens of kelvin. We note that our mean-field theory overestimates the tendency towards ordering, as fluctuations like spin waves are neglected in the estimation of $T_{\text{FM}}$. The ferromagnetic state at half filling can be an insulator at zero temperature, when the spin majority bands are fully filled and separated from the empty spin minority bands by an energy gap $\Delta_{\text{FM}}$. Away from half filling, the ferromagnetic state is generically metallic with spin dependent Fermi surfaces.

IV Superconductivity

The metallic state away from half filling can be susceptible to superconducting instability due to enhanced electron-phonon interaction in moiré flatband systems. Here we study superconductivity mediated by electron-acoustic phonon interactions. The in-plane acoustic longitudinal phonon modes mediate effective electron attraction as follows

[TABLE]

where $\hat{\psi}_{\tau\sigma\ell s}(\boldsymbol{r})$ is the electron field operator at the coarse-grained position $\boldsymbol{r}$ associated with valley $\tau K$, sublattice $\sigma=A,B$, layer $\ell=1,2,3,4$ and spin $s=\uparrow,\downarrow$. In Eq. (9), we only retain attractive interactions that pair electrons from opposite valleys. The coupling constant $g_{0}$ is given by $D^{2}/(\rho_{m}v_{s}^{2})$, where $D$ is the deformation potential, $\rho_{m}$ is the mass density of monolayer graphene, and $v_{s}$ is the velocity of acoustic longitudinal phonon. Using $D=30$ eV, $\rho_{m}=7.6\times 10^{-8}$ g/cm2, $v_{s}=2\times 10^{6}$ cm/s, we estimate $g_{0}$ to be 474 meV nm2. Here we neglect retardation effects in the phonon mediated electron attraction for simplicity.

As we showed previously in Ref. Wu et al., 2019, the attraction in Eq.(9) can be decomposed into four different pairing channels that are distinguished by their orbital and spin characters: (1) intrasublattice spin-singlet $s$-wave pairing, i.e., $(is_{y})_{ss^{\prime}}\hat{\psi}_{+\sigma\ell s}^{\dagger}\hat{\psi}_{-\sigma\ell s^{\prime}}^{\dagger}$; (2) intersublattice spin-triplet $p$-wave pairing, e.g., $\mathcal{F}_{ss^{\prime}}\hat{\psi}_{+A\ell s}^{\dagger}\hat{\psi}_{-B\ell s^{\prime}}^{\dagger}$, where $\mathcal{F}$ can be any one of the three symmetric tensors $(s_{0}\pm s_{z})/2$ and $s_{x}$; (3) intersublattice spin-singlet $d$-wave pairing, e.g., $(is_{y})_{ss^{\prime}}\hat{\psi}_{+A\ell s}^{\dagger}\hat{\psi}_{-B\ell s^{\prime}}^{\dagger}$; and (4) intrasublattice spin-triplet $f$-wave pairing, i.e., $\mathcal{F}_{ss^{\prime}}\hat{\psi}_{+\sigma\ell s}^{\dagger}\hat{\psi}_{-\sigma\ell s^{\prime}}^{\dagger}$. The $s$-wave and $f$-wave pairings are only distinguished by their spin characters, and the same is true for $p$ and $d$ pairings. The angular momenta of intersublattice Cooper pairs arise from the valley-contrast sublattice chirality under $C_{3}$ rotation Wu et al. (2018, 2019). In $AB$ bilayer graphene, one of the sublattices in each layer [$A_{1}$ and $B_{2}$ sites in Fig. 1(b)] is pushed to higher energy by interlayer tunneling. Therefore, intersublattice pairing is energetically less favorable compared to intrasublattice pairing in TDBG. In the following, we only consider interactions that pair electrons on the same sublattice, and project such interactions onto the first moiré conduction band. The projected pairing Hamiltonian is

[TABLE]

where we only keep interactions that pair electrons with opposite momenta, i.e., momentum $\boldsymbol{k}$ in $+K$ valley and momentum $-\boldsymbol{k}$ in $-K$ valley. The pairing Hamiltonian $H_{\text{p}}$ also has the SU(2)$\times$SU(2) symmetry, and supports both spin singlet $s$-wave and spin triplet $f$-wave pairings. Because of the ferromagnetism induced by Coulomb repulsion, equal spin pairing is more favored compared to spin singlet pairing. Therefore, we consider intervalley pairing between electrons with the same spin, which leads to the following Bardeen-Cooper-Schrieffer (BCS) mean-field Hamiltonian

[TABLE]

By combining the BCS Hamiltonian $H_{\text{BCS}}$ and the effective single-particle Hamiltonian $H_{\text{FM}}$ [Eq. (7)] that is renormalized by the Coulomb interaction, we obtain the superconducting linearized gap equation

[TABLE]

where $\mu$ is the chemical potential, and $E_{\boldsymbol{k},\tau,s}=\varepsilon_{\boldsymbol{k},\tau}+\Sigma_{\boldsymbol{k},\tau,s}$ is the effective band energy including the self energy. We have used the spinless time-reversal symmetry, which implies $E_{\boldsymbol{k},\tau,s}=E_{-\boldsymbol{k},-\tau,s}$, to simplify the superconducting susceptibility $\chi^{(s)}$. Because of this symmetry, ferromagnetism does not lead to depairing effect for superconductivity with intervalley equal-spin pairing. In Eq. (12), spin up and down channels have independent gap equations. The superconducting transition temperature $T_{\text{SC}}$ is reached when the largest eigenvalue of $\chi^{(s)}$ is 1. Fig. 6 plots $T_{\text{SC}}$ as a function of filling factor, and shows two superconducting domes respectively on the two sides of the half-filled ferromagnetic state. In Fig. 6, we take the value of $g_{0}$ (the attractive interaction strength) to be three times of 474 meV nm2 (the value obtained from the above electron-acoustic phonon coupling parameters) in order to get a value of $T_{\text{SC}}$ on the order of 1 K. We note that $T_{\text{SC}}$ is exponentially sensitive to $g_{0}$ as well as the moiré band flatness. A quantitative study of $T_{\text{SC}}$ is beyond the scope of this paper. In any case, the experimental parameters are not known with sufficient accuracy for a quantitative estimate of $T_{\text{SC}}$ at this stage of development of the field. The main purpose of this section is to point out the possibility of phonon-mediated spin triplet pairing in a ferromagnetic system.

We discuss the effect of an in-plane magnetic field $\boldsymbol{B}_{\parallel}$ on $T_{\text{SC}}$ in the $f$-wave channel. If the parent state for superconductivity is spin unpolarized, then $T_{\text{SC}}$ can be slightly enhanced by $\boldsymbol{B}_{\parallel}$ in the low-field regime, because Zeeman energy leads to an effective spin dependent chemical potential shift Lee et al. ; Wu and Das Sarma ; Scheurer and Samajdar . On the other hand, if the parent state already has maximum spin polarization allowed by a given filling factor, then an externally applied $\boldsymbol{B}_{\parallel}$ field can no longer change the amount of spin polarization, and $T_{\text{SC}}$ is reduced by $\boldsymbol{B}_{\parallel}$ due to orbital effect Lee et al. ; Wu and Das Sarma . In Ref.Liu et al., a, $T_{\text{SC}}$ is found to be slightly enhanced by weak $\boldsymbol{B}_{\parallel}$ field, indicating that the superconducting state has spin triplet pairing but with no spin polarization. Our mean-field phase diagram in Fig. 6 likely overestimates the filling range for ferromagnetism. We emphasize that there is always a superconducting instability in a partially filled band regardless of the presence or absence of ferromagnetism in our theory, where the superconductivity is mediated by electron-phonon interactions and ferromagnetism is driven by Coulomb repulsion. An additional signature of electron-acoustic phonon interaction in TDBG is that phonon scattering can lead to large linear-in-$T$ resistivity in transport above some crossover temperatures Li et al. .

Finally, We note that experimental signatures of SC in TDBG are not yet conclusive, as discussed in detail in Ref. He et al., .

V Conclusion

In conclusion, we have presented a theory of ferromagnetism induced by Coulomb repulsion and superconductivity mediated by electron-acoustic phonon interactions in moiré bands of TDBG. In our theoretical phase diagram, there can be a ferromagnetic correlated insulator at half filling, and superconducting domes on both the electron and hole sides of the half-filled insulator. Ferromagnetism and superconductivity are two prototypical orders that can occur in moiré flat bands, while there are many other possible competing and/or intertwined orders, such as nematicity that breaks rotational symmetry and density wave state that breaks moiré translation symmetry Hsu et al. . In TDBG, there is experimental evidence that states with both spin and valley polarization are possibly stabilized at 1/4 and 3/4 fillings by a finite in-plane magnetic field Liu et al. (a); Cao et al. (b); He et al. . CIs at these factors could also be spin and/or valley polarized states. Because of the valley contrast Chern numbers in the non-interacting moiré bands, valley polarized CIs can also display quantum anomalous Hall effects. Our work should be viewed as a step towards a full quantitative theory of the potentially very rich TDBG phase diagram. A note-worthy qualitative feature of the current work is the possibility, already apparent at the mean field level, that SC and FMCI phases, although they arise from different interactions (electron-phonon for SC and electron-electron for FMCI), could compete with each other in TDBG moiré flatband with the FM phase centered around half-filling and the SC domes manifesting on both electron- and hole-doped sides of half-filling. The fact that this could be the experimental TDBG situation may indicate that our theory captures some essential qualitative aspect of moiré interaction physics although our use of mean field theory (and many other approximations, e.g., neglect of higher bands) exaggerates the quantitative stability of the symmetry-broken phases compared with experiments.

TDBG and other related moiré systems represent a highly tunable platform, where moiré band structure can be effectively controlled by the out-of-plane displacement field, as revealed by our theoretical study. This feature allows in situ control of band structure, and provides unprecedented opportunities to study many-body physics.

VI acknowledgment

F. W. thanks Y.-T. Hsu, X. Li, and R.-X. Zhang for discussions. This work is supported by Laboratory for Physical Sciences.

Appendix A Short-Range Interactions

We show that short-range interactions, in particular, the atomic scale on-site Hubbard repulsion, explicitly break the SU(2)$\times$SU(2) symmetry down to spin SU(2) symmetry, and favors ferromagnetic states in which spins in the two valleys are polarized along the same direction.

The on-site Hubbard repulsion on monolayer graphene honeycomb lattice is given by

[TABLE]

where $U_{0}$ is the on-site Hubbard interaction, $\boldsymbol{R}$ is the lattice vector, $\sigma$ is the sublattice index ($A$, $B$), $N$ is the number of unit cells in the monolayer, and the prime on the summation of the second line is the momentum conservation constraint, i.e., $\boldsymbol{k}_{1}+\boldsymbol{k}_{2}$ is equivalent to $\boldsymbol{k}_{3}+\boldsymbol{k}_{4}$ modulo reciprocal lattice vectors. To obtain a continuum model, we only keep states near $\pm K$ points:

[TABLE]

where $\tau=\pm 1$ is the valley index, and the prime on the summation of $\tau_{i}$ implies the valley conservation $\tau_{1}+\tau_{2}=\tau_{3}+\tau_{4}$ due to momentum conservation. In the operator $b^{\dagger}_{\boldsymbol{k}\tau\sigma s}$, the momentum $\boldsymbol{k}$ is measured relative to $\tau K$. We make a Fourier transformation to introduce the coarse-grained real-space position $\boldsymbol{r}$:

[TABLE]

where $\mathcal{A}$ is the system area. The on-site repulsion can then be transformed to a continuum Hamiltonian with local interaction:

[TABLE]

where $u_{0}=U_{0}\mathcal{A}_{0}$, and $\mathcal{A}_{0}=\sqrt{3}a_{0}^{2}/2$ is the area per unit cell in the monolayer. The local repulsion in Eq. (16) can swap the valley indices of a pair of electrons, and therefore, break the SU(2)$\times$SU(2) symmetry down spin SU(2) symmetry.

We project $H_{2}$ in Eq. (16) to the first moiré conduction band, perform Hartree-Fock decomposition using the ansatz given in the main text, and obtain the following mean-field Hamiltonian:

[TABLE]

[TABLE]

It is clear from Eq. (17) that the local repulsion favors a valley unpolarized but spin polarized ferromagnetic state with spins associated with the two valley polarized to the same direction.

The energy scale for the short-range repulsion is $E_{U}=U_{0}a_{0}^{2}/a_{M}^{2}$, which is about 2 meV using $U_{0}=5$ eV, and is an order-of-magnitude weaker compared to the long-range Coulomb interaction. The atomic-scale on-site Hubbard interaction acts a weak anisotropy that breaks SU(2)$\times$SU(2) down to SU(2) symmetry, and selects a particular set of states out of an SU(2)$\times$SU(2) multiplet. For examples, the short-range Hubbard interaction aligns spins associated with the two valleys in the ferromagnetic phase, and suppresses $s$-wave but not $f$-wave pairing in the superconducting phase.

Appendix B Theory for magnon excitations

In this Appendix, we present the theory for magnon excitations, which has been discussed in moiré systems in Refs. Wu and Das Sarma, 2020 and Alavirad and Sau, . The spin magnon states can be parametrized as follows

[TABLE]

where $|\text{FM}\rangle$ is the half-filled FM insulating state in which spin $\downarrow$ bands in both valleys are empty, $z_{\boldsymbol{k},\boldsymbol{Q}}$ are variational parameters, and $\boldsymbol{Q}$ is the momentum of the magnon. In the magnon state $|\boldsymbol{Q}\rangle_{S}$, we make a single spin flip from the occupied spin $\uparrow$ band to unoccupied spin $\downarrow$ band within the same $+K$ valley. Variation of the magnon energy with respect to $z_{\boldsymbol{k},\boldsymbol{Q}}$ leads to the following Bethe-Salpeter equation

[TABLE]

where the first part in $\mathcal{H}_{\boldsymbol{k}\boldsymbol{k}^{\prime}}^{(\boldsymbol{Q})}$ is the quasiparticle energy cost of the particle-hole transition, and the second part represents the electron-hole attraction. The eigenvalue $\mathcal{E}_{S}(\boldsymbol{Q})$ represents the energy of spin magnons. We note that there is another spin wave mode in $-K$ valley, which can be formulated in a similar way as Eq. (19).

In addition to spin magnon states, there are also valley magnon states with a valley flip

[TABLE]

The corresponding Bethe-Salpeter equation is

[TABLE]

which gives rise to the valley magnon spectrum in Fig. 5(b).