Moir\'e Flat Bands in Twisted Double Bilayer Graphene

TL;DR

This paper explores twisted double bilayer graphene, revealing a narrow magic angle region with flat bands and establishing a fundamental link to twisted bilayer graphene, advancing understanding of moiré superlattice physics.

Contribution

The study introduces a continuum model for TDBG, demonstrating a fundamental connection to TBG and identifying intrinsic flat bands near the magic angle without external fields.

Findings

Identification of a narrow flat band region around 1.3° twist angle.

Discovery of an intrinsic energy gap at charge neutrality due to symmetric polarization.

Establishment of a mathematical link between TDBG and TBG models.

Abstract

We investigate twisted double bilayer graphene (TDBG), a four-layer system composed of two AB-stacked graphene bilayers rotated with respect to each other by a small angle. Our ab initio band structure calculations reveal a considerable energy gap at the charge point neutrality that we assign to the intrinsic symmetric polarization (ISP). We then introduce the ISP effect into the tight-binding parameterization and perform calculations on TDBG models that include lattice relaxation effects down to very small twist angles. We identify a narrow region around the magic angle $\theta^\circ = 1.3^{\circ}$ characterized by a manifold of remarkably flat bands gapped out from other states even without external electric fields. To understand the fundamental origin of the magic angle in TDBG, we construct a continuum model that points to a hidden mathematical link to the twisted bilayer graphene (TBG) model, thus indicating that the band flattening is a fundamental feature of TDBG, and is not a result of external fields.

Tables

Table 1. Table 1: Comparison of the band gaps (in meV) for TDBG models characterized by different commensurate parameters m 𝑚 m twist angles θ 𝜃 \theta calculated using DFT, conventional tight-binding model (TB) and the TB model that includes the intrinsic symmetric polarization (TB+ISP).

$m$	$1$	$2$	$3$	$4$	$5$	$6$
$\theta$ (∘)	$21.79$	$13.17$	$9.43$	$7.34$	$6.01$	$5.09$
DFT	$30.08$	$34.41$	$31.45$	$29.88$	$23.03$	$25.26$
TB	$1.16$	$3.31$	$2.54$	$1.41$	$0.02$	$1.58$
TB+ISP	$30.94$	$33.06$	$32.34$	$31.20$	$29.75$	$28.05$

Equations

$${\mathcal{H}}_{\text{AB}}=\left(\begin{array}{cc}{v}_F\mathbf{\sigma }\cdot \mathbf{k}& T_0\\ T_0^{†}& v_F\mathbf{\sigma }\cdot \mathbf{k}\end{array}\right),T_0=\left(\begin{array}{cc}{\gamma }_{4}k& {\gamma }_3{k}^{\ast }\\ {\gamma }_1& {\gamma }_{4}k\end{array}\right),$$

$$T_{\theta }(\mathbf{r})=\sum _{n=1,2,3}{T}_n{e}^{-i{\mathbf{q}}_n\mathbf{r}},$$

$$T_n=e^{-i{\mathcal{G}}_{\theta }^{(n)}\mathbf{d}}{\hat{\Omega }}_{\varphi }^{n-1}\left(\begin{array}{cc}{w}_{AA}& w_{AB}\\ w_{AB}& w_{AA}\end{array}\right){\hat{\Omega }}_{\varphi }^{1-n},$$

$$\mathcal{H}=\left(\begin{array}{cccc}{\mathcal{H}}_1& T_0& 0& 0\\ T_0^{†}& {\mathcal{H}}_2& T_{\theta }(\mathbf{r})& 0\\ 0& T_{\theta }^{†}(\mathbf{r})& {\mathcal{H}}_3& T_0\\ 0& 0& T_0^{†}& {\mathcal{H}}_4\end{array}\right),$$

$$\mathcal{H}=\left(\begin{array}{cc}0& D^{†}(\mathbf{r})\\ D(\mathbf{r})& 0\end{array}\right),D(\mathbf{r})=\left(\begin{array}{cc}-2i\bar{\partial }& \alpha {\hat{A}}_1(\mathbf{r})\\ \alpha {\hat{A}}_2(\mathbf{r})& -2i\bar{\partial }\end{array}\right),$$

