Flatbands in twisted double bilayer graphene

TL;DR

This paper investigates the electronic structure of twisted double bilayer graphene, showing how flatbands with narrow bandwidths can be achieved and tuned by parameters like twist angle, pressure, and electric fields, with implications for novel electronic states.

Contribution

It provides a detailed analysis of flatband emergence in twisted bi-bilayer graphene, highlighting how system parameters influence bandwidth and gap properties, and compares it to twisted bilayer graphene.

Findings

Flatbands are generally flatter in twisted bi-bilayer graphene than in twisted bilayer graphene.

Vertical pressure can increase the first magic angle from 1.05° to about 1.5° at 2.5 GPa.

Moderate electric fields can lift degeneracy and enhance gaps near the Dirac point.

Abstract

Flatbands with extremely narrow bandwidths on the order of a few mili-electron volts can appear in twisted multilayer graphene systems for appropriate system parameters. Here we investigate the electronic structure of a twisted bi-bilayer graphene, or twisted double bilayer graphene, to find the parameter space where isolated flatbands can emerge as a function of twist angle, vertical pressure, and interlayer potential differences. We find that in twisted bi-bilayer graphene the bandwidth is generally flatter than in twisted bilayer graphene by roughly up to a factor of two in the same parameter space of twist angle $\theta$ and interlayer coupling $\omega$, making it in principle simpler to tailor narrow bandwidth flatbands. Application of vertical pressure can enhance the first magic angle in minimal models at $\theta \sim 1.05^{\circ}$ to larger values of up to $\theta \sim 1.5^{\circ}$ when $ P \sim 2.5$~GPa, where $\theta \propto \omega/ \upsilon_{F}$. Narrow bandwidths are expected in bi-bilayers for a continuous range of small twist angles, i.e. without magic angles, when intrinsic bilayer gaps open by electric fields, or due to remote hopping terms. We find that moderate vertical electric fields can contribute in lifting the degeneracy of the low energy flatbands by enhancing the primary gap near the Dirac point and the secondary gap with the higher energy bands. Distinct valley Chern bands are expected near $0^{\circ}$ or $180^{\circ}$ alignments.

Tables

Table 1. Table 1: External pressure P ​ ( ω s ) 𝑃 subscript 𝜔 𝑠 P(\omega_{s}) and associated interlayer coupling ω s ​ ( P ) subscript 𝜔 𝑠 𝑃 \omega_{s}(P) for commensurate stacking geometries ( s 𝑠 s = AA, AB, BA) illustrated in Fig.  8 . We list the fitting coefficients for P ​ ( ω ) 𝑃 𝜔 P(\omega) in Eq. ( 10 ) A 𝐴 A (GPa/eV 2 ), for B 𝐵 B (GPa/eV), and C 𝐶 C (GPa), and the inverse fit for ω ​ ( P ) 𝜔 𝑃 \omega(P) in Eq. ( 11 ) A ′ superscript 𝐴 ′ A^{\prime} (eV), for B ′ superscript 𝐵 ′ B^{\prime} (eV 2 ), and C ′ superscript 𝐶 ′ C^{\prime} (eV 2 /GPa).

		LDA			EXX+RPA
Stacking (s)		$A$	$B$	$C$	$A$	$B$	$C$
AA		473.6	-22.17	-0.9011	543.3	-61.01	0.735
AB		306.5	-36.56	-0.04339	324.7	-35.47	-0.4671
Stacking (s)		$A^{\prime }$	$B^{\prime }$	$C^{\prime }$	$A^{\prime }$	$B^{\prime }$	$C^{\prime }$
AA		0.0234	0.0025	0.0021	0.0561	0.0018	0.0018
AB		0.0596	0.0037	0.0033	0.0546	0.0044	0.0031

Table 2. Table 2: The Chern numbers table of tBBG for twist angles near the 0 ∘ (Type-I) and 180 ∘ (Type-II) alignment with and without remote hoping parameters at the angles θ = 1.5 ∘ , 1.05 ∘ 𝜃 superscript 1.5 superscript 1.05 \theta=1.5^{\circ},1.05^{\circ} and 0.7 ∘ whose band structures are represented in Fig.  10 . The Chern numbers presented here are for low energy conduction (C) and valence bands (V), and one higher energy bands in each conduction (C + + 1) and valence (V − - 1) bands. The interlayer potentials are respectively Δ ′ = 0.03 superscript Δ ′ 0.03 \Delta^{\prime}=0.03  eV for θ ∼ 0 ∘ similar-to 𝜃 superscript 0 \theta\sim 0^{\circ} and Δ ′ = 0.05 superscript Δ ′ 0.05 \Delta^{\prime}=0.05  eV for θ ∼ 180 ∘ similar-to 𝜃 superscript 180 \theta\sim 180^{\circ} that are large enough to isolate the low energy bands.

		Minimal				Remote
Bands		V ($-$1)	V	C	C ($+$1)	V ($-$1)	V	C	C ($+$1)
		$\theta ={1.5}^{\circ }$
Type-I		1	$-$2	2	$-$1	2	$-$2	2	$-$1
Type-II		0	1	1	0	0	1	1	0
		$\theta ={1.05}^{\circ }$
Type-I		0	$-$2	2	0	0	$-$2	2	0
Type-II		$-$1	1	1	$-$1	$-$1	1	1	$-$1
		$\theta ={0.7}^{\circ }$
Type-I		$-$2	0	0	2	$-$2	0	0	$\dots$
Type-II		0	0	0	0	0	0	0	0

Equations

$$H_{\mathrm{tBBG}}(\theta )=\left(\begin{array}{cccc}{h}_t^{+}+{\bar{\Delta }}_1& t_s^{+}& 0& 0\\ t_s^{+†}& h_b^{+}+{\bar{\Delta }}_2& T(\mathbf{r})& 0\\ 0& T^{†}(\mathbf{r})& h_t^{-}+{\bar{\Delta }}_3& t_s^{-}\\ 0& 0& t_s^{-†}& h_b^{-}+{\bar{\Delta }}_4\end{array}\right),$$

