Band structure and topological property of twisted double bilayer graphenes

TL;DR

This paper investigates the electronic band structure and topological properties of twisted double bilayer graphene, revealing how different stacking arrangements influence their insulating and topological behaviors under various conditions.

Contribution

It provides a comparative analysis of AB-AB and AB-BA twisted double bilayer graphene, highlighting their distinct topological properties and responses to electric fields using continuum Hamiltonian calculations.

Findings

AB-AB is a trivial insulator without electric field.

AB-BA exhibits valley Hall insulator behavior.

Electric field opens a gap and causes electron-hole asymmetry.

Abstract

We study the electronic band structure and the topological property of the twisted double bilayer graphene, or a pair of AB-stacked bilayer graphenes rotationally stacked on top of each other. We consider two different arrangements, AB-AB and AB-BA, which differ in the relative orientation. For each system, we calculate the energy band and the valley Chern number using the continuum Hamiltonian. We show that the AB-AB and the AB-BA have similar band structures, while the Chern numbers associated with the corresponding bands are completely different. In the absence of the perpendicular electric field, in particular, the AB-AB system is a trivial insulator when the Fermi energy is in a gap, while the AB-BA is a valley Hall insulator. Also, the lowest electron and hole bands of the AB-AB are entangled by the symmetry protected band touching points, while they are separated in the AB-BA. In both cases, the perpendicular electric field immediately opens an energy gap at the charge neutral point, where the electron branch becomes much narrower than the hole branch, due to the significant electron-hole asymmetry.

Equations

$$H_{\text{AB-AB}}=\left(\begin{array}{cccc}{H}_0({\mathbf{k}}_1)& g^{†}({\mathbf{k}}_1)& & \\ g({\mathbf{k}}_1)& H_0^{\prime }({\mathbf{k}}_1)& U^{†}& \\ & U& H_0({\mathbf{k}}_2)& g^{†}({\mathbf{k}}_2)\\ & & g({\mathbf{k}}_2)& H_0^{\prime }({\mathbf{k}}_2)\end{array}\right)+V,$$

$$H_0(\mathbf{k})=\left(\begin{array}{cc}0& -\hslash v{k}_{-}\\ -\hslash v{k}_{+}& {\Delta }^{\prime }\end{array}\right),H_0^{\prime }(\mathbf{k})=\left(\begin{array}{cc}{\Delta }^{\prime }& -\hslash v{k}_{-}\\ -\hslash v{k}_{+}& 0\end{array}\right)$$

$$g(\mathbf{k})=\left(\begin{array}{cc}\hslash v_4{k}_{+}& {\gamma }_1\\ \hslash v_3{k}_{-}& \hslash v_4{k}_{+}\end{array}\right),$$

$$U$$

$$+\left(\begin{array}{cc}u& u^{\prime }{\omega }^{\xi }\\ u^{\prime }{\omega }^{-\xi }& u\end{array}\right)e^{i\xi ({\mathbf{G}}_1^{\mathrm{M}}+{\mathbf{G}}_2^{\mathrm{M}})\cdot \mathbf{r}},$$

$$V=\left(\begin{array}{cccc}\frac{3}{2}\Delta \hat{1}& & & \\ & \frac{1}{2}\Delta \hat{1}& & \\ & & -\frac{1}{2}\Delta \hat{1}& \\ & & & -\frac{3}{2}\Delta \hat{1}\end{array}\right),$$

$$H_{\text{AB-BA}}=\left(\begin{array}{cccc}{H}_0({\mathbf{k}}_1)& g^{†}({\mathbf{k}}_1)& & \\ g({\mathbf{k}}_1)& H_0^{\prime }({\mathbf{k}}_1)& U^{†}& \\ & U& H_0^{\prime }({\mathbf{k}}_2)& g({\mathbf{k}}_2)\\ & & g^{†}({\mathbf{k}}_2)& H_0({\mathbf{k}}_2)\end{array}\right)+V.$$