Full text

Twisted bilayer graphene (TBG) has recently attracted considerable attention following the discovery of the correlated insulator and superconducting phases Cao et al. (2018a, b); Yankowitz et al. (2019), if tuned to the so-called magic angle of 1.1∘, at which almost dispersionless (flat) electronic bands emerge near the charge neutrality point Trambly de Laissardière et al. (2010); Bistritzer and MacDonald (2011); San-Jose et al. (2012); Gargiulo and Yazyev (2017); Tarnopolsky et al. (2019); Carr et al. (2019). The appearance of flat bands at magic angles in TBG is not just a coincidence in material properties engineering, but a fundamental feature of the TBG-like Hamiltonians Tarnopolsky et al. (2019); Khalaf et al. (2019). In such models, the flat band in the electronic spectrum appears if two Dirac cones are brought sufficiently close in reciprocal space via the moiré interlayer potential, so they hybridize Lopes dos Santos et al. (2007) and degenerate into remarkably flat bands Bistritzer and MacDonald (2011); Tarnopolsky et al. (2019). In this regard, one can consider an extension of this picture to the hybridization of touching (gapless) pairs of parabolic bands that are primarily flatter at the touching point. Such scenario can be realized, for example, in the Bernal-stacked (AB-stacked) bilayer graphene characterized by touching parabolic bands at the $K$ and $K^{\prime}$ points in the Brillouin zone Castro Neto et al. (2009). This idea has proved to be fruitful as witnessed by the recent experimental reports of unconventional superconductivity and spin-polarized insulating phases in twisted double bilayer graphene (TDBG) at twist angle $1.24^{\circ}$ Liu et al. (2019a); Cao et al. (2019); Shen et al. (2019); Lee et al. (2019). Given the current progress in manufacturing TDBG, and strong sensitivity of the superconducting phase to the factors that modify the band structure (electric displacement and in-plane magnetic fields), one of the key questions is understanding the very nature of band flattening in TDBG and the evolution of the band structure upon changing the twist angle.

In this Letter, we address the moiré flat bands in TDBG using ab initio and tight-binding approximation calculations as well as an effective continuum model. Our DFT calculations performed on TDBG models with twist angles down to 2.5∘ reveal the presence of band gaps as large as 35 meV. These persistent band gaps are explained by the intrinsic polarization of individual TBG components due to the proximity with the complementary TBG counterpart. By including this intrinsic symmetric polarization (ISP), an effect that has not been considered previously, into the tight-binging (TB) model we are able to perform accurate band structure calculations on TDBG models that include lattice relaxation effects down to very small twist angles. In addition gap opening, the ISP flatten the electronic bands near the charge neutrality point. We also find a well-defined magic angle of $1.3^{\circ}$, at which both electron and hole gaps are maximized. Outside the narrow range of twist angles around $1.3^{\circ}$ the bands become highly dispersive, in contrast to the case of twisted transition metal dichalcogenides bilayers Wu et al. (2019); Naik and Jain (2018) where the magic angle is hard to define. To understand the origin of magic angle in TDBG, we further construct an effective continuum model describing the behavior at large moiré periodicities. This allows us to understand the mechanism of band flattening in TDBG through the exact mapping onto the continuum TBG Hamiltonian Tarnopolsky et al. (2019) if the particle-hole asymmetry is excluded.

