Gate-Tunable Fractional Chern Insulators in Twisted Double Bilayer Graphene

We predict twisted double bilayer graphene to be a versatile platform for the realization of fractional Chern insulators readily targeted by tuning the gate potential and the twist angle. Remarkably, these topologically ordered states of matter, including spin singlet Halperin states and spin polarized states in Chern number $\mathcal C=1$ and $\mathcal{C}= 2$ bands, occur at high temperatures and without the need for an external magnetic field.

Introduction. Following the remarkable discovery of superconductivity and correlated states in magic angle twisted bilayer graphene [1][2][3][4][5][6], understanding the phase diagram resulting from electron-electron interactions in different Moiré heterostructures have attracted the interest of many.
As a natural progression to the study of the acclaimed twisted bilayer graphene, recent investigations have geared towards twisted double bilayer graphene (TDBG) systems [7][8][9][10][11][12], formed when two bilayers of graphene-instead of monolayers-are rotated with respect to each other. Experiments on TDBG have established it as a promising platform for interaction-driven states. In particular it has been reported [13,14] that spin polarized correlated insulator states arise in the narrow flat conduction band, and superconductivity emerges upon doping [15].
An important class of states that could arise due to the strong electron-electron interactions in an isolated flat band is fractional Chern insulators (FCI) [16,17]. These are lattice analogs of the fractional quantum Hall effect (FQH) in a Chern band with the advantage that they occur at zero magnetic field and with gaps that are potentially significantly larger than the gaps in the conventional fractional quantum Hall states. There has been a growing interest in investigating different kinds of FQH states in Moiré systems [18][19][20][21][22][23]. Most saliently, recent theoretical studies pointed at the existence of valley polarized FCI states at fractional filling of the topological C = 1 valence or conduction band of twisted bilayer graphene aligned with boron nitride [18][19][20]. While very encouraging, this setup generically exhibits strong particle-hole symmetry breaking terms and accompanying gapless competing phases making the necessary finetuning of experimental parameters challenging [18].
Motivated by these findings, here we explore the premises for FCI states in TDBG. TDBG servers as a natural platform to probe FCI states for a number of reasons. First, its layer configuration breaks the C 2 symmetry by default, hence removing the band touching at the Dirac points without needing to add a substrate layer. This results in separated bands around charge neutrality that could be individually studied. Secondly, while realistic models of TDBG that include effects of trigonal wrapping and particle-hole asymmetry terms generically lack the existence of magic angles [24], the TDBG bandwidth is nevertheless controlled by tuning the twist angle θ and more significantly via the application of an electric field U , providing an extra degree of tunability in experimental setups that are absent in twisted bilayer graphene systems in the regime of twist angles θ = 0.6 • − 1.4 • that we consider. This is to be contrasted with twisted bilayer graphene in the regime of very tiny twist angles θ 1 • where electric fields have a qualitatively different impact on the properties [25][26][27][28][29][30]. In TDBG the topology of the bands is also controlled by varying U and θ. In particular, the first conduction band, which we focus on, has Chern number values ranging between C = 0 and C = 3 [7][8][9] allowing us access not only to Landau levels alike bands with C = 1 but also to bands with C > 1 that are topologically distinct from continuum Landau levels. The appearance of an easily accessed C = 2 flat band opens the possibility of realizing novel C = 2 FCIs [31][32][33][34][35][36][37] in realistic materials, which was not reported yet in Moiré systems without an external magnetic field.
In this Letter, we provide compelling evidence for the existence of FCI states in TDBG through a detailed mi- croscopic study of projected Coloumb interactions onto the first conduction band for both C = 1 and C = 2 regimes. We find a variety of ferromagnetic and spin singlet FCI states at different filling factors that are highly tunable i.e, could be accessed by moving around in the U − θ space. Setup. We consider electrons interacting via the screened Coulomb potential in the TDBG Moiré superlattice. We choose the Yukawa potential V (q) = e 2 4π r 0S 2π √ |q| 2 +κ 2 to describe the screening of the Coulomb interaction, where e is the electron charge, 0 is the dielectric constant of vacuum, r ≈ 4 is the relative dielectric constant of the material [9], S is the area of the Moiré superlattice, and κ measures the screening strength which we set as κ = 1/a M with a M the lattice constant of TDBG [38]. Unless otherwise stated, we assume both polarized spin and valley degrees of freedom for electrons. When the electrons are doped above the charge neutrality point, the first conduction band is partially filled. If the first conduction band is isolated from other bands below and above, it is fair to project the total Hamiltonian to this active band, leading to  (1) can be labeled by a two-dimensional momentum (K 1 , K 2 ). Band isolation. One of the most attractive advantages of TDBG is the high tunability of low-energy bands via the twist angle θ and the vertical voltage bias U . In particular, the first conduction band can be easily isolated and shows a rich phase diagram of band Chern number C when the twist angle θ and vertical voltage U are varied [8]. Before studying the interaction driven physics, we first explore the U −θ space to find the regions where the first conduction band is isolated such that the projected Hamiltonian (1) applies. In order to obtain a realistic model of TDBG, we choose w 1 = 100meV and w 0 /w 1 = 0.7 throughout this work, where w 0 and w 1 are the AA and AB tunneling strengths between two sheets of bilayer graphene, respectively [38]. With this choice of parameters, we calculate the indirect gap ∆ ≡ min(∆ e , ∆ v ) of the first conduction band, where ∆ e and ∆ v are the indirect gaps to the next excited band and the valence band, respectively. We observe three main regions with significant ∆ in the U − θ space. The isolated first conduction band has C = 0, 1 and 2 in these three regions, respectively ( Fig. 1) (similar results were reported in Ref. [8] for different model parameters). In the following, we will focus on the C = 1 and C = 2 regions in the phase diagram, and probe the existence of robust FCIs in both regions.
FCIs in the C = 1 region. Now we address the possibility of FCIs stabilized by the screened Coulomb interaction in the C = 1 region. We first examine ν = 1/3, where one may in general expect robust FCIs as the lattice analogs of the celebrated Laughlin state [39]. Indeed, we observe clear three-fold ground-state degeneracies in the energy spectrum of the Hamiltonian (1) for U = 20meV − 30meV and θ = 0.6 • − 0.8 • [ Fig. 2(a)], where the three approximately degenerate ground states have momenta (K 1 , K 2 ) which are consistent with the prediction of Haldane statistics for the ν = 1/3 FCIs in a C = 1 band [40][41][42]. The three-fold ground-state degeneracy persists during the magnetic flux insertion through the handles of the toroidal system, thus further confirming the robustness of finite size fingerprint of a topological degeneracy [ Fig. 2(b)]. The energy gap separating the three ground states from excited states becomes much larger than the ground-state splitting as the system size grows, and is very likely to survive in the thermodynamic limit as suggested by a finite-size scaling [ Fig. 2(c)]. Remarkably, the gap corresponds to a temperature of about 10 Kelvin, which is at least an order of magnitude higher than required by the conventional fractional quantum Hall states in two-dimensional electron gases. Importantly this energy scale is still significantly smaller than the band gap, ∆, confirming the validity of the band projection. The entanglement spectroscopy of the ground-state manifold with the particle-cut entanglement spectrum (PES) [40,43,44] further corroborate that the ground states are topologically nontrivial. By dividing the whole system into N A and N − N A electrons and labeling each PES level by the total momentum (K A 1 , K A 2 ) of those N A electrons, we find a clear entanglement gap ∆ξ ≈ 0.3 separating the low-lying PES levels from higher ones [ Fig. 2(d)]. The number of levels below the entanglement gap exactly matches the pertinent counting of quasihole excitations in the ν = 1/3 Abelian FCIs [40][41][42], which rules out competing possibilities such as charge density waves.
We have also considered other filling factors. At ν = 2/3, we find the particle-hole conjugate of the ν = 1/3 FCIs reported above [38]. Moreover, we observe vestiges of a possible five-fold ground-state degeneracies at ν = 2/5 and ν = 3/5 [38], which may suggest the ν = 2/5 Jain state and its particle-hole conjugate. However, these states are fragile against band dispersion and compete with spinful states (as discussed below).
FCIs in the C = 2 region. Previous studies have reported novel FCIs residing in flat bands with higher Chern numbers [31][32][33][34][35][36][37]. Unlike FCIs in |C| = 1 flat bands, these high-C FCIs do not have usual continuum FQH states as counterparts. In particular, application of composite fermion theory to Bloch bands predicts a series of high-C FCIs at ν = r/(rk|C| + 1) [36], where k > 0 is the number of flux attached to each particle in the composite fermion theory (k should be even for fermions) and r is the number of fully filled compositefermion bands. Note that r can be negative which corresponds to negative flux attachment. The energy gap usually decays with the decreasing of ν-specifically with increasing denominators-which makes the states more fragile and hard to observe both in finite size calculations in experiments. In the following we consequently consider relatively large fillings in the ν = r/(2rk + 1) branch for the C = 2 region of TDBG.
We first examine ν = 1/3 by setting k = 2 and r = −1. Remarkably, we find clear three-fold ground-state degeneracies for relatively large systems (N ≥ 8 electrons) with U = 50meV − 70meV and θ = 1.2 • − 1.4 • [ Fig. 3(a)]. While the energies vary more during the flux insertion than at ν = 1/3 in the C = 1 region, the three-fold topological degeneracy still persists in the spectral flow [ Fig. 3(b)]. The finite-size scaling for the three available system sizes within our computational limit suggests that an energy gap of ∼ 5 Kelvin probably survives in the thermodynamic limit [ Fig. 3(c)]. Moreover, we observe an entanglement gap ∆ξ ≈ 0.2 in the PES [ Fig. 3(d)]. Albeit being small, this gap gives a low-energy PES manifold in which the number of levels is the same as that for the ν = 1/3 Laughlin FCIs in |C| = 1 bands. This coincidence of the PES counting between FCIs at the same filling in |C| = 2 bands and |C| = 1 bands have been noticed previously for other |C| = 2 FCIs [35,36], thus further confirming the topological nontrivial property of the ν = 1/3 states in the C = 2 region of TDBG.
Spinful FCIs. So far we have assumed both valley and spin polarization. Now motivated by experiments we keep valley polarization [14], but bring the spin degree of freedom back. In this case, as both the band dispersion and the interaction are independent on the spin flavor if one neglects Hund's coupling terms that are generically much weaker than the Coulomb interactions, the spinful many-body Hamiltonian [38] is SU(2) symmetric within each valley, which allows us to label each energy level by the total spin S and its z-component S z = (N ↑ − N ↓ )/2, with N ↑ and N ↑ the number of spin-up and spin-down electrons, respectively. At ν = 1/3 in both the C = 1 and C = 2 regions, we find three-fold ground-state degeneracies in all S z sectors and the ground energies with different S z are identical (Fig. 4). This confirms the assertion that the ν = 1/3 FCIs observed in both C = 1 and C = 2 regions are indeed ferromagnetic with total spin S = N/2.
On the other hand, we find that such ferromagnetism can disappear at other filling factors, leading to spinsinglet ground states. For example, at ν = 2/5 in the C = 1 region with U around 20meV and θ = 0.65 • −0.75 • , the ground states are in the S z = 0 sector, thus having total S = 0. Remarkably, five-fold ground-state degeneracies appear in this case (Fig. 5), strongly suggesting the Halperin (332) state as the ground state. As it was found in other graphene based Moiré systems that ferromagnetic states are favored with the increasing of w 0 /w 1 [19], there could be a phase transition from the Halperin (332) state to the ν = 2/5 Jain state at w 0 /w 1 larger than the value 0.7 chosen in this work.
Discussion. In this work we have established twisted double bilayer graphene as a flexible platform for a plethora of fractional Chern insulators which are stabilized by tuning a gate voltage and the twist angle. Remarkably, these exotic states occur at high temperatures and in parameter regimes that are directly experimentally accessible. This may likely enable the experimental discovery of states with no continuum analogue such as the C = 2 fractional Chern insulator at filling ν = 1/3, and provides an intriguing path toward the possibility of technological applications. In the latter context, the investigation of possible non-Abelian states provides a particularly intriguing outlook.

