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

###### 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 double bilayer graphene (TDBLG), consisting of 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 . We construct a phase diagram as a function of twist angle and , incorporating interactions via a Hartree-Fock approximation. At half filling, ferromagnetic insulators are stabilized, typically with valley Chern number . 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.

^{†}

^{†}J.Y. Lee and E. Khalaf contribute equally to this work.

## I Introduction

Recent experiments have reported interaction induced insulating states and superconductivity in twisted bilayer graphene and related systems Cao *et al.* (2018a, b); Chen *et al.* (2019a, b); Yankowitz *et al.* (2019); Lu *et al.* (2019).
In all these settings, long wavelength Moiré patterns are believed to give rise to narrow electronic bands Lopes dos Santos *et al.* (2012); Bistritzer and MacDonald (2011a); Morell *et al.* (2010); Chen *et al.* (2019a); Zhang *et al.* (2019) that enhance the effect of electron-electron interactions. Indeed the observed phase diagram of correlated insulators with proximate superconducting states is reminiscent of strongly correlated materials such as the cuprates Lee *et al.* (2006). On the other hand, a different route to correlated insulators is observed in quantum Hall systems, for instance, when the spin and valley degeneracy of graphene’s Landau levels are spontaneously broken by interactions, leading to ferromagnetic insulators. In contrast, ferromagnetic insulators are relatively rare in correlated solids, where antiferromagnetic order is the norm. A key question in the field of Moiré materials is which paradigm is better suited to capturing the physics of TBLG. Topological aspects of the nearly flat bands of TBLG have recently been emphasized Po *et al.* (2018a, b); Ahn *et al.* (2018); Lian *et al.* (2018a) along with the connection to quantum Hall wavefunctions Tarnopolsky *et al.* (2019); Bultinck *et al.* (2019), although time reversal remains a symmetry. As a consequence of these complexities, and despite much theoretical work Xu and Balents (2018); Po *et al.* (2018a); Isobe *et al.* (2018); You and Vishwanath (2018); Kang and Vafek (2018); Xie and MacDonald (2018); Lin and Nandkishore (2019); Dodaro *et al.* (2018); Padhi *et al.* (2018); Wu *et al.* (2018); Lian *et al.* (2018b); Zou *et al.* (2018), the origin and nature of the insulators and superconductors in TBLG remains disputed.

Here we study the problem of a pair of bilayer graphene sheets, twisted with respect to one another with AB-AB stacking structure. The essential difference from the monolayer on monolayer problem is a difference in symmetry. The rotation symmetry of TBLG about the vertical axis is broken by the Bernal stacking of the individual bi-layers in twisted double BLG (TDBLG). This lower symmetry setting simplifies the problem in some ways - for example in TBLG, the conduction and valence bands touch at Dirac points. In contrast, these bands can be separated once is broken. Furthermore, the bands can now acquire a net Chern number (opposite valleys have opposite Chern numbers owing to time reversal symmetry). This allows us to draw intuition from the related and well studied problem of electrons in a Landau level. For example, in graphene, an appropriate strain pattern produces a uniform ‘psuedo’-magnetic field with opposite orientation for the two valleys Levy *et al.* (2010). This problem, in the absence of spin degeneracy, was numerically studied in Ghaemi *et al.* (2012). At neutrality, when the lowest Landau levels for the two valleys are half filled, time reversal is spontaneously broken leading to a valley ferromagnet - which is also an integer quantum Hall insulator. Indeed a related ground state with spontaneous quantum Hall response, although metallic, was observed in TBLG in the presence of breaking substrate potentials Sharpe *et al.* (2019); Xie and MacDonald (2018); Zou *et al.* (2018); Zhang *et al.* (2019); Bultinck *et al.* (2019). Thus, ferromagnetic ground states are to be expected in this regime.