Twisted double bilayer graphene (TDBG) is a very special configuration among the large family of graphene multilayers, which has gained a well-deserved attention due to the characteristic moiré physics (see e.g. Refs Khalaf et al., 2019; Cea et al., 2019; Mora et al., 2019; Zhang et al., 2019; Koshino, 2019; Choi and Choi, 2019). In the most general setting, the geometric overlay of four lattices leads to a sophisticated pattern of three interfering moiré superlattices, which may or may not be periodic and commensurate. A more fruitful and simpler situation for both experiment and theory is the four-layer graphene configuration with only one twist parameter. This can be achieved in two ways, either twisting an AB-stacked bilayer on another AB-stacked bilayer by angle $\theta$, or working in the twist-alternate setting Khalaf et al. (2019). We focus on the former situation shown in Fig. 1a that addresses recent experimental works Liu et al. (2019a); Cao et al. (2019); Shen et al. (2019); Lee et al. (2019). The corresponding mini Brillouin zone of the moiré superlattice (Fig. 1b) is the same as in the case of TBG. Below, we consider TDBG models that correspond to the principal moiré branch with a discrete set of periods $\lambda(\theta_{m})=a/2\sin(\theta_{m}/2)$, see Refs. Lopes dos Santos et al., 2007; Shallcross et al., 2008a, b. For every commensurate twist angle $\theta_{m}$, the supercell periodicity vectors are defined as $\mathbf{t_{1}}=m\mathbf{a_{1}}+(m+1)\mathbf{a_{2}}$ and $\mathbf{t_{2}}=-(m+1)\mathbf{a_{1}}+(2m+1)\mathbf{a_{2}}$, where $a_{1,2}=a(\sqrt{3}/2,\pm 1/2)$ are the lattice vectors of graphene with lattice constant $a=2.46$ Å and $m$ is commensuration condition parameter, $\cos\theta_{m}=(3m^{2}+3m+{1}/{2})/(3m^{2}+3m+1)$.

DFT calculations.— We performed DFT calculations 111Our DFT calculations have been performed using the Vienna ab initio simulation package (VASP) Kresse and Furthmüller (1996) within the GGA approximation. The cutoff energy of $400$ eV was chosen for the plane wave basis. The Brillouin zone was sampled using meshes of $N\times N$ k points, where $N=36,24,12,9,6,6$ for the TDBG models described by $m=1,..,6$, respectively. The total energy convergence threshold was set to $10^{-7}$ eV. on TDBG models characterized by commensuration parameters $m$ up to $m=13$, which corresponds to $\theta=2.45^{\circ}$ and includes 4376 atoms per supercell. Calculations on larger models for smaller twist angles are unaccessible for computational cost reasons. Figure 2 shows the band structures of $m=6$ ($\theta=5.09^{\circ}$) and $m=13$ ($\theta=2.45^{\circ}$) TDBG models. One can immediately notice a gap opening where the two parabolic bands characteristic of AB-stacked BLG are supposed to touch. Table 1 further confirms the systematic presence of few tens meV band gaps in $m=1,..,6$ TDBG models when DFT calculations are considered. Interestingly, these band gaps are not qualitatively reproduced by the conventional tight-binding (TB) model Trambly de Laissardière et al. (2012) widely used for studying twisted multilayer graphene, even though gap opening in a certain range of twist angles has already been pointed out Choi and Choi (2019). This allows us to conclude that the mechanism responsible for this gap opening is not accounted for by the standard TB model. Indeed, each of the AB-stacked BLG counterparts in TDBG is placed in an asymmetric environment with one of the graphene layers facing vacuum and another facing adjacent graphene layer. This produces a slight potential difference between the two graphene layers which opens a band gap in otherwise gapless AB-stacked BLG. Below, we will refer to this effect as the intrinsic symmetric polarization (ISP). Fig. 1a depicts it as effective electric fields that have opposite orientation in the two BLG subsystems. We find that introducing a single universal value of 9 mV/Å for the electric field into the tight-binding model allows to reproduce the DFT results for all investigated TDBG models (TB+ISP in Table 1 and Fig. 2).