$$C_n=\frac{1}{2\pi }{\int }_{\mathrm{MBZ}}{\mathcal{F}}_{n,\mathbf{k}}{d}^{2}k,$$

$${\mathcal{F}}_{n,\mathbf{k}}=\frac{\partial a_{n,\mathbf{k}}^{(y)}}{\partial k_x}-\frac{\partial a_{n,\mathbf{k}}^{(x)}}{\partial k_y},a_{n,\mathbf{k}}^{(\mu )}=\frac{1}{i}⟨u_{n,\mathbf{k}}|\frac{\partial }{\partial k_{\mu }}|u_{n,\mathbf{k}}⟩,$$

$${\Sigma }^{-1}{H}_{\text{AB-AB}}\Sigma =-H_{\text{AB-AB}}^{\ast },$$

$$\Sigma =\left(\begin{array}{cccc}& & & {\sigma }_x\\ & & -{\sigma }_x& \\ & {\sigma }_x& & \\ -{\sigma }_x& & & \end{array}\right).$$

$$({\Sigma }^{\prime -1}\tilde{P})H_{\text{AB-BA}}(\tilde{P}{\Sigma }^{\prime })=-H_{\text{AB-BA}},$$

$${\Sigma }^{\prime }=\left(\begin{array}{cccc}& & & \hat{1}\\ & & -\hat{1}& \\ & \hat{1}& & \\ -\hat{1}& & & \end{array}\right),$$

Full text

I Introduction

The electronic property in a stack of two-dimensional (2D) materials sensitively depends on the relative twist angle $\theta$ between the adjacent layers, and we often have dramatic angle-dependent phenomena which are never observed in an isolated layer. The best known example is the twisted bilayer graphene (twisted BLG), or a rotationally stacked pair of monolayer graphenes, where a long-period moiré interference pattern significantly modifies the Dirac dispersion Lopes dos Santos et al. (2007); Mele (2010); Trambly de Laissardière et al. (2010); Shallcross et al. (2010); Morell et al. (2010); Bistritzer and MacDonald (2011); Moon and Koshino (2012); de Laissardiere et al. (2012); Moon and Koshino (2013); Weckbecker et al. (2016). Recently, the superconductivity and correlated insulating states are discovered in the magic-angle twisted BLG with extremely flat bands, Cao et al. (2018a, b); Yankowitz et al. (2019) and it is followed by a number of theoretical studies on the detailed properties of the flat bands and the possible mechanism of the superconductivity. Yuan and Fu (2018); Po et al. (2018); Xu and Balents (2018); Kang and Vafek (2018); Koshino et al. (2018); Ochi et al. (2018); Isobe et al. (2018); Dodaro et al. (2018); Padhi et al. (2018); Wu et al. (2018); Tarnopolsky et al. (2019); Zou et al. (2018) Graphene on hexagonal boron nitride (hBN) also exhibits the moiré-induced physics such as the formation of the secondary Dirac bands and the miniband structure. Kindermann et al. (2012); Wallbank et al. (2013); Mucha-Kruczyński et al. (2013); Jung et al. (2014); Moon and Koshino (2014); Dean et al. (2013); Ponomarenko et al. (2013); Hunt et al. (2013); Yu et al. (2014) A recent experiment reported the correlated insulating states in ABC-trilayer graphene on hBN, which is tunable by the external gate electric field. Chen et al. (2019) Controlling the twist angle in a stack of 2D materials provides powerful means to manipulate quantum properties of the electronic systems.

In this paper, we study a different type of moiré system, the twisted double bilayer graphene, which is composed of a pair of AB-stacked BLGs rotationally stacked on top of each other. The AB-stacked BLG is the most stable form of bilayer graphene which has the stacking structure of graphite. McCann and Koshino (2013) Here we consider two different arrangements, AB-AB and AB-BA, as illustrated in Figs. 1(a) and 1(b), respectively, where the AB-BA is obtained just by 180∘ rotation of the second BLG in the AB-AB. For each case, we derive the continuum Hamiltonian by extending the approach for the twisted BLGLopes dos Santos et al. (2007); Bistritzer and MacDonald (2011); Kindermann and First (2011); Lopes dos Santos et al. (2012); Moon and Koshino (2013); Koshino (2015); Koshino and Moon (2015); Weckbecker et al. (2016), and calculate the energy bands as well as the valley Chern numbers. Here we include the interlayer asymmetric potential $\Delta$ induced by the gate electric field.