In this work we will study twisted double BLG as follows: (i) We model the single particle physics of TDBLG, including the effects of additional terms like trigonal warping which play little role in the monolayer case, but significantly modify the dispersions here. A finite vertical electric field is then required to obtain isolated and relatively flat bands, and a parameter regime where this occurs for the conduction band is identified and its topological properties are computed. (ii) Next, we study the effect of interactions within a self consistent Hartree Fock approximation. At half-filling a spin ferromagnetic insulator is predicted, which is assisted by the Hund’s interaction as anticipated in Zhang *et al.* (2019) and consistent with the numerical studies of graphene in a psuedomagnetic field Ghaemi *et al.* (2012) (iii) We study possible instabilities of the doped ferromagnet, in particular the possibility of pairing induced by ferromagnetic fluctuations, within a simplified model. A spin triplet, inter-valley paired superconductor with uniform pairing amplitude on the Fermi surface is predicted (a spin triplet ‘s-wave’, where antisymmetry is enforced by valley quantum numbers). (iv) Finally, we study experimental signatures of these states in a magnetic field, applied either parallel or perpendicular to the plane, which compare well with a recent experiment Liu *et al.* (2019). In particular, we identify an unexpected orbital effect of an in-plane magnetic field, which competes with the Zeeman splitting.

## Ii Single-particle physics

### ii.1 Hamiltonian and symmetries

We consider a system consisting of two layers of AB-stacked bilayer graphene as follows. Starting from an initial configuration with ABAB stacking, we twist the two bilayers relative to each other by a small angle . Each bilayer graphene (BLG) has a tight-binding Hamiltonian given by

(1) |

which is labelled in the order of , , , . Here, we consider a realistic model of BLG. AB-stacking means that the -site of the first layer () sits on top of the -site of the second layer (). This gives the small on-site energy for these sites. For a detailed explanation, see Appendix. A. Here, , where , , . One can expand near as

(2) |

Here, is the distance between carbon atoms. Throughout, we will use the phenomenological value extracted from Ref. Jung and MacDonald (2014)

(3) |

Additionally, the potential difference between the top and bottom graphene layer, is an important parameter in the experiment, which is controlled by the gate voltage difference. For a displacement field strength , ABAB system’s dielectric constant and the thickness of the BLG/BLG system , .

The bottom layer of the top BLG and the top layer of the bottom BLG are coupled via the Moiré hoping term, given by

(4) |

where the hopping matrix couples the Bloch states of the top and bottom BLG with the momentum difference , where . In the original Bistritzer-Macdonald model, Bistritzer and MacDonald (2011b). However, in a realistic twisted model, the ratio is not equal to one which accounts for the effects of lattice relaxation. This is shown to be important for obtaining the observed bandgaps in twisted single BLG, where is taken to be around 0.75 for the first magic angle Nam and Koshino (2017); Koshino *et al.* (2018).
In the case of twisted single BLG, the value of is crucial for the gap between first and second conduction (valence) bands. This is also true for the twisted double BLG in our model; in fact, the model without relaxation cannot explain the insulating state at the filling state in twisted double BLG Liu *et al.* (2019). Indeed, a recent study by Choi et al. Choi and Choi (2019) showed significant energy gap between first and second conduction bands at , which can only be explained in a continuum model with relaxation.
Although it is reported that the relaxation parameter is a monotonically decreasing function of the twist angle Carr *et al.* (2019), we will use the fixed value of for entire twist angles. Here, we choose , meV and meV for the plots presented in this paper. Note that the qualitative features do not change even if we enhance the relaxation. For numerical results with different relaxation parameters, see Appendix. B.

In an idealized model, the band becomes almost perfectly flat at the angle as in Fig. 2(a). However, in a realistic model we include trigonal warping () and particle-hole asymmetry () within each BLG which cannot be ignored, hence the numerical results obtained here differs from those in Ref. Zhang *et al.* (2019), as shown in Fig. 2(b). In untwisted systems, these terms are unimportant because their energy scale is a tiny fraction of the relevant band width. However, in the Moiré band of the twisted systems, these terms can drastically change the physics, giving rise to a noticeable effect on
band widths or separations in real experiments Liu *et al.* (2019). Once these terms are included, the bands acquire a sizable dispersion (10-20 meV) and we do not obtain the vanishing bandwidth behavior (‘magic angle’) anymore.