Flat bands in the tight-binding model with lattice relaxation effects.— In order to extend the scope of our models to twist angles in the range relevant to the flat-band physics, we perform the tight-binding model calculations. Our tight-binding model 222We consider a tight-binding model Hamiltonian for carbon atom $p_{z}$ orbitals $H=\sum_{i\neq j}V_{ij}a^{{\dagger}}_{i}a_{j}$, where $a^{{\dagger}}_{i},a_{j}$ are the creation and annihilation operators. The $\sigma$-type and $\pi$-type Slater-Koster parameters contribute to $V_{ij}=V_{pp\pi}\sin^{2}\theta+V_{pp\sigma}\cos^{2}\theta$, where $\theta$ is the angle between the orbital axes and $\mathbf{R}_{ij}=\mathbf{R}_{i}-\mathbf{R}_{j}$ connects the two orbital centers Slater and Koster (1954). For $\theta=\pi/2$ ($\theta=0$) which corresponds to the pair of atoms in the same layer (the pair of atoms on top of each other) $V_{ij}=V_{pp\pi}$ ($V_{ij}=V_{pp\sigma}$). The Slater-Koster parameters depend on the distance $r$ between two orbitals as $V_{pp\pi}(r)=V^{0}_{pp\pi}e^{q_{\pi}(1-r/a_{\pi})}F_{c}(r)$ and $V_{pp\sigma}(r)=V^{0}_{pp\sigma}e^{q_{\sigma}(1-r/a_{\sigma})}F_{c}(r)$, where $a_{\pi}$ is the first nearest neighbor distance in the plane and $V^{0}_{pp\pi}$ is the corresponding coupling value. In our work we assume $V^{0}_{pp\pi}=-2.81$ eV, which is larger than the conventional value of $-$2.7 eV Trambly de Laissardière et al. (2012), in order to provide a better fit to our DFT calculations. The second nearest neighbour coupling of $0.1\times V^{0}_{pp\pi}$ fixes the value of $q_{\pi}/a_{\pi}$. Here, $a_{\sigma}$ is the interlayer distance in AB-stacking BLG. In our relaxed structures the interlayer distance between the surface and inner layers is $3.364$ Å and between two inner layers is $3.410$ Å. Choosing the same exponential decay for both hopping parameters gives $q_{\sigma}/a_{\sigma}=q_{\pi}/a_{\pi}$. $F_{c}$ is a smooth cutoff function that takes into account the distance between orbitals $F_{c}(r)=(1+e^{(r-r_{c})/l_{c}})^{-1}$ Trambly de Laissardière et al. (2012), where $l_{c}=0.265$ Å, $r_{c}=2.5a=6.165$ Å. For $r\ll r_{c}$, $F_{c}(r)\simeq 1$ and for $r\gg r_{c}$, $F_{c}(r)\simeq 0$. The on-site energy of $p_{z}$ orbitals is set to $\epsilon_{i}=-0.78$ eV to adjust the reference. All tight-binding parameters being defined, the band structures were calculated using WannierTools Wu et al. (2018). is based on the one described in Ref. Trambly de Laissardière et al. (2012) with the ISP being taken into account as described above. Figure 2 evidences excellent agreement with the DFT results. Furthermore, we explicitly include lattice relaxation which becomes crucial at small twist angles using the methodology established by us previously Gargiulo and Yazyev (2017). Important details are summarized in Ref. 333 The lattice relaxation of TDBG models was performed in the classical force-field approach using the LAMMPS package LAM . The second generation REBO potential Brenner et al. (2002) and the Kolmogorov-Crespi (KC) potential were used to describe the intra-layer and inter-layer interactions, respectively. The initial inter-layer distance was set as 3.35 Å and the $D_{3}$ symmetry was preserved during the relaxation. The relaxation was performed using the FIRE algorithm Bitzek et al. (2006) until the total force acting on each atom becomes less than 10*-6* eV/atom. Details of the relaxation results are presented in Ref. SI .