$$t_{\mathrm{AB}}^{\pm }=\left(\begin{array}{cc}-{\upsilon }_4{\pi }^{\pm †}& -{\upsilon }_3{\pi }^{\pm }\\ t_1& -{\upsilon }_4{\pi }^{\pm †}\end{array}\right),t_{\mathrm{BA}}^{\pm }=t_{\mathrm{AB}}^{\pm †}$$

$$T(\mathbf{r})=\sum _{j=0,\pm }{e}^{-i{\mathbf{Q}}_j\mathbf{r}}{T}_{l,l^{\prime }}^j,$$

$$T^0=\left(\begin{array}{cc}{\omega }^{\prime }& \omega \\ \omega & {\omega }^{\prime }\end{array}\right),T^{\pm }=\left(\begin{array}{cc}{\omega }^{\prime }& \omega e^{\mp i2\pi /3}\\ \omega e^{\pm i2\pi /3}& {\omega }^{\prime }\end{array}\right)$$

$${\omega }_{l{l}^{\prime }}={\int }_{\mathrm{cell}}\mathrm{d}\vec{d}{H}_{l{l}^{\prime }}(K,\vec{d})\simeq \sum _s\frac{H_{l{l}^{\prime }}(K,{\vec{d}}_s)}{3}$$

$$\omega$$

$${\omega }^{\prime }$$

$$U_{\mathrm{eff}}=\frac{e^2}{4\pi {\epsilon }_r{\epsilon }_0{l}_M}\exp (-l_M/{\lambda }_{\mathrm{D}})$$

$${\theta }_n^{\circ }$$

$$P$$

$$\omega$$

$$C_{\upsilon }={\int }_{\mathrm{mBZ}}{\mathrm{d}}^2\vec{k}{\Omega }_n(\vec{k})/(2\pi )$$

$${\Omega }_n(\vec{k})=-2\sum _{n^{\prime }\ne n}\mathrm{Im}\left[\frac{⟨u_n|\frac{\partial H}{\partial k_x}|u_{n^{\prime }}⟩⟨u_{n^{\prime }}|\frac{\partial H}{\partial k_y}|u_n⟩}{{\left(E_{n^{\prime }}-E_n\right)}^2}\right]$$

$$U=\frac{e^2}{4\pi {\epsilon }_r{\epsilon }_0{l}_M}$$

$$H_{eff}(k)=\left(\begin{array}{cccc}{\Delta }_1& h^{+}(k)& 0& 0\\ {h^{+}k}^{†}& -{\Delta }_2& h_t& 0\\ 0& h_t^{†}& {\Delta }_2& h^{-}(k)\\ 0& 0& {h^{-}(k)}^{†}& -{\Delta }_1\end{array}\right),$$

$$h^{\pm }(k)=\frac{{\upsilon }_{\mathrm{F}}^2}{3\omega }{\left[(k_x-i{k}_y)e^{i\theta /2}\right]}^2$$

$$H_t=\omega \left[\begin{array}{cc}0& 0\\ 1& 0\end{array}\right],\omega \left[\begin{array}{cc}0& 0\\ e^{-i2\pi /3}& 0\end{array}\right],\omega \left[\begin{array}{cc}0& 0\\ e^{i2\pi /3}& 0\end{array}\right]$$

$$H=\sum _{a=t,b}{\Psi }_a^{†}(\mathbb{k})H^a(\mathbb{k},\theta ){\Psi }_a(\mathbb{k})+\sum _{j=0,\pm }\left({\Psi }_t^{†}(\mathbb{k})T_j{\Psi }_b(\mathbb{k}+{\mathbb{Q}}_j)+h.c\right)$$

$$H^t(k,\theta )=\left(\begin{array}{cccc}{\Delta }_1& {\upsilon }_0(\xi k_x-i{k}_y)e^{-i\theta /2}& -{\upsilon }_4(\xi k_x-i{k}_y)e^{-i\theta /2}& {\upsilon }_3(\xi k_x+i{k}_y)e^{-i\theta /2}\\ {\upsilon }_0(\xi k_x+i{k}_y)e^{i\theta /2}& {\Delta }_1& t_{\perp }& -{\upsilon }_4(\xi k_x-i{k}_y)e^{-i\theta /2}\\ -{\upsilon }_4(\xi k_x+i{k}_y)e^{i\theta /2}& t_{\perp }& -{\Delta }_2& {\upsilon }_0(\xi k_x-i{k}_y)e^{-i\theta /2}\\ {\upsilon }_3(\xi k_x-i{k}_y)e^{i\theta /2}& -{\upsilon }_4(\xi k_x+i{k}_y)e^{i\theta /2}& {\upsilon }_0(\xi k_x+i{k}_y)e^{i\theta /2}& -{\Delta }_2\end{array}\right)$$

$$H^b(k,\theta )=\left(\begin{array}{cccc}{\Delta }_2& {\upsilon }_0(\xi k_x-i{k}_y)e^{i\theta /2}& -{\upsilon }_4(\xi k_x-i{k}_y)e^{i\theta /2}& {\upsilon }_3(\xi k_x+i{k}_y)e^{i\theta /2}\\ {\upsilon }_0(\xi k_x+i{k}_y)e^{-i\theta /2}& {\Delta }_2& t_{\perp }& -{\upsilon }_4(\xi k_x-i{k}_y)e^{i\theta /2}\\ -{\upsilon }_4(\xi k_x+i{k}_y)e^{-i\theta /2}& t_{\perp }& -{\Delta }_1& {\upsilon }_0(\xi k_x-i{k}_y)e^{i\theta /2}\\ {\upsilon }_3(\xi k_x-i{k}_y)e^{-i\theta /2}& -{\upsilon }_4(\xi k_x+i{k}_y)e^{-i\theta /2}& {\upsilon }_0(\xi k_x+i{k}_y)e^{-i\theta /2}& -{\Delta }_1\end{array}\right)$$

$$T_j\simeq \left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 3\omega & 0& 0& 0\end{array}\right)$$