The Hamiltonian (71) has the following symmetries (i) three-fold rotation symmetry , and (ii) time-reversal symmetry . Furthermore, in the absence of a vertical electric field one obtains (iii) mirror reflection about the -axis , The first and third act within each valley while (ii) exchanges valleys. Mirror symmetry also permutes A and B sublattice sites and exchanges top and bottom layers, and . A non-zero perpendicular electric field breaks since . Since reverse the sign of Berry curvature, the Chern number of each band satisfies . In the absence of a vertical electric field, this implies . Note that time reversal symmetry flips momentum and exchange valleys at the same time. Therefore, symmetries constrain band dispersions by the following:

(5) |

Given these symmetries, we only present the band structure along the specific cut –––. When , Mirror symmetry maps the system backs to itself while flipping . Thus, and are degenerate along when . Finally, we have assumed that in the small angle limit, valley quantum number is well defined implying a symmetry arising from the decoupling of Moiré and atomic lattice scale physics.

We note here a crucial difference between twisted double bilayer graphene compared to its single bilayer counterpart which is the absence of two-fold rotational symmetry (i.e. vs rotation symmetry). This additional symmetry when combined with time reversal leads to a symmetry that leaves the momentum and valley index invariant. This combined symmetry ensures vanishing Berry curvature and Chern number and is responsible for the existence of a Dirac point which cannot be removed in the TBLG system, which implies that the minimal set of isolated bands at neutrality consists of two bands per spin per valley. In contrast, the physics of the double bilayer system is controlled by a single narrow band (per spin per valley).

### ii.2 Band structure

Unlike the twisted monolayer-monolayer system, Bistritzer and MacDonald (2011b), the twisted bilayer-bilayer system exhibits a significant overlap of the conduction and valence band due to the trigonal warping () and particle-hole asymmetry () terms. (See Fig. 2). Such band overlap would tend to reduce interaction effects; thus, one needs to tune the system so that the band of interest is isolated and relatively flat. This can be achieved by applying a gate voltage between top and bottom layers, which introduces a vertical electric field that separates the two bands. Using numerical simulation, we checked the conditions under which such band isolation is possible. We found that this occurs on tuning the displacement field, but only for the first conduction band. This behavior is originated from the particle-hole asymmetry term McCann and Koshino (2013), which significantly reduces the isolated region for the valence band. To characterize the isolation of the first conduction band, we define

(6) |

where represents a dispersion of -th conduction/valence band.

In Fig. 3, we plot the region where the conduction band is isolated. We notice that although the bandwidth is not as flat as that reported in numerical studies of magic-angle twisted bilayer graphene, correlation effects can still play an important role. Indeed, there is some debate regarding the bandwidth of magic angle TBLG itself, with reported bandwidths ranging from 10-40 meV Choi *et al.* (2019). In the next section, we will demonstrate that the twisted double BLG system can still host a correlation induced insulating phase by using a self-consistent Hartree-Fock calculation.

Indeed, consistent with these observations, the recent experimental data on twisted double BLG Liu *et al.* (2019) shows 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 Liu *et al.* (2019), consistent with our observation regarding band isolation and band width of the conduction versus the valence band.

### ii.3 Chern Number and Wannier Obstruction

One of the interesting features in Moiré systems is the appearance of the non-zero Chern number . In ordinary condensed matter systems, time-reversal forbids the existence of Chern bands. However, in Moiré systems, Chern bands carrying opposite Chern numbers for different valleys can arise at a finite gating voltage due to the valley decoupling. The overall system still satisfies time reversal symmetry which exchanges the two valleys. Therefore, if the spontaneous valley polarization occurs, a Chern insulator can appear without explicitly breaking time reversal symmetry Zhang *et al.* (2019); Ghaemi *et al.* (2012); Xie and MacDonald (2018); Bultinck *et al.* (2019).

First, at , the reflection symmetry enforces for both valleys. ( maps the system back to itself without exchanging valleys, but so Berry curvature flips its sign),
In the idealized case Zhang *et al.* (2019), as we increase , band inversion between conduction and valence bands occurs at the Moiré -point ( for negative ) with a quadratic touching. Thus, Chern number of is exchanged.