Figure 3a shows the band structure of TDBG at twist angle $\theta=2.45^{\circ}$, which is above the magic angle $\theta^{\star}$=1.3∘, for both rigidly twisted and relaxed models. In this case, the lattice relaxation increases the band gap $\Delta_{gap}$ produced by the ISP. Further decrease of the twist angle closes this band gap and results in a manifold of very narrow bands separated from the rest of bands by gaps $\Delta_{e}$ and $\Delta_{h}$ (Fig. 3b). Importantly, taking into account the relaxation effects increases these gaps dramatically, resulting in close-to-maximum values of $\Delta_{e}=37.6$ meV and $\Delta_{h}=38.8$ meV at $\theta$=1.3∘ (Fig. 3b). At the same time, the width of the flat-band manifold $E_{W}=11$ meV is achieved. Figures 3c-g summarize the dependence of crucial band structure parameters $\Delta_{gap}$, $\Delta_{e}$, $\Delta_{h}$ and $E_{W}$ as well as the $\Gamma$-point gap $\Delta_{\Gamma}$ for the entire ensemble of investigated TDBG models with twist angles down to 0.8∘. The following picture of the twist-angle dependence emerges. Firstly, as $\theta^{\star}$=1.3∘ is approached from above ($\theta>\theta^{\star}$), the ISP-induced energy gap $\Delta_{gap}$, which is of order 30 meV outside of the magic angle region (Fig. 3c), rapidly collapses to zero when the flattened bands are pre-formed and reopens only at $\theta<1^{\circ}$. Secondly, the bandwidth of the flat-band manifold $E_{W}$ is dramatically suppressed at the magic angle $\theta^{\star}=1.3^{\circ}$, where the bands become remarkably flat (Figs. 3b,f). Thirdly, within the same narrow region around $\theta^{\star}$, energy gaps $\Delta_{e}$ and $\Delta_{h}$ are maximized. The flat-band manifold is separated from the rest of bands within a twist angle region $1^{\circ}<\theta<2^{\circ}$ (Figs. 3d,e). Yet the magic angle 1.3*∘* is a well-defined quantity at which the band width is minimized, while the gap to excited states is close-to-the-maximum value of 37.6 meV (the maximum is offset towards 1.35*∘* degrees at which the $\Delta_{e}$=42.7 meV). Figures 3c-g also show that the relaxation effects are pivotal at the magic angle, not only flattening the entire band manifold, but also increasing by an order of magnitude the gap to excited states (Fig. 3b), making the TDBG half-bandwidth/bandgap ratio $14.7\%$ as favorable as in the parent TBG superlattice Carr et al. (2019).

Origin of band flattening.— We can address the origin of the band flattening at the magic angle $\theta^{\star}=1.3^{\circ}$ by constructing a minimalistic continuum model based on Refs. Tarnopolsky et al., 2019; San-Jose et al., 2012. At small twist angles, the distinction between commensurate and incommensurate moiré structures is experimentally not relevant, and we can treat the physically important moiré branch $L(\theta)=a/2\sin(\theta/2)$ as a continuum function. To construct the minimal TDBG Hamiltonian, we start from an effective model for the Bernal-stacked bilayer graphene Castro Neto et al. (2009)

[TABLE]

where $k=k_{x}+ik_{y}$. The terms $\gamma_{3}$ and $\gamma_{4}$ represent trigonal warping and particle-hole asymmetry. When two graphene bilayers are twisted to form TDBG, the inner graphene layers (2nd and 3rd) are coupled through the twist-induced interlayer coupling $T_{\theta}(\mathbf{r})$. To treat both $T_{0}$ and $T_{\theta}$ consistently, for now we neglect the $\gamma_{3}$ and $\gamma_{4}$ terms, which gives

[TABLE]

with a moiré three-fold star of $\mathbf{q}_{i}$, $|\mathbf{q}_{j}|=2k_{D}\sin\nicefrac{{\theta}}{{2}}$, equirotated by $\phi=2\pi/3$, and

[TABLE]

with $\hat{\Omega}_{\phi}=\cos\phi\,\sigma_{x}-\sin\phi\,\sigma_{y}$. Here, $\bm{\mathcal{G}}^{(0)}_{\theta}=0$, $\bm{\mathcal{G}}^{(1)}_{\theta}=\mathbf{q}_{2}-\mathbf{q}_{1}$, $\bm{\mathcal{G}}^{(2)}_{\theta}=\mathbf{q}_{3}-\mathbf{q}_{1}$ are the moiré reciprocal cell vectors and $\mathbf{d}$ is the relative displacement of one bilayer with respect to another one. The effective TDBG Hamiltonian thus reads

[TABLE]

where $\mathcal{H}_{1,2}=-i\bm{\sigma}_{-\theta/2}\bm{\nabla}+\Delta_{1,2}$ and $\mathcal{H}_{3,4}=-i\bm{\sigma}_{+\theta/2}\bm{\nabla}+\Delta_{3,4}$, where $\Delta_{i}$ is electric potential on the i${}^{\text{th}}$ layer.