$$H_{eff}=\left(\begin{array}{cc}{\Delta }_1& \frac{-1}{(3\omega {)}^3}\left[(3\omega {)}^2{\upsilon }_3^2{q}^2+{\upsilon }_F^4{q}^4-6\omega {\upsilon }_3{\upsilon }_F^4{q}^4\right]\\ \frac{-1}{(3\omega {)}^3}\left[(3\omega {)}^2{\upsilon }_3^2{q^{†}}^2+{\upsilon }_F^4{q^{†}}^4-6\omega {\upsilon }_3{\upsilon }_F^4{q^{†}}^4\right]& -{\Delta }_1\end{array}\right)$$

$${\lambda }^{\pm }=\pm \sqrt{{\Delta }_1^2+{\left|\frac{1}{(3\omega {)}^3}\left[(3\omega {)}^2{\upsilon }_3^2{q}_m^2+{\upsilon }_F^4{q}_m^4-6\omega {\upsilon }_3{\upsilon }_F^2{q}_m^4\right]\right|}^2},$$

$$W=({\lambda }^{max}-{\lambda }^{min})\simeq \frac{1}{(3\omega {)}^3}\left[(3\omega {)}^2{\upsilon }_3^2\frac{|\Delta K{|}^2}{4}+{\upsilon }_F^4\frac{|\Delta K{|}^4}{16}-6\omega {\upsilon }_3{\upsilon }_F^2\frac{|\Delta K{|}^4}{16}\right]-{\Delta }_1$$

$$(6\omega {\upsilon }_3{\upsilon }_F^2-{\upsilon }_F^4)|\Delta K{|}^4-4(3\omega {)}^2{\upsilon }_3^2|\Delta K{|}^2+16(3\omega {)}^3{\Delta }_1=0$$

Full text

I Introduction

Research on twisted hybrid van der Waals 2D materials has seen recently a new surge of interest following experimental observations of exotic quantum phases due to strong electron correlations Cao ; kyounghwan and especially signatures of unconventional superconductivity pjarillo_superconductivity ; Yankowitz ; planckian in twisted bilayer graphene (tBG), raising hopes of finding new clues for understanding analogous behaviors seen in more complex systems. complexoxides In tBG the spatial variation of interlayer coupling modifies the intrinsic Dirac cone band structure of graphene in such a way that the band dispersion is almost completely suppressed at the so called magic twist angles. bistritzer2011 When the bandwidth $W$ of these low energy bands is sufficiently narrow it is possible to achieve the $U/W\gtrsim 1$ condition that makes the effective Coulomb repulsion $U$ more dominant. A considerable body of literature has formed recently on the Coulomb interaction driven broken symmetry phases Ochi ; Pizarro ; Venderbos ; Phillips2018 ; koshino_wannier ; Kang1 ; Kang2 and superconductivity in tBG flatbands Fidrysiak ; Roy ; Huang ; Tanmoy ; Liu ; Peltonen ; Kennes ; Isobe ; Vishwanath ; Vishwanath2 ; Wu ; Guinea ; Gonzalez ; Su ; Lian ; Laksono ; Tang ; Dante1 in an effort to elucidate the nature of the superconducting phases. Analogous observations of Coulomb interaction driven correlated phases and superconductivity have been observed in ABC trilayer graphene (TG) on hexagonal boron nitride (hBN) guorui_mott2018 ; guorui_superconductivity2019 ; guorui2019b where the flattening of the low energy bands is facilitated by the presence of a vertical electric field that introduces a band gap at the primary Dirac point of a chiral two-dimensional electron gas (2DEG), prljung ; prlzhang or in twisted gapped Dirac materials. Dante2 ; javvaji ; mos2bilayers It was suggested that the low energy flatbands footnote could have well defined valley Chern numbers and give rise to spontaneous quantum Hall phases when the band degeneracy is lifted by Coulomb interactions. todadri ; prljung

The proposals of flatbands in several types of multilayer graphene materials is suggesting that they can arise in a large variety of 2D material combinations provided that we choose the appropriate intrinsic electronic structure of each layer and their interlayer coupling. todadri In this work we study the flatband bandwidth phase diagram of twisted BG/BG system, that we refer to as twisted bi-bilayer graphene (tBBG) or twisted double bilayer graphene, which consists of two bilayer graphene units with a twist. We assess for this system the effect of the interlayer coupling strength and the interlayer potential differences between the layers in the resulting bandwidth of the low energy flatbands. It is expected that the smaller parabolic band dispersion slopes at low energy in a BG can favor the formation of flatbands upon interlayer hybridization. This manuscript is structured as follows. In Sec. II we introduce the theoretical details of the continuum model Hamiltonian used to formulate the problem. In Sec. III we present the phase diagram of $U/W$, the flatband bandwidth and gaps as a function of different system parameters, such as the twist angle, the interlayer coupling strength, and the interlayer potential differences due to a vertical electric field. In Sec. IV we discuss the valley Chern number phase diagrams, and then close the paper in Sec. V with the summary and conclusions.

II Model Hamiltonian for twisted bi-bilayer graphene

Models proposed in the literature to capture the electronic structure of tBG relied either on tight-binding calculations lopes ; shallcross ; koshino often based on the distance dependent two center approximation models for the hopping terms between interlayer carbon atoms, or by using other more sophisticated parametrizations. jung2014 ; kaxirasparam The successful formulation of a rigorous moire bands theory on the basis of the moire pattern superlattice bistritzer2011 allows to obtain accurate continuum models for the Hamiltonian informed from first principles calculations. jung2014 In the present work we extend the continuum model of Bistritzer-MacDonald model for the tBG bistritzer2011 to the case of tBBG. The Hamiltonian of tBBG at valley $K$ that we use captures the interlayer coupling between the twisted layers through a first harmonic stacking-dependent interlayer tunneling function. Schematic representations of commensurate and twisted BG/BG structures are shown in Fig. 1. We write the Hamiltonian of twisted top ($+$) and bottom ($-$) bilayer graphene Hamiltonian subject to $\Delta_{i}$ intralayer potentials as

[TABLE]