The energy band structures are found to be similar between the AB-AB and the AB-BA, but the topological nature is different. In the absence of $\Delta$, the lowest electron and hole bands of the AB-AB are entangled by the symmetry protected band touching points, while they are separated in the AB-BA due to the different space symmetry. In both cases, the asymmetric potential $\Delta$ immediately opens an energy gap at the charge neutral point. We find that the graphite band parameters such as $\gamma_{3}$ and $\gamma_{4}$ play an important role in the electron-hole asymmetry, where the electron band becomes much narrower than the hole band as increasing $\Delta$.

The crucial difference between AB-AB and AB-BA is found in the Chern number. In the absence of $\Delta$, in particular, the AB-AB double bilayer becomes a trivial insulator because the symmetry requires all the Chern numbers to vanish, while the AB-BA is a valley Hall insulator with finite Chern number. We demonstrate the evolution of the Chern numbers as a function of $\Delta$, where we see that the energy bands of AB-AB and AB-BA carry completely different topological numbers, even though the band structures are similar. The difference in the Chern number would be observed by the measurement of the valley Hall conductivity Mak et al. (2014); Shimazaki et al. (2015), and also by the Landau level structure in the magnetic field.

This paper is organized as follows: In Sec. II, we define the lattice structures of AB-AB and AB-BA double bilayers, and then introduce a continuum Hamiltonian for each system in Sec. III. In Sec. IV, we study the band structures and the evolution of Chern numbers as a function of the twist angle and the asymmetric potential, where we discuss in detail about similarity and difference between the two systems. A brief conclusion is presented in Sec. V.

II Atomic structure

The AB-stacked BLG is composed of a pair of monolayer graphenes, with four atoms in the unit cell, labeled $A_{1}$, $B_{1}$ on the layer 1 (upper layer) and $A_{2}$, $B_{2}$ on the layer 2 (lower layer). McCann and Koshino (2013) The two graphene layers are arranged so that $B_{1}$ and $A_{2}$ are vertically located. We refer to these two atomic sites as dimer sites because the electronic orbitals on them are strongly coupled. The other two atoms, $A_{1}$ and $B_{2}$ are directly above or below the hexagon center of the other layer, and are referred to as non-dimer sites.

We compose the twisted AB-AB double bilayer graphene by stacking the first AB-stacked BLG (layers 1 and 2) on top of the second AB-stacked BLG (layers 3 and 4), as Fig. 1(a). We start from the non-rotated geometry where $B_{1}$, $A_{2}$, $B_{3}$ and $A_{4}$ are vertically aligned at the origin $O$, and then rotate the first and the second BLGs around $O$ by $-\theta/2$ and $+\theta/2$, respectively. The system has a three-fold in-plane rotation $C_{3z}$ symmetry along the $z$-axis (perpendicular to the layer), and a two-fold rotation $C_{2x}$ along the $x$-axis. The twisted AB-BA double bilayer can be defined just by rotating the second BLG (layer 3 and 4) of the AB-AB by 180∘ as in Fig. 1(b), where we flip the definition of $A$ site and $B$ site for layer 3 and 4. The system is symmetric under a three-fold in-plane rotation $C_{3z}$ and a two-fold rotation $C_{2y}$ along the $y$-axis.

