Higher-Order Topological Insulator in Twisted Bilayer Graphene

TL;DR

This paper proposes that twisted bilayer graphene at large angles can act as a higher-order topological insulator, exhibiting robust corner states due to bulk topology and symmetry considerations.

Contribution

It introduces twisted bilayer graphene as a new material platform for higher-order topological insulators with angle-independent corner states.

Findings

Corner states occur at half filling in twisted bilayer graphene.

Bulk gap opens due to intervalley scattering at large angles.

Corner states are robust as long as symmetries are preserved.

Abstract

Higher-order topological insulators are newly proposed topological phases of matter, whose bulk topology manifests as localized modes at two- or higher-dimensional lower boundaries. In this work, we propose the twisted bilayer graphenes with large angles as higher-order topological insulators, hosting topological corner charges. At large commensurate angles, the intervalley scattering opens up the bulk gap and the corner states occur at half filling. Based on both first-principles calculations and analytic analysis, we show the striking results that the emergence of the corner states do not depend on the choice of the specific angles as long as the underlying symmetries are intact. Our results show that the twisted bilayer graphene can serve as a robust candidate material of two-dimensional higher-order topological insulator.

Equations

$$\nu ={\nu }_{+}={\nu }_{-}(\mathrm{mod}2),$$

$${\nu }_{\pm }=\frac{1}{i\pi }\log \det [{\mathcal{U}}_{\pm }(\Gamma ,M_3){\mathcal{U}}_{\pm }(M_3,\Gamma )].$$

$$(-1{)}^{{\omega }_2}=\prod _{{\Gamma }_i\in \mathrm{TRIM}}(-1{)}^{[N_{\mathrm{occ}}^{-}({\Gamma }_i)/2]},$$

$$\det [{\mathcal{U}}_{\pm }(\Gamma ,M_3){\mathcal{U}}_{\pm }(M_3,\Gamma )]$$

$$=\det [{\mathcal{U}}_{\pm }(\Gamma ,M_3)C_{2z}(C_{2z}^{-1}{\mathcal{U}}_{\pm }(M_3,\Gamma )C_{2z})C_{2z}^{-1}]$$

$$=\prod _{n\in \mathrm{occ},C_{2x}=\pm 1}{\zeta }_n(\Gamma ){\zeta }_n(M_3),$$

$$(-1{)}^{{\nu }_{+}}$$