where $h^{\pm}_{t/b}=h_{t/b}(\pm\theta/2)$ such that the relative twist angle between the bilayers is $\theta$. The top and bottom BG are labeled through the positive/negative ($+/-$) rotation signs, while in turn we have top/bottom ($t/b$) graphene layers within each BG that are coupled through the matrices $t_{s}^{\pm}$ that we define later. The site potentials for each graphene layer $\Delta_{i}$ are mapped on its sublattices through $\bar{\Delta}_{i}=\Delta_{i}\mathbb{1}$ where $i=1,\,2,\,3,\,4$ are the layer labels from top to bottom, and $\mathbb{1}$ is a $2\times 2$ identity matrix. Potential differences can give rise to band gaps at the primary Dirac point of each BG and shift the associated band edges. We will discuss later on the effects of these intralayer potentials in the electronic structure of the flatbands. The Hamiltonian of graphene is given by $h^{\pm}_{l}(\theta)=h^{\pm}(\theta)+\delta\left(\mathbb{1}-ls\sigma_{z}\right)/2$ where the second term adds a $\delta=0.015$ eV sublattice potential at the higher energy dimer sites at the $t/b$ layers $l=\pm$, accuratebilayer that depends on AB or BA stacking $s=\pm$ respectively. The Dirac Hamiltonian given by $h(\theta)=\upsilon_{F}\hat{R}_{-\theta}{\bf p}\cdot\sigma_{xy}$ includes a phase shift due to a rotation $\hat{R}_{-\theta}$ such that $e^{\pm i\theta_{\bf p}}\rightarrow e^{\pm i(\theta_{\bf p}-\theta)}$, where $\sigma_{xy}=(\sigma_{x},\sigma_{y})$ and $\sigma_{z}$ are the graphene sublattice pseudospin Pauli matrices, and the momentum is defined in the $xy$-plane ${\bf p}=(p_{x},p_{y})$, where we assume $K$ valley unless stated otherwise. The Fermi velocity $\upsilon_{\rm F}=\upsilon_{0}$ defined from $\upsilon_{i}=\sqrt{3}\left|t_{i}\right|a/2\hbar$ is related to the intralayer nearest neighbor hopping term $t_{0}=-2.6$ eV within the local density approximation (LDA), ldahopping while an enhanced $t_{0}=-3.1$ eV and ab initio interlayer tunning captures better the experimental moire band features. dillon2015 The interlayer coupling model of a bilayer graphene is given by

[TABLE]

satisfying $t_{s=+}=t_{s=-}^{{\dagger}}$ for AB or BA ($s=\pm 1$) stacking dependent interlayer coupling that consists of a minimal coupling term $t_{s}=t_{1}(\sigma_{x}-is\sigma_{y})/2$ plus remote hopping contributions through $t_{3}=0.283$, $t_{4}=0.138$ terms, giving rise to trigonal warping and electron-hole asymmetry. The $\pi^{\pm}$ operators include the phases due to $\pm\theta/2$ layer rotation. The type II AB/BA bi-bilayers near 180*∘* alignment can be modeled by controlling the stacking of bottom BG from AB to BA by using $s=-1$ for the bottom BG. The interlayer tunneling is defined as the Hamiltonian matrix element at the Dirac point $t_{1}=H_{BA^{\prime}}(K,\vec{d}_{\rm AB})$ between $B$ to $A^{\prime}$ sites from bottom to top layer for AB stacking when both atomic sites are vertically aligned and assume $t_{1}=0.361$ eV at zero pressure within LDA. accuratebilayer We can identify the interlayer tunneling with the first harmonic expansion coefficient of the interlayer coupling such that $t_{1}=3\omega$, jung2014 and for simplicity we use the same AB stacking tunneling within each Bernal BG and the twisted interfaces. The minimal model approximation uses $\delta=t_{3}=t_{4}=0$ in Eq. 15. The presence of remote hopping terms will lead to broadening of the low energy flatbands and enhancemenet in electron-hole asymmetry. This behavior is not strange since the $t_{3}$ trigonal warping widens the range of band touching points at three points away from the Dirac point at directions connecting the $K$ points with $\Gamma$, accuratebilayer and the $t_{4}$ term breaks the intrinsic electron hole symmetry of bilayer graphene. mccann

The moire Brillouin zone (mBZ) orientation is preserved when the top and bottom graphene layers rotate symmetrically in opposite senses. In the small angle approximation the interlayer coupling Hamiltonian is given by

[TABLE]

where the three ${\bf Q}_{j}$ vectors ${\bf Q}_{0}=K_{D}\theta(0,-1)$ and ${\bf Q}_{\pm}=K\theta(\pm\sqrt{3}/2,1/2)$ are proportional to twist angle $\theta$ and $K_{D}=4\pi/3a$ is the Brillouin zone corner length of graphene, whose lattice constant is $a=2.46~{}\AA$, and here the indices $l,\,l^{\prime}$ label the sublattices of neighboring twisted surface layers. The interlayer coupling matrices between the two rotated adjacent layers are given by

[TABLE]

using a form that distinguishes interlayer tunneling matrix elements $\omega=\omega_{BA^{\prime}}$ and $\omega^{\prime}=\omega_{AA^{\prime}}$ for different and same sublattice sites between the layers. The convention taken here for the $T^{j}$ matrices jung2014 assume an initial AA stacking configuration $\tau=(0,0)$ and differs by a phase factor with respect to the initial AB stacking $\tau=(0,a/\sqrt{3})$. bistritzer2011 The greater interlayer separation $c$ compared to the carbon-carbon distances $a_{\rm CC}$ lead to slowly varying interlayer tunneling function $T({\bf r})$ and the moire patterns can often be accurately described within a first harmonic approximation bistritzer2011 ; jung2014 . In this limit, and assuming no corrugation effects, the interlayer coupling strength can be well approximated by a single parameter $\omega=\omega^{\prime}$ whose value was calculated within the local density approximation (LDA) to be $\omega\sim 0.113$ eV when averaged for every stacking at a fixed interlayer distance $c_{\rm AB}=3.35~{}\AA$ of AB stacking. A somewhat weaker $\omega\sim 0.098$ eV is expected when the interlayer relaxations for farther AA interlayer distance $c_{\rm AA}=3.57~{}\AA$ within LDA is accounted for in the averaging process. jung2014 The interlayer tunneling matrix elements $H_{ll^{\prime}}(K,\vec{d})$ are evaluated at the Dirac point $K$ for a commensurate system with stacking sliding vector $\vec{d}$ through the lattice Fourier transform of the distant real space hopping terms connecting the sites $l$ and $l^{\prime}$. The tunneling matrix elements $\omega$, $\omega^{\prime}$ for twisted systems are obtained averaging over all possible commensurate stacking configurations given by the integral