We define $\mathbf{a}_{1}=a(1,0)$ and $\mathbf{a}_{2}=a(1/2,\sqrt{3}/2)$ as the lattice vectors of the initial BLGs before the rotation, where $a\approx 0.246\,\mathrm{nm}$ is the lattice constant of graphene. The corresponding reciprocal lattice vectors are $\mathbf{a}^{*}_{1}=(2\pi/a)(1,-1/\sqrt{3})$ and $\mathbf{a}^{*}_{2}=(2\pi/a)(0,2/\sqrt{3})$. After the rotation, the lattice vector of the $l$-th BLG is given by $\mathbf{a}_{i}^{(l)}=R(\mp\theta/2)\mathbf{a}_{i}$ with $\mp$ for $l=1,2$, respectively, where $R(\theta)$ represents the rotation matrix by $\theta$. The reciprocal lattice vectors become $\mathbf{a}_{i}^{*(l)}=R(\mp\theta/2)\mathbf{a}^{*}_{i}$. In a small $\theta$, the reciprocal lattice vectors for the moiré pattern is given by $\mathbf{G}^{\rm M}_{i}=\textbf{a}^{*(1)}_{i}-\textbf{a}^{*(2)}_{i}\,(i=1,2)$, and the real-space lattice vectors $\mathbf{L}^{\rm M}_{j}$ can then be obtained from $\mathbf{G}^{\rm M}_{i}\cdot\mathbf{L}^{\rm M}_{j}=2\pi\delta_{ij}$. A moiré unit cell is spanned by $\mathbf{L}^{\rm M}_{1}$ and $\mathbf{L}^{\rm M}_{2}$. The lattice constant $L_{\rm M}=|\mathbf{L}^{\rm M}_{1}|=|\mathbf{L}^{\rm M}_{2}|$ is $L_{\rm M}=a/[2\sin(\theta/2)]$. Figure 1(c) illustrates the Brillouin zone folding, where two large hexagons represent the first Brillouin zones of the first and the second BLGs, and the small hexagon is the moiré Brillouin zone of the twisted double BLG. The graphene’s Dirac points (the band touching points) are located at $\mathbf{K}^{(l)}_{\xi}=-\xi[2\mathbf{a}^{(l)*}_{1}+\mathbf{a}^{(l)*}_{2}]/3$ for the $l$-th BLG, where $\xi=\pm 1$ is the valley index. We label the symmetric points of the moiré Brillouin zone as $\bar{\Gamma}$, $\bar{M}$, $\bar{K}$ and $\bar{K^{\prime}}$ as in Fig. 1(c).

III Continuum Hamiltonian

To describe the electronic band structure of the twisted double bilayers, we adopt the continuum method based on the Dirac Hamiltonian.Lopes dos Santos et al. (2007); Bistritzer and MacDonald (2011); Kindermann and First (2011); Lopes dos Santos et al. (2012); Moon and Koshino (2013); Koshino (2015); Koshino and Moon (2015); Weckbecker et al. (2016) The validity of the continuum model was verified for twisted BLG by the direct comparison to the tight-binding model. Moon and Koshino (2013); Weckbecker et al. (2016) We define the Bloch bases of $p_{z}$ orbitals at sublattice $X=A_{1},B_{1},\cdots,A_{4},B_{4}$ as $|\mathbf{k},X\rangle=N^{-1/2}\sum_{\mathbf{R}_{X}}e^{i\mathbf{k}\cdot\mathbf{R}_{X}}$, where $|\mathbf{R}_{X}\rangle$ is the atomic $p_{z}$ orbital at the site $\mathbf{R}_{X}$, $\mathbf{k}$ is the two-dimensional Bloch wave vector and $N$ is the number of same sublattices in the system. The continuum Hamiltonian for twisted AB-AB double bilayer graphene at small twist angle $\theta(\ll 1)$ is written in 8 $\times$ 8 matrix for the Bloch bases of $(A_{1},B_{1},A_{2},B_{2},A_{3},B_{3},A_{4},B_{4})$ as

[TABLE]

where $\mathbf{k}_{l}=R(\pm\theta/2)({\mathbf{k}}-\mathbf{K}^{(l)}_{\xi})$ with $\pm$ for $l=1$ and 2, respectively, and

[TABLE]