$$Q^{(6)}=\frac{e}{2}{w}_2+\frac{e}{6}(\#C_3)(\text{mod}e),$$

$${\mathbf{X}}_{m,n}\equiv \mathbf{X}+m{\mathbf{b}}_{M,1}+n{\mathbf{b}}_{M,2},$$

$$\displaystyle\hat{C}_{2z}=\left(\begin{array}[]{cc}\zeta({\bm{X}})&0\\ 0&\sigma_{x}\otimes I_{N}\\ \end{array}\right),$$

$$(-1{)}^{{\omega }_2}=\zeta (M)\zeta (\Gamma ).$$

$${\nu }_{\pm }$$

$$=\det [⟨u_i(\Gamma )|{\mathcal{U}}_{\pm }(\Gamma ,M_3){\mathcal{U}}_{\pm }(M_3,{\Gamma }^{\prime })|u_j(\Gamma )⟩]$$

$$=\det [⟨u_i(\Gamma )|{\mathcal{U}}_{\pm }(\Gamma ,M_3){\hat{C}}_{2z}({\hat{C}}_{2z}{\mathcal{U}}_{\pm }(M_3,{\Gamma }^{\prime }){\hat{C}}_{2z}){\hat{C}}_{2z}|u_j(\Gamma )⟩]$$

$$=\det [\sum _{a,b,c\in occ.}⟨u_i(\Gamma )|{\mathcal{U}}_{\pm }(\Gamma ,M_3)|u_a(M_3)⟩⟨u_a(M_3)|{\hat{C}}_{2z}|u_b(M_3)⟩⟨u_b(M_3)|{\mathcal{U}}_{\pm }(M_3,\Gamma )|u_c(\Gamma )⟩⟨u_c(\Gamma )|{\hat{C}}_{2z}|u_j(\Gamma )⟩]$$

$$=\det [\sum _{a\in occ.}⟨u_i(\Gamma )|{\hat{C}}_{2z}|u_a(\Gamma )⟩⟨u_a(M_3)|{\hat{C}}_{2z}|u_j(M_3)⟩]$$

$$=\prod _{n\in \mathrm{occ},C_{2x}=\pm 1}{\zeta }_n(\Gamma ){\zeta }_n(M_3),$$

$${\nu }_{\pm }=\prod _{n\in \mathrm{occ},C_{2x}=\pm 1}{\zeta }_n(\Gamma ){\zeta }_n(M_1){\zeta }_n(M_2){\zeta }_n(M_3).$$

$${\nu }_{\pm }=\prod _{{\Gamma }_i\in \text{TRIM}}(-1{)}^{N_{occ}^{-}({\Gamma }_i)/2}\equiv (-1{)}^{{\omega }_2}.$$

Full text

Introduction - Generalizing the concept of topology in diverse systems has been one of the most important topics in condensed matter physicsQi and Zhang (2011); Hasan and Kane (2010). One recently proposed interesting path of extending the knowledge of the topological phases of matter is to consider higher-order generalizationBenalcazar et al. (2017a). That is to say, the band topology of $d$-dimensional insulator manifests as non-trivial $d-2$ or lower-dimensional boundary states, which is known as the higher-order topological insulator(HOTI)Benalcazar et al. (2017b); Langbehn et al. (2017); Song et al. (2017); Schindler et al. (2018a); Ezawa (2018a); Serra-Garcia et al. (2018); Peterson et al. (2018); Imhof et al. (2018); Khalaf (2018); Geier et al. (2018); Khalaf et al. (2018); Kunst et al. (2018); Yan et al. (2018); Ezawa (2018b); Xue et al. (2019); van Miert and Ortix (2018); Trifunovic and Brouwer (2019); Ezawa (2018c); Goren et al. (2018); Lin and Hughes (2018); Liu et al. (2019); Ezawa (2018d); Serra-Garcia et al. (2019); Araki et al. (2019); Franca et al. (2018); Kooi et al. (2018); Wheeler et al. (2018); Kang et al. (2018); Ono et al. (2019). There have been several theoretical proposals of the higher-order topological insulators in both two- and three- dimensional solid state systems. For example, in three dimensions, the higher-order topological insulators harbor one-dimensional gapless modes, which is referred to as the hinge modes. Examples include a bismuth crystalSchindler et al. (2018b), SnTe, surface-modified $\text{Bi}_{2}\text{TeI}$, BiSe, BiTeSchindler et al. (2018a), phosphoreneEzawa (2018b), van der Waals multilayeresKooi et al. (2018), and transition metal dichalcogenides $\text{XTe}_{2}$ Wang et al. (2018). In two dimensions, [math]-D localized modes appear at the corners of the two-dimensional materials, known as the corner states. The proposed examples are two-dimensional phosphoreneEzawa (2018b), monolayer graphdiyneLee et al. (2019); Sheng et al. (2019), and effective model of the twisted bilayer grapheneAhn et al. (2019a); Kindermann (2015).

Although there already exist a few candidate materials for the two-dimensional HOTIs, experimental signature that unambiguously identifies such phases remains elusive to date. To this end, it is desirable to identify higher-order topological materials that are readily available and highly controllable such as bismuth and other two-dimensional materials uncovered in Ref. Kim et al. . Here, we propose that a large angle twisted bilayer graphene is a generic higher-order topological insulator, characterized by the topologically protected corner states. Twisted bilayer graphene has clear advantages over other candidates for the experimental detection of the corner charge: (i) In our proposal, the HOTI is realized at half filling, thus corner state occurs without any fine-tuning of chemical potential, which is distinct from other casesAhn et al. (2019a); Song et al. (2018). (ii) In graphene, ultra clean transfer techniques are available that enables isolation from undesirable substrate effectsLee et al. (2017). (iii) The HOTI in graphene is very unique system realized by intervalley scattering with negligible spin-orbit coupling. Thus, one can avoid any ambiguity arising from the spin-orbit coupling effect.

In this work, for the first time, we propose that the twisted bilayer graphenes with large commensurate angles can be generically the higher order topological insulators. We first present the first-principles calculations of the twisted bilayer graphene at the commensurate angle, $\theta=21.78^{\circ}$, as a representative example. In this case, we find that the sizable gap ($\sim$9 meV) exists which is originated from the absence of the $U(1)$ valley symmetry, $U(1)_{v}$. We then show that twisted bilayer graphene has non-trivial higher-order band topology, characterized by the occurrence of topological corner states with the fractional electron charge $e/2$ in the mirror symmetric corners and $C_{6}$ symmetric corners. In addition, we generalize our discussion by showing that the non-trivial band topology is guaranteed to exist regardless of the specific commensurate twist angle or the microscopic details of the atomic structure as long as the underlying symmetries are preserved.

Lattice model and Symmetries - We begin our discussion by introducing the atomic structure and the associated crystalline symmetries of twisted bilayer graphene. We consider a specific set of twisted bilayer structures that belong to the hexagonal space group # 177 (point group $D_{6}$). The atomic configuration can be readily constructed by twisting AA-stacked bilayer graphene with respect to the collinear axis at the hexagonal center [See Figs. 1(a) and 1(b)]. The twist preserves both $C_{6z}$ and $C_{2x}$ about the out-of-plane $z$- and in-plane $x$-axes, respectively, which generate another important rotation $C_{2y}$ about the $y$-axis. The original translation symmetry is broken by the twist, but we can define the Moiré translational symmetry, depending on the twist angle. For any coprime integers $p$ and $q$, a twist by $\theta_{p,q}\!=\!\arccos\frac{3p^{2}+3pq+q^{2}/2}{3p^{2}+3pq+q^{2}}$ results in an enlarged Moiré unit cell with the lattice constant $L\!=\!a\sqrt{\frac{3p^{2}+3pq+q^{2}}{{\rm gcd}(q,3)}}$ Zou et al. (2018); Lopes dos Santos et al. (2007), where $a$ is the original lattice constant and gcd represents the greatest common divisor.

Electronic energy bands - In the limit where $\theta\!\lesssim\!1^{\circ}$ without the lattice distortions, the Moiré potential has long periodicity in real space, resulting in negligible interaction between valleys. In this limit, including the so-called magic angles where the Fermi velocity vanishesBistritzer and MacDonald (2011), each valley is decoupled and the $U(1)$ valley symmetry $U(1)_{v}$ is approximately preserved and it, together with $C_{2z}\mathcal{T}$ symmetry, provides topological protection of four Dirac points associated with the $\mathbb{Z}_{2}$-quantized Berry phases $\pi$ Ahn and Yang (2017); Song et al. (2018). Here, $\mathcal{T}$ represents time-reversal symmetry, where $\mathcal{T}^{2}=1$ without spin-orbit coupling.

However, in the large angle limit $\theta\rightarrow 30^{\circ}$, which is the concern of this paper, the $U(1)_{v}$ symmetry can be broken and the fourfold-degenerate Dirac points can split into two pairs of massive Dirac points, which here we claim to induce a HOTI phase. The intervalley coupling, and thus the gap opening between Dirac points, has a tendency to increase as the twist angle $\theta$ increases. Eventually, a sizable band gap of around 10 meV order opens at $\theta=21.78^{\circ}$ ($p=1$, $q=3$)Shallcross et al. (2008). Focusing on this particular angle first, below we demonstrate the HOTI phase in twisted bilayer graphene using first-principles calculations.

Figure 2(a) shows the density functional theory (DFT) electronic band structures of the twisted bilayer graphene with $\theta$ = 27.78∘. The results show that the band structure features a narrow gap of $\sim$ 9 meV, which is more or less similar to the previous result of $\sim$ 7 meV Shallcross et al. (2008). While a global direct band gap appears throughout the entire Moiré BZ between the conduction and valence bands, the minimum value occurs in the vicinity of $K$ (as well as $K^{\prime}$) as shown in the right panel of Fig. 2(a). Using all the occupied 56 bands, we explore the higher-order band topology by directly calculating a mirror-winding number and the second Stiefel-Whitney number.

Mirror-winding number - Whereas we find that the Zak phase along any closed loop is trivial, a finer topological classification can be made along the mirror invariant line utilizing the presence of the mirror symmetry. The BZ has a mirror-invariant line $\Gamma$-$M_{3}$-$\Gamma$ preserving $C_{2x}$ symmetry [See Fig. 2(b)]. Along this mirror invariant line, we can decompose the Hamiltonian into two distinct sub-sectors characterized by the mirror eigenvalues $\pm 1$. Then, one can define the $\mathbb{Z}_{2}$ mirror winding number Chiu et al. (2016, 2013) as the mirror-resolved Zak phase, $\nu_{\pm}$, calculated along $\Gamma$-$M_{3}$-$\Gamma$ line for $C_{2x}=\pm 1$ subsector. The $\mathbb{Z}_{2}$ mirror winding number $\nu$ can be evaluated by

[TABLE]

where

[TABLE]

Here, $\mathcal{U}_{\pm}({\bm{k}}_{1},{\bm{k}}_{2})\equiv P\exp[i\int_{{\bm{k}}_{1}}^{{\bm{k}}_{2}}\bm{\mathcal{A}_{\pm}}({\bm{k}})\cdot d{\bm{k}}]=\tilde{P}_{\pm}({\bm{k}}_{1})[\prod_{\bm{k}}\tilde{P}_{\pm}({\bm{k}})]\tilde{P}_{\pm}({\bm{k}}_{2})$, $\mathcal{A}_{\pm}$ is the non-Abelian Berry connection evaluated from the $\pm$-sector, respectively. $P$ indicates the path ordering, and $\tilde{P}_{\pm}$ is the projection operator to the occupied mirror $\pm$-subspace. Using all the occupied DFT bands, we calculate the eigenphase $\theta_{\pm,i}\equiv{\rm Im}\log(u_{\pm,i})$, where $u_{\pm,i}$ is the $i$-th eigenvalue of the Wilson loop matrix $\mathcal{U}_{\pm}(\Gamma,M_{3})\mathcal{U}_{\pm}(M_{3},\Gamma)$. The results, presented in Fig. 2(c), show that $\nu=\nu_{\pm}=+1$, confirming the non-trivial band topology.

Stiefel-Whitney number - The physical manifestation of the non-trivial winding number $\nu=1$ is the emergence of the topological corner states as we discuss later. Such corner states can be further established by the second Stiefel-Whitney number Ahn et al. (2018). The twisted bilayer graphene preserves the two-dimensional inversion symmetry $C_{2z}$, whose product with time-reversal symmetries $C_{2z}\mathcal{T}$ imposes that the Hamiltonian is real symmetric matrix, thus being characterized by the Stiefel-Whitney classes. While the first Stiefel-Whitney number $\omega_{1}$, which is equivalent to the Zak phase, is turned out to be trivial, the non-trivial second $\mathbb{Z}_{2}$ Stiefel-Whitney number $\omega_{2}$ is expected from our DFT calculations. As proposed by Ahn et al. Ahn et al. (2018, 2019b), we evaluate $\omega_{2}$ from the parity eigenvalues of occupied bands at time-reversal invariant momenta (TRIM) $\Gamma_{i}\in\{\Gamma,M_{1},M_{2},M_{3}\}$ using

[TABLE]

where $N_{\rm occ}^{-}(\Gamma_{i})$ is the number of occupied bands with an odd parity at $\Gamma_{i}$. Figure 2(c) shows the calculated parity eigenvalues at TRIM. We find that $\Gamma$ ($M_{1,2,3}$) has 24 (30) odd bands out of the 56 occupied bands. Therefore the parities at $(\Gamma,M_{1},M_{2},M_{3})=(+,-,-,-)$, leading to $\omega_{2}=1$.

The mirror winding number $\nu$ and the second Stiefel-Whitney number $\omega_{2}$ are, in general, irrelevant topological invariants that depend on different symmetries ($C_{2x}$ and $C_{2z}\mathcal{T}$, respectively). However, in the twisted bilayer graphene,where both symmetries are present, we can formally equate these numbers. We decompose the Wilson loop into the two pieces of the Wilson lines related by the inversion symmetry [See Fig. 4(a)]. In this case, one can show the exact cancellation between the two Wilson lines up to the parity eigenvalues of $C_{2z}$ as following,

[TABLE]

where $\zeta_{n}({\bm{k}})$ is the eigenvalue of $C_{2z}$ at ${\bm{k}}$ on $n$-th band. Since $C_{3z}$ guarantees the same parity structure of the bands at $M_{1}$, $M_{2}$, and thus $\zeta_{n}(M_{1})\zeta_{n}(M_{2})=1$ for any band index $n$. Therefore, we finally equate the mirror winding number $\nu=\nu_{\pm}$ to the second Stiefel-Whitney number $\omega_{2}$,

[TABLE]

In the last line we used the absence of Zak phase (the first Stiefel-Whitney number $\omega_{1}=1$) along $\Gamma$-$M_{3}$-$\Gamma$, which enforces the same parity between the mirror-even and mirror-odd sectors 111See supplementary material at http://XXXXXXX.XXXX.XXX for the detailed proof of Eq. (5). .

Topological corner states-

After presenting the non-trivial topology of the bulk using the first-principles calculations, we utilize the tight-binding modelMoon and Koshino (2013) to demonstrate the topological corner states in a large size twisted bilayer graphene (more than $20,000$ atoms). We take the open boundary condition preserving the mirror symmetry. Figure 3(a) shows the tight-binding model calculations of the localized fractional e/2 charge existing at the mirror symmetric corners. In addition, one can also consider the boundary termination with $C_{6}$ symmetry. In such case, we can define $\mathbb{Z}_{6}$ valued $C_{6}$ quadrupole moment. This value is explicitly given as Benalcazar et al. (2018),

[TABLE]

where $\#C_{3}$ counts the difference in the number of the occupied bands with $\hat{C}_{3}$ eigenvalue $1$ between $K$ and $\Gamma$. We find that the second term in Eq. (6) vanishes. Accordingly, $e/2$ localized corner charge occurs at each $C_{6}$ symmetric corners. Fig. 3 shows such localized corner states derived from the tight binding model.

Generalization to an arbitrary twist angle - It is important to note that the Stiefel-Whitney number, as well as the mirror winding number, is determined by the eigenstates at TRIM, whereas the low-energy physics of the twisted bilayer graphene is described in the vicinity $K$ and $K^{\prime}$ points in the Moiré BZ. Therefore, the interlayer coupling does not alter the occupations of the TRIM points as compared to the non-interacting monolayer graphenes. As a consequence, the eigenvalues of $C_{2z}$ and the Stiefel-Whitney number remains robust irrespective on the microscopic details of the interlayer coupling as long as the associated symmetries are preserved. Therefore, we can extend our analysis to arbitrary commensurate angles without the microscopic band structure calculations, assuming that all the low energy physics is described near $K$ and $K^{\prime}$ points. Using this special property, we now show that the twisted bilayer graphene in all other commensurate angles (arbitrary $p,q$) are higher-order topologically non-trivial at half-filling. To do so, we first need to count the eigenvalues of inversion symmetry for the non-interacting monolayer graphene at the TRIM points in the Moiré BZ. We first notice that $\frac{3p^{2}+3pq+q^{2}}{gcd(q,3)}\equiv 2N+1$ ($N\in\mathbb{Z}$) is always an odd integer and the area of the monolayer BZ fills $2N+1$ numbers of the Moiré BZ. For example, Fig. 4(b) shows that the monolayer BZ fills $7$ times of Moiré BZ at $\theta=21.78^{\circ}$. In generic commensurate angles, there exist $2N+1$ occupied states at each TRIM point, say ${\bm{X}}$, in the Moiré BZ. In the original BZ of a monolayer, these states correspond to the occupied state of the monolayer graphene at distinct ‘non-TRIM’ points related by the Moiré reciprocal vectors [See Fig. 4(b)]:

[TABLE]

where $m,n$ are some integers. Among all the possible ‘independent’ $2N+1$ numbers of momentum points ${\bm{X}}_{m,n}$, only ${\bm{X}}_{0,0}$ is the inversion symmetric point in the monolayer BZ. Other $2N$ distinct momenta, ${\bm{X}}_{m,n}$, are mapped to ${\bm{X}}_{-m,-n}$ under inversion symmetry. Explicitly, one can choose a gauge such that $\hat{C}_{2z}$ symmetry operator acts in the occupied eigenstates space as, $\hat{C}_{2z}|{\bm{X}_{m\neq 0,n\neq 0}}\rangle=|{\bm{X}}_{-m,-n}\rangle$. If we choose the basis of the occupied states at the TRIM points in the Moiré BZ as $(|{\bm{X}}_{0,0}\rangle,|{\bm{X}}_{0,1}\rangle,|{\bm{X}}_{-0,-1}\rangle,...,|{\bm{X}}_{m,n}\rangle,|{\bm{X}}_{-m,-n}\rangle)^{T}$, $\hat{C}_{2z}$ operator in the Moiré BZ is explicitly represented as,

[TABLE]

where $\zeta({\bm{X}})$ is the inversion eigenvalue of the occupied state at ${\bm{X}}$ in the monolayer graphene. Applying the eigenvalues of Eq. (10) to the second Stiefel-Whitney number in Eq. (3), the contributions from ${\bm{X}}_{m\neq 0,n\neq 0}$ cancels out and we find that

[TABLE]

In the monolayer graphene, the inversion eigenvalues at $M$ point and $\Gamma$ point always differ by the sign, so the second Stiefel-Whitney number must be always non-trivial, $\omega_{2}=1$. As a consequence, the twisted bilayer graphene at all commensurate angles must be higher-order topologically non-trivial at half filling, as long as the gap is present at $K$ point. In principle, the topological transition to the trivial state is possible if the additional gap closing occurs at $M$ point. However, such gap closing contradicts with the previous first-principle calculationsShallcross et al. (2010).

Conclusion - In summary, we have studied the higher-order topological properties of the twisted bilayer graphene. When the twisted angle is large enough, the Moiré pattern becomes short ranged, the intervalley scattering plays an important role. It results in $U(1)_{v}$ symmetry breaking and the twisted bilayer graphene system opens up the gap even at charge neutrality. Based on the first principle DFT calculation, we have demonstrated the existence of gap. We have also shown that the twisted bilayer graphene naturally hosts topological corner states that are protected by the mirror winding number. It is worth emphasizing the emergence of the HOTIs in twisted bilayer graphene at any commensurate angles is very generic properties induced by intervalley scattering and exists in any commensurate twist angles as long as the system opens a gap. In addition, the corner states occur at the half filling such that it does not require any fine tuning of chemical potential. Our work provides important guidance for the search of the higher-order topological materials and paves the way for future experiments in the twisted bilayer graphenes.

I Supplementary Material for “Higher-Order Topological Insulator in Twisted Bilayer Graphene”

I.1 Methods for first-principles calculations

Our first-principles calculations were carried out based on density functional theory (DFT) as implemented in the Quantum Espresso package Giannozzi et al. (2009). We employed local density approximation (LDA) with the Perdew-Zunger type exchange-correlation energy Perdew and Zunger (1981). A norm–conserving, optimized, designed nonlocal pseudopotential is generated for describing the core electrons of carbon using the Opium package Rappe et al. (1990); Ramer and Rappe (1999). The electronic wave functions were expanded in terms of a discrete set of plane-waves basis within the energy cutoff of 680 eV. The 30$\times$30$\times$1 ${\bm{k}}$-points were sampled from the first Brillouin zone (BZ) using the Monkhorst-Pack scheme Monkhorst and Pack (1976). The atomic structures were fully relaxed within a force tolerance of 0.005 eV/Å. The van der Waals interaction described described based on the semiempirical dispersion-correction DFT (DFT-D) method Grimme (2006).

I.2 Detailed Calculation of Mirror winding number

In this section, we provide the detailed calculation of the mirror winding number. We start with the Wilson lines in Eq. (4). We can related the Wilson loop and the inversion eigenvalues as following,

[TABLE]

where $u_{n}(\mathbf{k})$ represents the eigenstate of $n$-th band at the momentum $\mathbf{k}$. We also use the notation $\Gamma^{\prime}\equiv\Gamma+\mathbf{b_{1}}+\mathbf{b_{2}}$ to avoid the confusion between $\Gamma-M_{3}$ and $M_{3}-\Gamma^{\prime}$ Wilson lines. Since $\sum_{n\in occ}\zeta_{n}(M_{1})\zeta_{n}(M_{2})=1$ due to $C_{3}$ symmetry, we can rewrite the winding number as,

[TABLE]

Since $\nu_{+}+\nu_{-}=0~{}({\rm mod}~{}2)$, we formally equate the mirror winding number to the second Stiefel-Whitney number as follows.

[TABLE]