[TABLE]

that in the first harmonic approximation can be approximated by taking the average of the sum over the three symmetric stacking configurations $s={\rm AA,\,AB,\,BA}$ at their respective equilibrium interlayer distances. jung2014 Because tunneling between interlayer sublattices in graphene on graphene vanish at symmetric stackings $s$ when they are not vertically aligned we have

[TABLE]

Using Eqs. (6-7) at zero pressure and using the EXX+RPA equilibrium distances for each stacking leconte2017 we get $\omega=0.12$ eV for $c_{\rm AB}$ such that $t_{1}=3\omega=0.36$ eV, close to the LDA interlayer coupling in Bernal BG, and $\omega^{\prime}=0.098$ eV for $c_{\rm AA}$ which is incidentally rather close to the interlayer tunneling from explicit integrations in $\vec{d}$ in Ref. [jung2014, ]. The effects of atomic relaxation in the moire patterns can have non-negligible effects in the details of the electronic structure for both intralayer potentials and interlayer coupling that can be captured with higher order harmonics in the moire $G$-vectors. jung2015 It was also noted that in tBG unequal interlayer coupling values $\omega\neq\omega^{\prime}$ enhances the gap between the low energy flatband and its neighboring higher energy band. koshino_wannier The band structures for type I AB/AB structures near $0^{\circ}$ alignment calculated for the minimal tBBG model are shown in Fig. 2 both for rigid unrelaxed $\omega=\omega^{\prime}$ and out of plane relaxed interlayer tunneling values $\omega\neq\omega^{\prime}$. In the case of tBG it was noted that the magic angle follows the $\theta\propto\omega/\upsilon_{\rm F}$ proportionality chittari2019 where the magic angles grow with increasing interlayer coupling strength. Hence, an enhanced $\upsilon_{\rm F}=1.0\times 10^{6}$ m/s (or $t_{0}=-3.1$ eV), together with ab initio tunneling $\omega=0.12$ eV, $\omega^{\prime}=0.098$ eV from Eqs. (6-7) leads to similar low energy bands as the LDA $\upsilon_{\rm F}=0.84\times 10^{6}$ m/s (or $t_{0}=-2.6$ eV) combined with weaker $\omega=0.1$ eV and $\omega^{\prime}=0.08$ eV values. In the calculations to follow we use the enhanced Fermi velocity of $\upsilon_{\rm F}=1.0\times 10^{6}$ m/s together with the ab initio interlayer tunneling $\omega$ and $\omega^{\prime}$ based on Eqs. (6-7) compatible with EXX+RPA equilibrium distances. Our calculations have used a configuration space with variable cutoff in momentum space of a radius of up to $6G_{1}=24\pi\theta/(\sqrt{3}a)$ using Hamiltonian matrices with sizes as large as 676$\times$676 such that $\theta\gtrsim\omega/(12\pi\left|t_{0}\right|)$ to obtain converged results in the limit of small $\theta$ and large $\omega$.

An important distinctive feature of tBBG with respect to tBG is that we have an additional control knob to change the electronic structure through a perpendicular external electric field that modifies the interlayer potential $\Delta_{i}$ values in Eq. (15). The potential drops introduced by an external electric field could be modeled through the parameter set $\Delta_{1}=-\Delta_{4}$, $\Delta_{2}=-\Delta_{3}$, redefined as $\Delta_{1}=(\Delta+\Delta^{\prime})/2$, $\Delta_{2}=(-\Delta+\Delta^{\prime})/2$ in terms of $2\Delta$, the interlayer potential difference within each BG, and $2\Delta^{\prime}$ the potential difference between the BG. We will use the relation $\Delta^{\prime}=2\Delta$ to introduce equal interlayer potential drops of $\Delta^{\prime}$ between the consecutive layers to model the effects of an electric field. A qualitatively different interlayer potential configuration consists in having the electric fields point in opposite directions at each BG. This could be done by grounding the tBBG device and using equipotential top/bottom gates to accumulate charges of the same sign at the outer layers. The potential distribution of this case can be modeled by $\Delta_{1}=-\Delta_{2}=-\Delta_{3}=\Delta_{4}=\Delta$, where the reversal of the relative mass sign between the top and bottom BG can modify the topology of the resulting flatbands. Distinct band topologies are thus expected for a system subject to a perpendicular electric field near 180*∘* (or $60^{\circ}$) alignment where for the same mass sign the chirality of the bands at $K$ are reversed.

III Flatbands as a function of twist angle, electric fields and pressure

In the following we discuss the electronic structure results and the moire flatband bandwidth in tBBG as a function of system parameters such as twist angle $\theta$, the inter-layer coupling $\omega$ tunable by pressure, and interlayer potential differences due to an electric field, in search of the optimal conditions for finding isolated low energy flatbands near the Fermi level. The interlayer potential differences due to an electric field are modeled combining potential differences between the layers within each BG ($\Delta$) and inter-BG potential offset ($\Delta^{\prime}$) which are related to each other through $\Delta^{\prime}=2\Delta$ that can lift the degeneracy of the flatbands.