with $k_{\pm}=\xi k_{x}\pm ik_{y}$. $H_{0}$ and $H^{\prime}_{0}$ are the Hamiltonian of monolayer graphene where $\Delta^{\prime}=0.050$eVMcCann and Koshino (2013) represents the on-site potential of dimer sites with respect to non-dimer sites. The parameter $v$ is the band velocity of monolayer graphene, and it is taken as $\hbar v/a=2.1354$ eV.Moon and Koshino (2013); Koshino et al. (2018) The matrix $g$ describes the interlayer coupling of the AB-stacked BLG, where $\gamma_{1}=0.4$eV is the coupling between dimer sites, and $v_{3}$ and $v_{4}$ are related to diagonal hoppings $\gamma_{3}=0.32$eV and $\gamma_{4}=0.044$ eV with the relation $v_{i}=(\sqrt{3}/2)\gamma_{i}a/\hbar\,(i=3,4)$. McCann and Koshino (2013) In the AB-stacked BLG, $v_{3}$ is responsible for the trigonal warping of the energy band and $v_{4}$ is for the electron-hole asymmetry.

The matrix $U$ is the moiré interlayer coupling between twisted layers given by Bistritzer and MacDonald (2011); Moon and Koshino (2013); Koshino et al. (2018)

[TABLE]

where $\omega=e^{2\pi i/3}$, $u=0.0797$eV and $u^{\prime}=0.0975$eV Koshino et al. (2018) are the amplitudes of diagonal and off-diagonal terms, respectively. The difference between $u$ and $u^{\prime}$ effectively describe the out-of-plane corrugation effect, which enhances the energy gaps between the lowest energy bands and the excited bands.Koshino et al. (2018); Nam and Koshino (2017); Tarnopolsky et al. (2019) Lastly, $V$ is the interlayer asymmetric potential,

[TABLE]

where $\hat{1}$ is $2\times 2$ unit matrix, and $\Delta$ represents the difference in the electrostatic energy between the adjacent layers. Here we simply assumed the perpendicular electric field is constant.

Noting that the lattice structure of the AB-AB double bilayer has $C_{2x}$ symmetry and also the valley degree of freedom $\xi=\pm$ is unchanged under $C_{2x}$, the Hamiltonian $H_{\textrm{AB-AB}}$ of each single valley commutes with $C_{2x}$, given that the asymmetric potential $\Delta$ is absent.

Similarly, the Hamiltonian of the twisted AB-BA double bilayer graphene is given by

[TABLE]

where $H_{0}(\mathbf{k}_{2})$ and $H^{\prime}_{0}(\mathbf{k}_{2})$ are interchanged and also $g(\mathbf{k}_{2})$ and $g^{\dagger}(\mathbf{k}_{2})$ are swapped in ${H}_{\textrm{AB-AB}}$. The lattice structure of the AB-BA double bilayer has $C_{2y}$ symmetry, and $C_{2y}$ interchanges the valleys $\xi=\pm$. As a result, the Hamiltonian $H_{\textrm{AB-BA}}$ with $\Delta=0$ commutes with $C_{2y}T$, where $T$ is the time reversal operator.

The calculation of the energy bands and the eigenstates is performed in the $k$-space picture. For a single Bloch vector $\mathbf{k}$ in the moiré Brillouin zone, the interlayer coupling $U$ hybridizes the graphene’s eigenstates at $\mathbf{q}=\mathbf{k}+\mathbf{G}$, where $\mathbf{G}=m_{1}\mathbf{G}^{\rm M}_{1}+m_{2}\mathbf{G}^{\rm M}_{2}$ and $m_{1}$ and $m_{2}$ are integers. The low-energy eigenstates can be obtained by numerically diagonalizing the Hamiltonian within the limited number of $\mathbf{q}$’s inside the cut-off circle $|\mathbf{q}-\mathbf{q}_{0}|<q_{c}$. Here $\mathbf{q}_{0}$ is taken as the midpoint between $\mathbf{K}^{(1)}_{\xi}$ and $\mathbf{K}^{(2)}_{\xi}$, and $q_{c}$ is set to $4|\mathbf{G}^{\rm M}_{1}|$. The calculation is done independently for each of $\xi=\pm$ as the intervalley coupling can be neglected in small twist angles.