SUPPLEMENTARY MATERIAL
In this supplementary material, we provide details about the single-particle model of twisted double bilayer graphene (TDBG), and more numerical results for FCIs in the C = 1 region.

TDBG model
We consider two sheets of AB stacked bilayer graphene (BLG) which are twisted with respect to one another, with the ABAB stacking pattern. Each AB stacked BLG sheet is modeled by the following single-particle Hamiltonian near the valley K + = 4π 3a (1, 0) of its Brillouin zone, where a ≈ 2.46Å is the lattice constant of graphene, and the basis is (ψ A1 (k), ψ B1 (k), ψ A2 (k), ψ B2 (k)) T with A, B the sublattice indices of monolayer graphene and 1, 2 the layer indices in the AB stacked BLG sheet. t 0 , t 1 , t 3 and t 4 are the A 1 − B 1 , A 1 − B 2 , B 1 − A 2 and B 1 − B 2 hopping strengths in the AB stacked BLG sheet, respectively, δ is the onsite energy difference between A and B sites, and U i 's are the gating voltage across the system. Throughout the work we adopt (t 0 , t 1 , t 3 , t 4 , δ) = (2610, 361, 283, 138, 15)meV [8,45].
When the two AB stacked BLG sheets are twisted, we focus on the Moiré Brillouin zone of TDBG formed near K + . The two primitive reciprocal lattice vectors of TDBG are chosen as , with θ the twist angle and a M = a/(2 sin θ 2 ) the lattice constant of TDBG. Let us denote K t + = R θ/2 K + and K b + = R −θ/2 K + , where R θ a counter-clockwise rotation around the z-axis in the momentum space. The single-particle Hamiltonian of TDBG for each spin favor can then be written as where ψ l=t,b (k) = (ψ A 1,l (k), ψ B 1,l (k), ψ A 2,l (k), ψ B 2,l (k)) T is the basis for electrons in top and bottom BLG sheets, respectively, h θ (k) = h(R θ k) with h(k) given in Eq. (S1), and q 0 = R −θ/2 K + − R θ/2 K + , q 1 = R 2π/3 q 0 and q 2 = R −2π/3 q 0 . We choose as U 1 = U/2, U 2 = U/6 for the top BLG sheet and U 1 = −U/6, U 2 = −U/2 for the bottom BLG sheet. The Moiré hoppings T j 's in TDBG are given by where w 0 and w 1 are the inter-sheet hopping strengths between AA and AB sites, respectively. In our numerics, we take w 1 = 100 meV and w 0 = 0.7w 1 .
For each k 0 in the MBZ, by writing k in the single-particle Hamiltonian H [Eq. (S2)] as k 0 + mG 1 + nG 2 and setting integers m, n = −d, ..., d, we can construct H(k 0 ) as a matrix of dimension 8(2d + 1) 2 . The eigenvalues and eigenvectors of this Hamiltonian matrix then give us the band energies and eigenvectors of TDBG at k 0 .
Here V (q) is the Fourier transform of the interaction potential. For the Yukawa potential, V (q) = e 2 4π S 2π √ |q| 2 +κ 2 , where e is the electron charge, is the dielectric constant of the material, S is the area of the Moiré superlattice, i N2 G 2 with k 1 i = 0, 1, · · · , N 1 − 1 and k 2 i = 0, 1, · · · , N 2 − 1. As the matrix elements of the single-particle Hamiltonian Eq. (S2) are identical for spin-up and spin-down electrons, we can generalize the projected many-body total Hamiltonian [Eq. (1) in the main text] to the spinful case: where c † k,σ (c k,σ ) creates (annihilates) an electron with momentum k and spin σ in the first conduction band (per valley). Note that both the band dispersion E(k) and the interaction matrix element V k1k2k3k4 are independent on the spin, leading to an SU(2) symmetric Hamiltonian. ν = 2/3, ν = 2/5 and ν = 3/5 in the C = 1 region Here we assume both spin and valley polarization. We first examine ν = 2/3 in the C = 1 region. With the same parameters as where we find the ν = 1/3 FCIs, we again observe nice three-fold ground-state degeneracies at ν = 2/3, as shown in Fig. S1(a). We identify these states as the particle-hole conjugates of the ν = 1/3 FCIs in the C = 1 region shown in the main text.
We then report the numerical results at ν = 2/5 and ν = 3/5 in the C = 1 region of TDBG. In both cases, our numerical results do not suggest well developed FCIs for finite systems within our computational limit and the model parameters which we choose. However, we still observe five lowest energy levels in the correct momentum sectors predicted for the ν = 2/5 Jain FCIs and their ν = 3/5 particle-hole conjugates, although the splitting between these states are significantly larger than the ν = 1/3 cases shown in the main text [Figs. S1(b) and S1(c)]. Therefore, the ν = 2/5 Jain FCIs and their ν = 3/5 particle-hole conjugates could be stabilized for larger system sizes or modified model parameters (especially w 0 and w 1 ).