In the band structures resulting from the minimal model shown in Fig. 2 for a twist angle of $\theta=1.5^{\circ}$, we can still distinguish features of the original BG band structure at the Dirac cones, and can identify the band structure near the magic angle $\theta=1.05^{\circ}$ of tBG. The first important observation is that the overall bandwidths of the low energy bands in tBBG are almost half of those corresponding to tBG for a similar range of $\theta$ and $\omega$ parameter values, suggesting that the tBBG system is generally more suitable for the generation of narrow bandwidth flatbands than in tBG. This is shown in Fig. 3 where we represent the bandwidth as a function of twist angle $\theta$ for fixed values of interlayer coupling $\omega$ and as a function of $\omega$ for fixed $\theta$ values. From the bandwidth versus $\theta$ dependence we can observe that the bandwidths remain below 10 meV for every twist angle below and around the first magic angle. Likewise, the bounce off of the bandwidth for increasing $\omega$ past the critical value at the first magic angle have maxima that are roughly half of those seen in tBG. chittari2019 We thus expect that in tBBG the twist angle control does not need to be as precise as in tBG to maintain a moderately narrow bandwidth on the order of $\sim$10 meV for twist angles smaller than $\sim 1^{\circ}$. Inclusion of remote hopping terms results in band gaps near charge neutrality for sufficiently large twist angles and widening of the bandwidths with respect to the minimal model, as shown in Fig. 4. The trigonal warping term in BG generates several band touching points in the vicinity of the Dirac point, and introduces particle-hole symmetry breaking of sufficient relevance especially when we apply an external electric field. In Fig. 4 we compare the band structures of the minimal model and the more complete model that includes the remote hopping parameters, both for type I near 0*∘* and type II near $180^{\circ}$ alignments. For the complete Hamiltonian we include the remote hopping terms $t_{3}$ and $t_{4}$, the site potential offset $\delta$ between the high energy dimer sites, and interlayer coupling $\omega$ and $\omega^{\prime}$ in Eqs. (6-7) evaluated at the equilibrium out of plane relaxed lattice distances.

The second important observation is the tunability of the primary and secondary band gaps accompanied by a variation in bandwidth due to an electric field, as illustrated in Fig. 4. Even for a twist angle $\sim$1.5*∘* considerably greater than the minimal model magic angle of $\sim 1^{\circ}$ the low energy bands can remain isolated thanks to the primary $\delta_{p}$ and secondary $\delta_{s}$ gaps. These gap values are obtained from the difference between the maximum (minimum) energy of the flatband and the minimum (maximum) energy of the neighboring higher (lower) band resulting in positive values when there is a gap and giving negative values when there is band overlap. Application of external fields contributes in changing the bandwidth and a relatively narrow bandwidth on the order of $\sim$25 meV or smaller is achievable for moderate electric fields that introduce interlayer potential differences typically of a few tens of meV. The bandwidths were obtained from the difference between the maximum and minimum band energy for a given band within the mBZ. Band gaps within each BG layer ($\Delta$) and band offsets between BG layers ($\Delta^{\prime}$) are simultaneously present when a perpendicular electric field is applied in the system. It is found that generally $\Delta^{\prime}$ contributes in widening the bandwidth of the low energy bands. An important factor for the onset of the interaction driven ordering is the isolation of the low energy bands that can be quantified from the primary gap $\delta_{p}$ near charge neutrality and the secondary gap $\delta_{s}$ near the $\Gamma$ point of mBZ, since greater band isolation reduces screening and strengthens the Coulomb interactions. Hence, the parameter space most likely to observe Coulomb driven ordered phases should have simultaneously smaller bandwidths $W$ and larger isolation gaps $\delta_{p}$ and $\delta_{s}$. We can estimate the ratio of $U_{\rm eff}/W$ from the effective 3D screened Coulomb potential

[TABLE]

where the moire length is $l_{M}=a/\theta$, and the Debye length $\lambda_{\rm D}=2\varepsilon_{0}/e^{2}D(\delta_{p},\delta_{s})$ uses the 2D density of states $D(\delta_{p},\delta_{s})=4\,\left(\left|\delta_{p}\right|u(-\delta_{p})+\left|\delta_{s}\right|u(-\delta_{s})\right)/(W^{2}A_{M})$ that assumes a value proportional to the band overlap ratio $\delta_{p/s}/W$ when $\delta_{p/s}<0$, where $u(x)$ is the heaviside step function, $\varepsilon_{r}=4$ and we counted four valley-spin degenerate electrons per moire unit cell area $A_{M}=\sqrt{3}\,l_{M}^{2}/2$ for each filled moire band. This ratio in Eq. (8) is used to find the parameter space region of twist angle and interlayer electric field with narrow bands and strong effective Coulomb interactions, see Fig. 15. While the parameter region near $\theta\sim 0.5^{\circ}$ show largest $U_{\rm eff}/W$ ratios compared to $\theta\gtrsim 1^{\circ}$ due to the greater flatness of the bands at small twist angles, the closer proximity of the neighboring energy bands in this regime may enhance the Coulomb screening in a way that is not captured in the screening model we have used. We expect that the electron-hole asymmetry resulting from the intrinsic asymmetry of BG with remote hopping terms can be further affected by the coupling with the substrate, for example by aligning BG with a hexagonal boron nitride substrate which in a heterojunction with single layer graphene it introduced a strong intrinsic particle-hole asymmetry in the electronic structure. dasilva-particle-hole

A phase diagram of the bandwidth, the primary gap, and secondary gap as a function of twist angle and external field strength is presented in Fig. 6. For simplicity we have assumed that the potential differences between contiguous layers are given by $\Delta^{\prime}=2\Delta$ and are the same for a given electric field, although the precise interlayer potential differences will depend on the screening between the layers. This phase diagram of $U_{\rm eff}/W$ ratios, bandwidths, primary, and secondary gaps both for electrons and holes illustrates the parameter space where the likelihood for finding ordered phases is higher. As a general trend we find a non-negligible asymmetry between electrons and holes, and the possibility of finding states with primary and secondary gaps both at small and large twist angles of $\sim 0.5^{\circ}$ and $\sim 1.5^{\circ}$ for sufficiently strong electric fields. Due to the almost linear increase of the low energy bandwidth with increasing twist angle there is a tradeoff between band isolation easier at larger twist angles and bandwidth increase to find the optimum parameter space where the Coulomb interaction effects will be strongest.

Application of pressure is also a useful control knob to tune the electronic structure of 2D materials. yankowitz1 ; Yankowitz ; kaxiras2018 ; chittari2019 In the bandwidth and gaps phase diagram as a function of twist angle and pressure shown in Fig. 7 both for $\Delta^{(\prime)}=0$ and $\Delta^{(\prime)}\neq 0$ we can observe that the enhancement of $\omega$ through pressure generally compresses the bandwidth of the flatbands. In the minimal model of tBBG the application of pressure leads to a magic angle line

[TABLE]