We calculate the Chern number of moiré subbands by the standard definition,

[TABLE]

where $n$ is the band index, MBZ represents the moiré Brillouin zone, and $F_{n,\mathbf{k}}$ is the Berry curvature defined by

[TABLE]

where $u_{n,\mathbf{k}}$ is the Bloch wave function of $n$-th subband. We numerically calculate the Chern numbers using the discretizing method. Fukui et al. (2005)

The symmetry imposes constraints on the Chern number. For the AB-AB double bilayer at $\Delta=0$, the $C_{2x}$ symmetry requires ${\cal F}_{n,(k_{x},-k_{y})}=-{\cal F}_{n,(k_{x},k_{y})}$, so that the Chern number of each single band must vanish. In the AB-BA double bilayer at $\Delta=0$, the $C_{2y}T$ symmetry requires ${\cal F}_{n,(k_{x},-k_{y})}={\cal F}_{n,(k_{x},k_{y})}$, and the Chern number can be finite.

IV Band structures and topological properties

IV.1 Minimal model

Before calculating the band structure with all the band parameters fully included, it is intuitive to consider the minimal model which neglects the relatively small parameters, $v_{3}$, $v_{4}$, $\Delta^{\prime}$ and the rotation matrix $R(\pm\theta/2)$ in the definitions of $\mathbf{k}^{(1)}$ and $\mathbf{k}^{(2)}$. Then the AB-AB Hamiltonian Eq. (1) has a fictitious particle-hole symmetry similar to TBGMoon and Koshino (2013),

[TABLE]

This immediately leads to the electron-hole symmetry in the energy bands, $E_{n,\mathbf{k}}=-E_{-n,-\mathbf{k}}$, and also the anti-symmetric relation in the Chern number, $C_{-n}=-C_{n}$, where $n$ and $-n$ stand for the band indexes of the corresponding electron and hole bands, respectively. Note that Eq. (9) holds even in the presence of the interlayer asymmetric potential $\Delta$.

The AB-BA Hamiltonian Eq. (6) has a different type of symmetry between the electron and the hole bands,

[TABLE]

where $\tilde{P}$ is a space inversion operator which works on envelop function as $\tilde{P}F_{X}(\mathbf{r})=F_{X}(-\mathbf{r})$, while it does not change the sublattice degree of freedom ($X=A_{1},B_{1},\cdots$). This again forces the electron-hole symmetry $E_{n,\mathbf{k}}=-E_{-n,-\mathbf{k}}$, but the Chern number becomes electron-hole symmetric, $C_{-n}=C_{n}$, because the operation lacks the complex conjugate.

Figures 2(a) and (b) show the minimal-model band structure, calculated for the AB-AB double bilayer and the AB-BA double bilayer, respectively, at the twist angle $\theta=1.33^{\circ}$ with $\Delta=0$ and 20 meV. The band structures of the two systems closely resemble each other. At $\Delta=0$, we have a pair of energy bands touching at Dirac point, which are isolated from the excited bands by energy gaps as in the twisted BLG. Koshino et al. (2018) A finite $\Delta$ immediately opens an energy gap at the charge neutral point. This is in a sharp contrast to the twisted BLG, where the perpendicular electric field never opens a gap at the charge neutral point, because the band touching is protected by $C_{2}T$ symmetry. Now the twisted double bilayer lacks $C_{2}$ symmetry.