This feature persists in the realistic Hamiltonian of Eq. 71.
The following scenario for the development of Chern band applies for most regimes in Fig. 4. In the non-idealized case with trigonal warping, 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 Moiré hopping, the degeneracy between four Dirac cones split, and the band inversion would happen first at three Dirac cones, exchanging Chern number by . Then, the band inversion would occur at the center Dirac cone, exchanging Chern number by . In total, it will still change the net Chern number by . At larger values of the gate voltage , the band inversion happens between first and second conduction band at point, and the Chern number then changes by ^{2}^{2}2It can change by for other parameter setting, decreasing the Chern number.

This can be further checked numerically by inspecting symmetry indicators of Chern insulators. There are three -invarinat momenta in the Brillouin zone, , , and . For a Bloch state with these momenta, rotation symmetry would map the state back to itself with a rotation eigenvalue:

(7) |

where is an angular momentum associated with the Bloch state .
Then, the Chern number of the th band can be determined modulo 3 by Fang *et al.* (2012); Turner *et al.* (2012); Matsugatani *et al.* (2018)

(8) |

Thus, by tracking how eigenvalues of the three momenta change with the gating voltage , we can understand how Chern number transition happens in the system. Indeed, the aforementioned scenario can be confirmed. For example, consider a Moiré first conduction band for valley at . At meV, we start with . At meV, Chern number changes by but it can be only captured by Berry curvature not by symmetry indicator. At meV, Chern number changes by , manifested by . At meV, Chern number again changes by , manifested by . See Fig. 11 in the Appendix for the detail.

### ii.4 Effect of magnetic field

Applying a magnetic field can give rise to two separate effects. The first is the Zeeman effect which splits bands with opposite spins by the amount with the Bohr magneton . The Zeeman splitting is independent on the direction of the applied field provided that the spins can align with it (there is no strong anisotropy or spin-orbit coupling). The second effect results from the coupling of the magnetic field to the orbital motion of the electron. The orbital effect is significantly larger for out-of-plane magnetic fields, but as we will discuss below, it cannot be neglected for in-plane fields. In this case, the magnetic field slightly alters the bands structure and can have subtle effects on physical observable.

#### ii.4.1 In-plane orbital effect

Let us now consider the effect of in-plane magnetic field . We can choose the gauge which does not depend on or , thus preserving the Moiré translation symmetry. The effect of a magnetic field on hopping terms can be easily included via Peierl’s substitution, where the hopping term from to is multiplied by the phase factor

(9) |

such that

(10) |

where

(11) |

where is the matrix element connecting layers and ( from bottom to top), is the interlayer distance, and is the unit vector in the direction.

Due to its small magnitude relative to the energy gap, it suffices to consider the in-plane orbital effect to first order in pertrubation theory. This amounts to adding the following in-plane orbital term to the single particle energies

(12) |

where is given by

(13) |

where is the valley index. Time-reversal symmetry implies that .

The orbital effect due to in-plane field amounts to a very small relative momentum shift and naively one might think it should be neglected. However, when compared to Zeeman effect, we find that they are of the same order of magnitude since (here is the Fermi velocity for single-layer graphene ). As a result, it cannot always be safely neglected as we will show later in this work.

The in-plane orbital -factor transforms under rotation as

(14) |

provided that the band is non-degenerate at . This implies that vanishes at any -invariant point.

In general, the in-plane orbital contributions affects the bands very differently from the Zeeman effect due to the dependence of the -factor. More specifically, the in-plane orbital effect changes the dispersion without moving the bands at high symmetry points whereas the Zeeman effect shifts the entire band uniformly. If the band minima/maxima occur at high symmetry point, then the in-plane orbital effect will have no influence (at least to linear order in field) on the energy gaps. However, 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 , as we will show later.

The in-plane orbital -factor for the isolated band is shown in Fig. 5(a). We can see that its magnitude can exceed the Zeeman effect in some areas in the Brillouin zone. We can also see that it has vanishing magnitude at the , and points which are -invariant.

#### ii.4.2 Out-of-plane orbital effect

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. Koshino (2011); Ju *et al.* (2017) whose -factor is given by

(15) |

Due to the time reversal symmetry, . The out-of-plane orbital -factor for the isolated band is given in Fig. 5. As we can see its value can be as large as 100 times the Zeeman effect. Due to the structure of the band dispersion in Fig. 5, one can control the bandwidth by applying the out-of-plane -field.

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