whose coefficient values of $C_{1}=27.5$, $C_{2}=10.5$, and $C_{3}=5.6$ agree within 10% with the results obtained for tBG, chittari2019 in keeping with the expected band scaling behavior proportional to $\alpha=\omega/(\theta\upsilon_{\rm F}K_{D})$, bistritzer2011 see supplemental materials. supplemental This observation suggests that increase of $\omega$, e.g. through external pressure $P$, yankowitz1 should allow to achieve narrow bandwidth features for larger twist angle $\theta$. This behavior holds both for the minimal model and when remote hopping terms are considered, although the remote hopping terms tend to broaden the bandwidth of the minimal model, see the supplemental information for a comparison of the phase diagrams. supplemental The phase diagram also illustrates how pressure can be used to enhance the secondary gaps both for valence and conduction bands, but it can at the same time suppress the primary gap. The phase diagram is modified substantially in the presence of interlayer potential differences $\Delta^{(\prime)}\neq 0$ triggered by a perpendicular electric field, generally widening the bandwidth in the phase diagram and shifting the weights between the primary and secondary gaps. The relationship between $P$ and $\omega$ values in the relevant range of pressures between 0 to $\sim$15 GPa if fitted by a second order polynomial and its positive root

[TABLE]

The fitting parameters are listed in Table 1 and they are found to be valid over a wide range of pressures stretching up to $\sim$30 GPa and also for negative values down to $\sim-1$ GPa. To obtain the above fitting parameters we have used the relationship between $P$ and interlayer separation $c$ for every stacking of Ref. [leconte2017, ], and the calculations of $\omega$ versus $c$ as detailed in Ref. [jung2014, ]. The explicit fitting functions for these quantities are presented in the supplemental material. supplemental This approach is different to that in Ref. [chittari2019, ] where the total pressure was obtained from the average of the local pressure values at the same interlayer distance $c$ for every stacking.

IV Topological flatbands in twisted bi-bilayer graphene

The presence of moiré superlattices give rise to avoided gaps between the bands in the moire mini Brilloun zone (mBZ) allowing them to have a well defined valley Chern number in a wide class of twisted multilayer systems, todadri ; senthil1 ; zaletel transition metal dichalcogenides, wumacdonald trilayer graphene on hexagonal boron nitride (TG/BN), prljung and for a variety of twisted gapped Dirac materials. javvaji For our tBBG system the possibility of opening a band gap $\delta_{p}$ near charge neutrality through an electric field together with the opening of a secondary gap $\delta_{s}$ with the higher energy bands makes tBBG an interesting platform to engineer flatbands with well defined valley Chern numbers that are tunable through an electric field like in TG/BN. todadri ; prljung The valley Chern number phase diagram in Fig. 9 as a function of interlayer potential differences and twist angle indicates the range of Chern numbers $C_{\upsilon}=0,\,\pm 1,\,\pm 2,\,\pm 3,\,\pm 4$ expected in these systems. The valley Chern numbers were calculated through

[TABLE]

by integrating in the moire Brillouin zone for each valley the Berry curvature for the $n^{th}$ band through rmp_berry

[TABLE]

where for every $k$-point we take sums through all the neighboring $n^{\prime}$ bands, the $|u_{n}\rangle$ are the moiré superlattice Bloch states and $E_{n}$ are the eigenvalues. There are clear qualitative differences between the valley Chern numbers for twist angles near $0^{\circ}$ for type-I AB/AB and those near $180^{\circ}$ for type-II AB/BA alignments. In the first case the valley Chern numbers between valence and conduction bands are generally opposite in value adding up to a total of zero, while in the second case they are the same number for both valence and conduction bands. These differences are naturally expected if we consider that the chirality of the massive bands at the top and bottom layers that couple to each other are interchanged depending on the alignment. Quantitative modifications in the valley Chern numbers are observed when we compare the phase diagrams of the minimal and remote hopping parameter models as a function of $\Delta$ where clear differences in particular for small $\Delta\sim 10$ meV region is observed, comparable in magnitude with the band distortions introduced by the remote hopping terms. The vicinity of $\theta\sim 1^{\circ}$ and low electric field $\Delta$ have spots where the valley Chern numbers for the valence and conduction bands differ in magnitude. For sufficiently large values of $\Delta$ the valley Chern numbers of the lowest energy flatbands agree for the minimal and remote hopping parameter models in a large parameter space of twist angles as illustrated in Fig. 10 by showing the band structure and the Berry curvatures used in the valley Chern number calculations. Hence, the topological properties of the bands remain overall relatively robust to small perturbations to the band structures introduced by the remote hopping parameters. The values of the low energy flatband Chern numbers and those of the neighboring higher energy bands are gathered in Table 2. When the valley Chern numbers between the minimal and remote hopping term models differ are comparable in a large parameter space.

V Summary and conclusions

We have extended the Bistritzer-MacDonald continuum model of twisted bilayer graphene (tBG) to investigate the electronic structure of twisted bi-bilayer graphene (tBBG) as a function of twist angle $\theta$, electric fields, and the pressure dependent inter-layer coupling. We have considered both the minimal model and also the effects of the remote hopping terms in the band structure calculation of bilayer graphene. The calculated bandwidth phase diagram for the low energy bands shows that the bandwiths are roughly a factor two narrower than those in twisted bilayer graphene (tBG) indicating that tBBG should be more forgiving in the twist angle precision required to access the strongly interacting regime, and for this reason we expect that the narrow band features in tBBG will be observed more simply than in tBG. The possibility of applying a perpendicular electric field is an interesting control knob that allows to enhance the separation between the flatbands and also influences the gaps with the higher energy bands favoring a more effective band isolation. At the same time, we find that interlayer potential differences can widen the bandwidths near the first magic angle of the minimal model and smoothen the bandwidth variation to give a continuous range of angles where the bandwidths are narrow. Within the minimal model, the bandwidth phase diagram for zero interlayer potential difference and small perturbations thereof is found to be closely similar to the case of tBG, maintaining the same linear dependence between $\theta$ and the interlayer coupling $\omega$ for the magic angles, and the inverse proportionality to the Fermi velocity of the graphene layers. Our calculations show that bi-bilayer graphene under a perpendicular electric field can host robust ordered phases for twist angles in the vicinity of $\sim 0.6^{\circ}$ and near $\sim 1.5^{\circ}$, with the parameter space for the conduction bands being generally favored over those of the valence bands. With proper electric fields we expect that practically all angles spanning the range between $0.4^{\circ}-1.6^{\circ}$ could host ordered phases. Application of pressure can also enhance the isolation of the bands when used in combination with appropriate electric fields. We have related the interlayer tunneling with external pressure through stacking dependent interlayer coupling parameters $\omega$ and $\omega^{\prime}$ compatible with the EXX+RPA interlayer potentials to capture the interdependence of corrugation and interlayer tunneling in a self-contained manner. A more detailed study that combines the effects in-plane moire strains will be addressed elsewhere. Comparison between electronic structure between the type I near $0^{\circ}$ aligned systems and type II $180^{\circ}$ aligned bi-bilayers indicate that $0^{\circ}$ aligned systems generally require weaker electric fields to achieve optimal flatband systems prone to interaction driven ordered phases, and give rise to distinct valley Chern number phase diagrams.