Although the band structures are pretty much similar between the AB-AB and the AB-BA cases, the properties of the Chern number are completely different. In Fig. 2, the black numbers indicate the Chern numbers of the central two bands in $\xi=+$ valley, and the blue numbers between the bands are the integrated Chern numbers summed over all the energy bands of $\xi=+$ below. Because of the time reversal symmetry, the Chern number of $\xi=-$ valley is opposite in sign to $\xi=+$. We actually see the expected relation $C_{-n}=-C_{n}$ for the AB-AB, and $C_{-n}=C_{n}$ for the AB-BA. In the absence of the asymmetric potential $\Delta$, the Chern numbers all vanish in the AB-AB because of the rigorous symmetry $C_{2x}$ mentioned in the previous section, while it is finite in the AB-BA. When the Fermi energy is inside one of those gaps, therefore, the AB-BA double bilayer is a valley Hall insulator, while AB-AB is a trivial insulator. The Chern numbers can be finite in the AB-AB once the asymmetric potential $\Delta$ is switched on, because it breaks $C_{2x}$. The integrated Chern numbers inside the central gap is $-2$ in the AB-AB, while 0 in the AB-BA. This is just equal to the sum of Chern numbers of two independent gapped BLGs, which is $-1$ for the AB stack while $+1$ for the BA stack. Martin et al. (2008); Koshino (2008); Zhang et al. (2013); Vaezi et al. (2013)

IV.2 Full parameter model

Inclusion of the additional band parameters neglected in the minimal model causes a significant change particularly in the low-energy sector. Figures 3(a) and (b) show the full-parameter band structure of the AB-AB and the AB-BA double bilayers, respectively, at the twist angle $\theta=1.33^{\circ}$ with $\Delta=0,5$ and 20 meV. The thin green lines indicate the energy bands of the minimal model [Fig. 2] for quantitative comparison. We see that the energy bands are now electron-hole asymmetric because the fictitious symmetry of Eq. (9) or Eq. (10) is broken. The band structures of AB-AB and AB-BA are still similar, but there are several important differences. At $\Delta=0$, in particular, the central energy bands of the AB-AB are touching at two points on the $\bar{\Gamma}-\bar{M}$ line, while they are anti-crossing in the AB-BA. The band touching of the former is protected by the $C_{2x}$ symmetry. Since the $k$-points on $\bar{\Gamma}-\bar{M}$ are invariant under $C_{2x}$ operation, the Bloch states on the line can be characterized by the eigenvalues of $C_{2x}$. The energy bands crossing at the center actually have the opposite eignevalues $C_{2x}=\pm 1$, so that they are never hybridized. The energy bands form a two-dimensional Dirac cone around each band touching point. Because of $C_{3}$ symmetry, we have six touching points in each single valley. Note that the energy bands of $\xi=-$ is just 180∘ rotation of $\xi=+$ band, so the band touching of $\xi=-$ are not seen in the figure.

In increasing $\Delta$, we see that the upper central band (the first conduction band) becomes much narrower than the lower central band (the first valence band), in both of the AB-AB and the AB-BA. As a result, the energy gap just above the upper band survives in relatively large $\Delta$, while the gap below the lower band is easily masked by the wide dispersion.

The properties of the Chern number are mostly carried over from the minimal model. A difference is seen in $\Delta=5$ meV in the AB-AB case, where the central bands have the Chern number $\pm 3$, unlike $\pm 2$ in the minimal model. This is attributed to the six Dirac points at $\Delta=0$, each of which contributes to the Berry curvature $\pi$ when gapped out. In even increasing $\Delta$, we have a band touching at $\bar{K}$ around $\Delta\sim 9$ meV, where the Chern number $+1$ is transferred from the lower central band to the higher central band. As a result, the Chern number of the central bands becomes $\pm 2$ as in the minimal model. A similar topological change is also observed in the AB-BA double bilayer, where the Chern numbers of the central two bands change from $(2,0)$ to $(1,1)$.