(16) |

where is the electron spin operator (which is for up/down spins) and . The valley orbital -factor is defined as

(17) |

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

## Iii Hartree-Fock mean field

Comparing Fig. 3 and 4, one can notice that the first conduction band carries a non-zero Chern number at most of the isolated regime. It is well known that non-zero Chern numbers prevent the existence of exponentially localized Wannier function Marzari and Vanderbilt (1997), and one cannot construct a Hubbard model unless valley-symmetry is broken or the number of bands in this model is enlarged to have the net Chern number zero. As a result, we cannot have a simple description for the system if we confine ourselves only to the first conduction band. Therefore, instead of seeking a real-space description, we discuss the many-body physics in -space Guinea and Walet (2018). As a consequence, it is impossible to state whether the emergent insulating phases would exhibit Mott-like physics arising from strong real-space repulsion. Thus, we will use the terminology, *correlated insulator* to refer to the insulating phase originating from the interaction.

In this section, we perform a self-consistent Hartree-Fock mean field theory to uncover the nature of the possible correlated insulating states at half and quarter filling observed in Ref. Liu *et al.* (2019). Our approach is similar to the one employed in Ref. Po *et al.* (2018a) with the main difference being the inclusion of form factors which, as we will show, significantly alters the conclusions.

### iii.1 Projecting the interaction onto the isolated band

The first step in our mean field calculation is to project the Coulomb interaction onto the relevant bands. To a first approximation, the valleys are completely decoupled leading to independent SU(2) rotation symmetry in each valley. Note that, unlike the quantum Hall effect in graphene, which has an effective SU(4) symmetry, here the opposite Chern number of the two valleys implies a lower SU(2) SU(2) SO(4) symmetry. A small intervalley Hund’s coupling term breaks this SU(2) SU(2) down to SU(2) U(1), the microscopic spin rotation symmetry. Let us sketch the derivation from the microscopic theory, with most of the details relegated to Appendix C.

The interacting Hamiltonian in momentum space is given by

(18) |

where and . is taken to be the screened Coulomb potential

(19) |

where denotes the inverse screening length and is the projected density on the isolated band. The summation over in (18) is not limited within the first Moiré BZ. We note that the density operator consists of two pieces: an intravalley part and an intervalley part . The latter only contributes when the valleys are completely decoupled and, by considering it, we can derive the intervalley Hund’s term.

The density is generally non-periodic in 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 is periodic under any reciprocal lattice translation for the original system. As a result, consists of several identical narrow peaks centered at reciprocal lattice vectors of the original bilayer graphene for the intravalley density or at for the intervalley density . This poses a problem since it implies that the summation over in (18) diverges.

To resolve this issue, we notice that the periodicity of 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 will decay for momenta larger than a certain cutoff which is given by the inverse size of the Wannier orbitals. Rather than attempting to precisely determine the value of from the graphene Wannier orbitals, we will consider 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 (18) to the vicinity of for the intravalley density and the vicinity of and for the intervalley density .

The resulting Hamiltonian consists of two parts

(20) |

contains the coupling between intravalley densities whereas contains the coupling between intervalley densities . They are given explicitly by

(21) |

(22) |

where the intravalley and intervalley form factors are defined as

(23) |

All momenta in (21) and (22) are measured in units of with denoting the dimensionless screened Coulomb interaction with a screening length . 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 . Rough estimations for and provide the scale of the two interaction terms and are given by

(24) |

Here, we used and used to denote the value of in degrees. Using a value of of about 5 at twist angles around yields meV and meV. We see that the term is significantly smaller than the term. It can be important, however, since it breaks the symmetry from SU(2) SU(2) to SU(2) 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 interaction. The term generally has the effect of favoring spin alignment and can be written in the form of inter-valley Hund’s coupling as in Zhang *et al.* (2019).

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

(25) | |||

(26) |

Time-reversal symmetry imposes an additional constraint on the gauge transformation, .

### iii.2 General setup: self-consistent Hartree-Fock

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

(27) |

which we will assume to be diagonal in and , . In the following, we will introduce the combined index such that is a matrix with components . Next, we expand the interaction (20) in the difference and neglect terms beyond linear order.