To understand the origin of band flattening at the magic angle, we further construct the minimal model reproducing perfectly flat bands. For this, we for a moment neglect $\Delta_{i}$ and switch off $w_{AA}$, which is a natural consequence of lattice relaxation at small twist angles Tarnopolsky et al. (2019). To be self-consistent in $T_{0}$ and $T_{\theta}$, we set $\gamma_{1}=3w_{AB}$, thus condition $T_{\theta=0}(\mathbf{d}_{AB})=T_{0}$ is satisfied. In this limit, the minimal TDBG Hamiltonian (3) has particle-hole and chiral symmetries. We see that upon these imposed conditions the bands become absolutely flat at the magic angle (Fig. 4). Similar to the TBG case, this absolutely flat band becomes highly dispersive both below and above the magic angle, making the value of the magic angle rigorously defined. One can show that the TDBG Hamiltonian maps directly onto the TBG case in the chirally-symmetric limit. Indeed, Hamiltonian (3) can be rewritten as

[TABLE]

which reminds the TKV model for TBG Tarnopolsky et al. (2019). One can further show that this Hamiltonian maps onto two TBG Hamiltonians SI , where the band flatness comes mathematically from the flatness of the lowest Landau level in the quantum Hall effect on torus Tarnopolsky et al. (2019); Liu et al. (2019b). Therefore, the band flatness in the magic-angle TDBG is of the same topological nature. More realistic description of the flat bands in TDBGs at magic angle $\theta^{\star}$ would require including the interplay between trigonal warping, particle-hole asymmetries and the ISP fields as a perturbation around the perfectly flat bands, but the well-defined magic angle would be preserved.

To our surprise, the introduced minimalistic model is qualitatively consistent with our atomistic calculations (Fig. 4), including the reasonable values of band gaps both at the magic angle and in its vicinity. We stress that the effect of the ISP field is more important at higher angles, opening a gap between otherwise touching parabolic bands. On contrary, at the magic angle, the lattice relaxation effects and particle-hole asymmetries suppress the ISP effect (Fig. 3b), so we neglect this intrinsic polarization in the continuum model results shown in Fig. 4. The effects of lattice relaxation on the band structure of TDBG are more complex than in the case of TBG, and cannot be qualitatively accounted for in the continuum model by simple reduction of $w_{AA}$. For example, taking trial parameters $\gamma_{3}$, $\gamma_{4}$ expected from the BLG case and estimative $w_{AA}\approx 0.8\,w_{AB}$ does not reproduce the band gap and band width at the magic angle, when compared to the tight-binding values. We anticipate the values of $\gamma_{3}$ and $\gamma_{4}$ to be significantly renormalized by the lattice relaxation effects, which we will address in a further study.

Conclusions.—With the help of atomistic calculations we show that twisted double bilayer graphene has a well-defined magic angle, at which it hosts isolated flat bands, gapped out from higher excited states by dramatic 38 meV. In terms of the relative flatness (bandwidth/bandgap ratio), the TDBG at the magic angle 1.3*∘* is very close to the parent TBG heterostructure. A direct algebraic mapping to the TBG case can be found in the continuum setting if particle-hole asymmetries are neglected. A surprising novel detail revealed by the DFT calculations is the intrinsic symmetric polarization of the TDBG layers, which has not been reported previously. This internal effect modifies the magic-angle band dispersion. Previously, it was considered that the bands in TDBG can be made relatively flat only by applying external electric fields Liu et al. (2019a); Koshino (2019). We however report that the bands in TDBG are naturally flat at the magic angle 1.3*∘* due to significant intrinsic polarization and lattice relaxation effects. Further application of electric fields provides a control mechanism over the bandwidth and flat-band structure, important for the fine-tuning with respect to the electron-electron interaction scale Liu et al. (2019a); Cao et al. (2019). Flat bands can promote exotic correlated states, such as unconventional superconductivity, fractional quantum Hall effect and ferromagnetism in the flat bands. These anticipated phenomena in TDBG remain to be addressed.

Acknowledgements. This work was supported by NCCR Marvel and the Swiss National Science Foundation (grant No. P2ELP2_175278). We thank Fernando Gargiulo for assistance.