Acknowledgements.

We acknowledge helpful discussions with S. Kahn, F. Wu and M. Koshino. This work was supoprted by Samsung Science and Technology Foundation under project no. SSTF-BA1802-06 for J.J. Financial support is acknowledged for N.R.C. from the Korean National Science Foundation through grant number NRF-2016R1A2B4010105, for B.L.C. by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2018R1A6A1A06024977) and NRF-2017R1D1A1B03035932. This research was supported in part by the US National Science Foundation under Grant No. NSF PHY-1748958.

Note added.$-$ Related experiments and theory are Refs. [philipkim, ; guangyu, ; pjarillo, ] and Refs. [choi, ; leey, ; koshi, ; liu, ].

VI Supplemental material

Here we provide additional figures to supplement the information in the main text for (1) the pressure dependent bandwidth and gaps phase diagram, (2) the parametrization of the interlayer distance dependent pressure and tunneling for the $P$ versus $\omega$ relationship, and (3) the low energy bands formulation of the Hamiltonian.

VI.1 Effect of interlayer coupling strength

The pressure dependent phase diagrams for valence (v), conduction (c) flatband bandwidth and primary $\delta_{p}$, secondary $\delta_{s}$ gaps between the low energy flatband with the surrounding bands with the variation of interlayer potential difference $\Delta^{\prime}$ for the 0*∘* (Type-I) and 60*∘* (Type-II) aligned BG/BG are shown to understand the effect of interlayer potential difference. We discuss the results obtained within the minimal model for two types of stacking in Fig. 11 and Fig. 12. The effects of remote hopping terms in the Hamiltonian are discussed in Fig. 13 and  14

VI.2 Coulomb interaction potential

Coulomb interaction driven instabilities arise in tBBG due to narrow bandwidth. The strong correlations are evident from the $U/W$ ratio, and here we take ratio of the effective coulomb interaction potential without screening:

[TABLE]

The Fig. 15 is obtained for the parameter space of $\theta$ and $\Delta$ .

VI.3 P vs $\omega$ relations

Here we present the various relations that are using in the main-text, Pressure versus interlayer distance is shown in Fig. 16, $\omega$ versus interlayer distance is given in Fig. 17. The relationship between $\omega^{\prime}$ and $\omega$ is shown in Fig. 18 for RPA, and in Fig. 19 is for the LDA.

VI.4 Low energy effective model

The low energy effective model formulation results in simpler expressions for the Hamiltonian and often allows to obtain solutions with compact analytical forms. In the main text we have presented the best results for the full bands model but here we present for sake of completeness some analysis of the low energy bi-bilayer models subject to an electric field. The discrepancies between the results indicates that the low energy model is not useful to guide experimental efforts with predictive accuracy. The low energy model hamiltonian for the twisted bi-bilayer with top ($+$) and bottom ($-$) bilayers is given by:

[TABLE]

where

[TABLE]

and the site potentials for each layer are given by the parameter $\Delta_{i}$ where $i=1,\,2$ from top and bottom layers in each Bernal stacked bilayer. The tunneling matrix are as below:

[TABLE]

Expression for bandwidth

Type-I twisted BGBG (AB stacked BG on AB stacked BG) Hamiltonian can be written as,

[TABLE]

where,

[TABLE]

[TABLE]

[TABLE]

where, $\xi=\pm 1$, for $K,\,K^{\prime}$ valleys, $\Psi_{b}=\{\phi_{B_{1}^{b}}\,,\phi_{A_{1}^{b}}\,,\phi_{B_{2}^{b}}\,,\phi_{A_{2}^{b}}\}^{\dagger}$, $\Psi_{t}=\{\phi_{A_{1}^{t}}\,,\phi_{B_{1}^{t}}\,,\phi_{A_{2}^{t}}\,,\phi_{B_{2}^{t}}\}^{\dagger}$ and $t_{\perp}=3\omega$ is the vertical hopping strength between the dimmer cites in bilayers and $T_{j}$ is the interlayer hopping matrix valid in the moderate twist angle limit. For analytical discussions, we have derived the Effective two band model Hamiltonian of tBGBG-I by integrating out the higher energy dimmer sites using the perturbation method,

[TABLE]

where, $q\,=\,\left[k^{2}-(\Delta K/2)^{2}\right]\,;\,\,|\Delta K|=K_{\theta/2}-K_{-\theta/2}=\dfrac{4\pi}{3a}2\,sin(\theta/2)$ is the momentum shift between the Dirac cones of top and bottom bilayers due to rotation and $a=\sqrt{3}a_{cc}$. For simplicity we have taken $\upsilon_{4}=0$. The eigen energies of the above effective Hamiltonian are given by,

[TABLE]

where, $q_{m}=|\Delta K|/2$ is the maximum value of $q$ at the Dirac points. The difference between the maximum and minimum values of the band energy is defined as the Bandwidth of the corresponding band. Here we have calculated the bandwidth,

[TABLE]

From the above expression one can observe that the band width is decreasing nonlinearly with the interlayer coupling $\omega$ and external electric field is reducing the band width and trigonal warping strength is enhancing it within the limit of moderate electric field.

VI.5 Magic twist angle

By equating the band width $(W)$ to zero we can get the magic angle, where it leads the formation of flatbands.

[TABLE]