The resulting mean field Hamiltonian has the form

(28) |

Here, is a column vector in the index , is a diagonal matrix containing the single particle energies and are 4 4 matrices in given by

(29) |

and

(30) |

The matrix simply contains the form factors defined in (82)

(31) |

and is the projector on the valley with denoting the Pauli matrices in the valley space.

In both (29) and (30), 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 Po *et al.* (2018a); Zhang *et al.* (2019) 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 . 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 Xie and MacDonald (2018). In that work, the Hartree-Fock corrections to the flat bands coming from all other () 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 Jung and MacDonald (2014); Carr *et al.* (2019). Thus, we only include the effects arising from filling the isolated band.

To write the self-consistency condition, we diagonalize by introducing the variables for some unitary . We then impose the constraint . In the following, we will only consider possible gapped phases at integer fillings . In this case, the self-consistency condition has the form

(32) |

where is a -independent matrix containing ones along the diagonal and zeroes everywhere else. Our assumption that the phase is gapped has to be checked self-consistently by computing the mean field band structure

(33) |

and ensuring that correlation induced gap for filling defined as

(34) |

is positive. Here, we assumed that the mean field bands are sorted in order of increasing energy. ()

We notice that is, in general, not gauge invariant. Instead it transforms as

(35) |

under the gauge transformation (26). In the following mean field analysis, we will choose the gauge such that which guarantees the gauge independence of .

### iii.3 Gapped symmetry-broken phases

In this section, we will investigate the different possible symmetry-broken gapped phases. The model has spatial and translation symmetries which we will assume to be unbroken. In addition, there is spinless time-reversal symmetry which may be spontaneously broken. Without the dispersion, the form factors or the Hund’s coupling term, the model has full SU(4) symmetry corresponding to an arbitrary rotation in spin and valley spaces. When we include the dispersion and form factors, this symmetry is broken down to SU(2) SU(2) corresponding to independent spin rotation in each valley. This symmetry is further broken by the weak intervalley Hund’s coupling term (22) to SU(2) U(1) which correspond to total spin rotation and valley charge conservation.

The structure of the order parameter in spin and valley spaces can be investigated as follows. There are 15 possible order parameters corresponding to the generators of SU(4). Such generators can generally by written as where and represent the Pauli matrices in spin and valley spaces respectively. This splits the generators into the following groups: (i) one generator corresponds to a valley polarized (VP) state which break time-reversal symmetry, (ii) 3 generators correspond to a spin-polarized (SP) state which breaks SU(2) spin symmetry, (iii) 2 generators correspond to inter-valley coherent order (IVC) which breaks valley U(1) charge conservation, (iv) 3 generators correspond to both spin and valley polarization breaking both SU(2) and time-reversal, and (v) 6 generators break both valley U(1) and SU(2) spin rotation, corresponding spin-polarized inter-valley coherent order.

#### iii.3.1 Numerical solution to self-consistency equation

At half-filling, there are three possible translationally symmetric gapped states given by spin polarization, valley polarization or intervalley coherent order which break SU(2) spin rotation, time-reversal, or U(1) valley charge conservation respectively. Relegating the details of the analysis to Appendix D, we summarize the results below. In the following, we will neglect the effect of the intervalley Hund’s coupling since it is much smaller than the main part of the interaction . We will discuss its effect later in III.3.3.

The self-consistency condition for spin or valley polarization can be satisfied by simply taking and which is diagonal in valley and spin spaces and -independent. This means that the self-consistency equation (32) is trivially satisfied (since ). Without the Hund’s coupling, the total energies of the two states are equal. For inter-valley coherent order, the self-consistency condition cannot be satisfied by taking to be -independent. Instead, we need to consider a more general form for the order parameter given by

(36) |

does not break time-reversal symmetry provided that and . The self-consistency equation for the IVC order has a relatively complicated form (provided in Appendix D) and can be only solved numerically. The total ground state energy can then be computed and compared to that of the VP/SP state. The results are shown in Fig. 6 for different values of the potential and the twist angle . From the Figure, we can see that for most of the phase diagram the VP/SP state has a lower energy than the IVC state. This is particularly pronounced in the band isolation region.