Finally, we present in Fig. 4 the twist angle dependence of the band structure in (a) the AB-AB and (b) the AB-BA double bilayers with $\Delta=0$. For the AB-AB, the central electron and hole bands are separated by an energy gap at $\theta=2^{\circ}$, and they get closer in decreasing $\theta$. At $\theta\sim 1.44^{\circ}$, there is a quadratic band touching on $\bar{\Gamma}-\bar{M}$ line, and a pair of the Dirac points are formed below that angle. Those band touching points are protected by $C_{2x}$ symmetry as already argued. In even lower angles, the central energy bands become narrower and narrower, and at the same time the higher excited bands collides with the central bands. At $\theta=1.1^{\circ}$ and $0.9^{\circ}$, we have an insulating gap between the third and the fourth valence bands. Because of $C_{2x}$ symmetry, the Chern number is zero everywhere as long as $\Delta=0$.

In the AB-BA case, we have a similar evolution of the band structure, while the Chern number is generally non-zero. At $\theta=2^{\circ}$, the charge neutral point is a valley Hall insulator with the Chern number 2. In decreasing $\theta$, we have a topological change at $\theta\sim 1.45^{\circ}$, where the Chern number 3 is transferred from the lower band to the higher band through the three touching points arranged in 120∘symmetry. Unlike the AB-AB, the band touching occurs only at the topological transition, and the bands are separated again after the transition. In smaller angles less than 1∘, the central bands touch with the excited bands, and there are complex topological changes between them. We have a new insulating gap between the third and the fourth valence bands, where the Chern number is $-2$.

In this work, we assumed different parameters $u$ and $u^{\prime}$ to describe possible corrugation effect, where the adopted values, $u=0.0797$eV and $u^{\prime}=0.0975$eV, are taken from the twisted BLG (monolayer-monolayer) Koshino et al. (2018). In the twisted double BLG, however, the corrugation would be reduced to some extent considering that bilayer graphene is stiffer than monolayer graphene, and then the difference between $u$ and $u^{\prime}$ should also decrease accordingly. To see this effect, we present the band structure assuming $u=u^{\prime}=0.0975$eV in Fig. 5, where thin blue lines represent the original calculation from Fig. 3. We see that the qualitative features are similar, while the energy gaps between the lowest bands and the excited bands are smaller in the $u=u^{\prime}$ model than in the $u\neq u^{\prime}$ model, as in the twisted BLG. Koshino et al. (2018); Nam and Koshino (2017); Tarnopolsky et al. (2019). The real situation should be somewhere between the two cases.

V Conclusion

We have studied the electronic band structure and the Chern numbers of AB-AB and AB-BA twisted double bilayer graphenes, and found that the two systems have similar band structures, but with completely different topological properties. In the absence of the asymmetric potential $\Delta$ (perpendicular electric field), in particular, the AB-BA double bilayer is a valley Hall insulator when the Fermi energy is in a gap, while the AB-AB is a trivial insulator due to the symmetry constraint. Also, the energy bands of the AB-AB in $\Delta=0$ are entangled by the symmetry protected band touching points, while they are all separated in the AB-BA. The common features shared by the two systems is that a pair of narrow bands at the charge neutral point are immediately gapped by applying the perpendicular electric field, unlike the twisted BLG (monolayer-monolayer). There the graphite band parameters such as $\gamma_{3}$, $\gamma_{4}$ play an important role in the electron-hole asymmetry, where the electron branch becomes much narrower than the hole branch in increasing the perpendicular electric field.

Acknowledgments

MK thanks the fruitful discussions with Jeil Jung, Pablo Jarillo-Herrero, Philip Kim, Eslam Khalaf, Jong Yeon Lee and Ashvin Vishwanath. MK acknowledges the financial support of JSPS KAKENHI Grant Number JP17K05496.

Note added: The recent preprints reported experimental observations of superconductivity and correlated insulating states in the twisted double bilayer graphenes. Shen et al. (2019); Liu et al. (2019); Cao et al. (2019) During the completion of this work, we became aware of recent theoretical works on the electronic properties of twisted AB-AB double BLG Chebrolu et al. (2019); Choi and Choi (2019); Lee et al. (2019). Just after the submission of this manuscript, we have come to notice a recent theoretical study on the electronic and topological properties on twisted multilayer graphene systems with various stacking configurations. Liu and Dai (2019)