The ground state energy for these symmetry-broken phases is only meaningful in our setting if the bands given in (33) are gapped at half-filling. The value of the correlated gap (34) for the VP/SP state is shown in Fig. 6. We can see that the gap is positive everywhere implying that an insulating phase is self-consistent within our analysis. For the band isolation region, we find a gap around 10 meV which depends very weakly on the twist angle and voltage .

At quarter-filling, a gapped insulating state has to be spin-polarized in addition to being either valley polarized or IVC. The competition between the latter two can be investigated similar to the half-filled case by solving the self-consistency condition and comparing the ground state energies. The problem is identical to the half-filed case except that the electrons are now effectively spinless. The number of spin species enters the mean field Hamiltonian (28) through the traces. The Hartree term in the kinetic part and the Fock term in the potential part have a single trace leading to a factor of , whereas the Hartree term in the potential part has two traces leading to a factor of . Thus, the results for quarter filling are obtained from the half-filling results by dividing the Hartree kinetic energy and Fock potential energy by 2 and the Hartree potential energy by 4 and we find that, similar to half-filling, VP is more energetically favorable than IVC for most of the parameter space.

The correlated gap between the split bands in the mean-field theory arises from the Fock term. One subtlety is that the Hartree term affect the band gap as well because it modifies the band dispersion, rendering the band minimum/maximum accordingly. As a result, the correlated gap would *increases* for SP+VP state compared to that of VP state.

#### iii.3.2 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

(37) |

where . This approximation can be justified as follows: the starting symmetry of the isolated band is SU(4) which is broken to SU(2) SU(2) due to the asymmetry between the two valley in energies and form factors (, ). In the following, we will assume that breaking SU(4) to SU(2) 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 . 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 and . This is equivalent to expanding in time-reversal symmetry breaking terms within each valley.

The variables and can be obtained by solving a linearized version of the self-consistency equation (see Appendix D for details). The resulting energy difference between the IVC state and the SP/VP state can be written (up to second order in ) is given by

(38) |

where and are given in Appendix D. The first term in (38) reproduces the non self-consistent Hartree-Fock energies obtained in Ref. Zhang *et al.* (2019) 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. Po *et al.* (2018a) which considers the simplified setting where all form factors are taken equal to 1. In addition, was taken equal to a constant which is cutoff at large momenta yielding the interaction strength and leading to . In this case,

(39) |

which implies that the IVC phase is energetically favored to the VP/SP phase in agreement with the conclusion of Ref. Po *et al.* (2018a) ^{4}^{4}4The result differs by a factor of 2 due to the incorrect way the large limit was implemented Po *et al.* (2018a)..

Our result (38) 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. Bultinck *et al.* (2019) that vortices in the IVC order parameter are unavoidable. The existence of vortices is neglected in the expansion (37) which assumes that is small everywhere. This implies that the expression (38) 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 Zhang *et al.* (2019)

(40) |

Here, is the Chern number for the valley and the parameter determines how quickly the form factor decays with which we take equal to to reproduce the Landau level form factors for . 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 Usov (1988). 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 for the valleys. For , we know that it is possible to write such a tight-binding model. For non-zero , 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 (40), the self-consistency equations can be solved by performing Fourier transform to real space. Following a series of straightforward steps, we get

(41) |

where is the modified Bessel function of the first kind and . The proportionality here indicates that we have dropped a constant positive factor given by which does not influence the competition between the two phases.

The expression (41) depends only on two dimensionless parameters: (i) the Chern number and (ii) 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 (41) 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 and is given in Fig. 7. For , 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 , 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 Bultinck *et al.* (2019)), we expect the transition from VP/SP to IVC to happen at even larger values of implying that VP/SP is the most energetically favorable insulator at half-filling whenever the bandwidth is relatively narrow.

#### iii.3.3 Effect of intervalley Hund’s coupling

We now consider the effect of intervalley Hund’s coupling term (22). This term is about two orders of magnitude smaller than the main part of the interaction so it can be neglected whenever two phases differ in the main part of the interaction energy. As a result, it is not expected to have any considerable effect on the competition between IVC and SP/VP orders or on the non-linear self-consistency equation for the IVC order.

However, it plays an important role in the competition between valley and spin polarized states which are degenerate on the level of . At half filling, the energy of the valley polarized state is unaffected by the Hund’s term whereas the energy of the spin polarized state is changed as

(42) |

which obviously reduces the energy and makes the spin-polarized state more energetically favorable.

In addition to settling the competition between spin and valley polarized phases, the Hund’s coupling affects the value of the insulating gap. Typically, the interaction-induced gaps are of the order of a few meV since the band width is of the same order as the interaction parameter (cf. Fig. 6). The inclusion of the Hund’s term for the spin polarized phase at half-filling reduces the energy of the lowest bands, thus increasing the gap. The Hund’s term has the opposite effect on the valley polarized phase since it reduces the energy of some of the empty bands, thereby, reducing the Hartree-Fock gap (although the ground state energy is unaffected) as shown in Fig. 10. Similarly, the Hund’s coupling term reduces the gap for the spin and valley polarized phase at quarter filling by lowering the energy of one of the excited states. The effect of the intervalley Hund’s coupling on the different gapped insulating states is shown schematically in Fig. 8.

We note here that the reduction of the correlated gap at quarter filling (valley and spin polarized) relative to that at half-filling (only spin polarized) may explain why the former is more difficult to observe experimentally compared to the latter. It is know that the Hartree-Fock mean field method overestimates the value of the gap. Thus, we expect the realistic gaps to be smaller than the ones computed here. In this case, Hund’s coupling could play an important role in stabilizing (destabilizing) the half (quarter) filled insulator by increasing (decreasing) the correlated gap by a few meVs.

## Iv Superconductivity

When the correlated insulator is doped away from half-filling, a superconducting phase is observed for temperature below about 3.5 K Liu *et al.* (2019). To understand how superconductivity develops in TDBLG, let us now consider what happens at a generic non-integer filling. Within a single valley, the Fermi surface is threefold symmetric but has neither time-reversal nor twofold symmetry. Instead, the Fermi surfaces of the two valleys are related by time-reversal. Hence, if superconducting pairing is to take place, we expect it between time-reverse partners residing in opposite valleys.

We can now speculate about the origin of the pairing interaction based on the knowledge that, both in theory and experiment, the parent insulating state is very likely a spin-polarized ferromagnet. Upon doping, we expect the magnitude of the spin-polarization to be reduced yielding a Ferromagnetic metal with spin-split Fermi surface. Similar to other ferromagnetic metals Saxena *et al.* (2000); Aoki *et al.* (2001); Huy *et al.* (2007), we expect Ferromagnetic spin fluctuations to be the primary source of pairing. This motivates the following simplification: rather than considering the projected Coulomb interaction (20), we instead consider the simplified Hamiltonian given by

(43) |

with the spin operator defined by

(44) |

The form of the Hamiltonian (43) can be justified in a more involved treatment by considering the RPA susceptibilities for different orders and identifying the ferromagnetic one as the leading instability. The ferromagnetic susceptibility is peaked at which implies a simple -independent pairing function within each valley ^{5}^{5}5Despite its independence on the momentum within each valley (-wave), the pairing function can change sign between valleys making it effectively -wave. This justifies the simplification in (43) where the interaction is taken to be -independent.

The interaction term can be rewritten as

(45) |

with . When performing the BCS decoupling, we restrict ourselves to pairing between time-reversed pairing which corresponds to and . In this case, we can define the gap function as . It satisfies the linearized BCS equation

(46) |

where is the Fermi velocity at point on the Fermi s surface . Choosing to be -independent, we can simplify (46)

(47) |

where is related to by some constant rescaling (coming from the Fermi surface integral), is the Pauli matrix vector in spin space and is a matrix in spin and valley spaces. Since pairing only takes place between opposite valleys, is proportional to or which corresponds to valley triplet or singlet respectively. Since the gap was chosen to be -wave within each valley, the overall antisymmetry of the gap function implies that former scenario corresponds to a spin-singlet whereas the latter corresponds to a spin-triplet . Here, is the vector which captures the direction of the spin state.

The symmetry of the superconducting order parameter is obtained by finding the pairing channel for which is positive and maximum. Substituting the spin-singlet and triplet gap functions in (47) yields

(48) | |||