Quasicrystalline electronic states in 30$^\circ$ rotated twisted bilayer graphene

TL;DR

This paper develops a theoretical framework to describe the electronic states in 30° twisted bilayer graphene, revealing quasicrystalline properties, resonant states, and flat bands with fractal wave function patterns.

Contribution

It introduces a symmetry-respecting theoretical model for 30° twisted bilayer graphene, capturing its quasicrystalline electronic states and resonant phenomena.

Findings

Identification of 12-fold quantized angular momentum states

Presence of nearly flat bands with fractal wave functions

Spatial localization of resonant states in finite samples

Abstract

The recently realized bilayer graphene system with a twist angle of $30^\circ$ offers a new type of quasicrystal which unites the dodecagonal quasicrystalline nature and graphene's relativistic properties. Here, we introduce a concise theoretical framework that fully respects both the dodecagonal rotational symmetry and the massless Dirac nature, to describe the electronic states of the system. We find that the electronic spectrum consists of resonant states labeled by 12-fold quantized angular momentum, together with the extended relativistic states. The resulting quasi-band structure is composed of the nearly flat bands with spiky peaks in the density of states, where the wave functions exhibit characteristic patterns which fit to the fractal inflations of the quasicrystal tiling. We also demonstrate that the 12-fold resonant states appear as spatially-localized states in a finite-size geometry, which is another hallmark of quasicrystal. The theoretical method introduced here is applicable to a broad class of "extrinsic quasicrystals" composed of a pair of two-dimensional crystals overlaid on top of the other with incommensurate configurations.

Figures

Equations

$$\displaystyle{\bf R}_{X}=n_{1}{\bf a}_{1}+n_{2}{\bf a}_{2}+\mbox{\boldmath$\tau$}_{X}$$

$$\displaystyle{\bf R}_{\tilde{X}}=\tilde{n}_{1}\tilde{{\bf a}}_{1}+\tilde{n}_{2}\tilde{{\bf a}}_{2}+\mbox{\boldmath$\tau$}_{\tilde{X}}$$

$$|\mathbf{k},X⟩=\frac{1}{\sqrt{N}}\sum _{{\mathbf{R}}_X}{e}^{i\mathbf{k}\cdot {\mathbf{R}}_X}|{\mathbf{R}}_X⟩(\text{mbox}layer1),$$

$$|\tilde{\mathbf{k}},\tilde{X}⟩=\frac{1}{\sqrt{N}}\sum _{{\mathbf{R}}_{\tilde{X}}}{e}^{i\tilde{\mathbf{k}}\cdot {\mathbf{R}}_{\tilde{X}}}|{\mathbf{R}}_{\tilde{X}}⟩(\text{mbox}layer2),$$

$$-T(\mathbf{R})=V_{pp\pi }\left[1-{\left(\frac{\mathbf{R}\cdot {\mathbf{e}}_z}{R}\right)}^2\right]+V_{pp\sigma }{\left(\frac{\mathbf{R}\cdot {\mathbf{e}}_z}{R}\right)}^2,$$

$$V_{pp\pi }=V_{pp\pi }^0{e}^{-(R-a/\sqrt{3})/r_0},V_{pp\sigma }=V_{pp\sigma }^0{e}^{-(R-d)/r_0},$$

$$⟨{\mathbf{k}}^{\prime },X^{\prime }|H_1|\mathbf{k},X⟩=h_{X,X^{\prime }}(\mathbf{k}){\delta }_{{\mathbf{k}}^{\prime },\mathbf{k}},$$

$$\displaystyle h_{X,X^{\prime}}({\bf k})=\sum_{{\bf L}}-T({\bf L}+\mbox{\boldmath$\tau$}_{X^{\prime}X})e^{-i{\bf k}\cdot({\bf L}+\mbox{\boldmath$\tau$}_{X^{\prime}X})},$$

$$\displaystyle\langle\tilde{{\bf k}},\tilde{X}|U|{\bf k},X\rangle=-\sum_{{\bf G},\tilde{{\bf G}}}{t}({\bf k}+{\bf G})e^{-i{\bf G}\cdot\mbox{\boldmath\scriptsize$\tau$}_{X}+i\tilde{{\bf G}}\cdot\mbox{\boldmath\scriptsize$\tau$}_{\tilde{X}}}\,\delta_{{\bf k}+{\bf G},\tilde{{\bf k}}+\tilde{{\bf G}}},$$

$$t(\mathbf{q})=\frac{1}{S}\int T(\mathbf{r}+z_{\tilde{X}X}{\mathbf{e}}_z)e^{-i\mathbf{q}\cdot \mathbf{r}}d\mathbf{r}$$

$${\mathcal{H}}_{\mathrm{ring}}({\mathbf{k}}_0)=\left(\begin{array}{cccccc}{H}^{(0)}& W^{†}& & & & W\\ W& H^{(1)}& W^{†}& & & \\ & W& H^{(2)}& W^{†}& & \\ & & \ddots & \ddots & \ddots & \\ & & & W& H^{(10)}& W^{†}\\ W^{†}& & & & W& H^{(11)}\end{array}\right),$$

$$H^{(n)}({\mathbf{k}}_0)=\left(\begin{array}{cc}{h}_{AA}^{(n)}& h_{AB}^{(n)}\\ h_{BA}^{(n)}& h_{BB}^{(n)}\end{array}\right),W=-t_0\left(\begin{array}{cc}\omega & 1\\ 1& {\omega }^{\ast }\end{array}\right),$$

$$E_m^{\pm }=t_0\cos q_m\pm \sqrt{3t_0^2{\sin }^2{q}_m+(h_0-2t_0\cos q_m{)}^2},$$

$$P^{-1}(\psi )=\frac{{\sum }_i^{\prime }|{\psi }_i{|}^4}{{\left({\sum }_i^{\prime }|{\psi }_i{|}^2\right)}^2}.$$

$${\mathbf{k}}_1$$

$${\tilde{\mathbf{k}}}_1$$

$$|\mathbf{q}(={\mathbf{k}}_0+{\mathbf{G}}_1+{\tilde{\mathbf{G}}}_1)|\le \mathcal{O}(|{\mathbf{a}}_i^{\ast }|).$$

$$|{\mathbf{k}}_2|\le \mathcal{O}(|{\mathbf{a}}_i^{\ast }|).$$

$${\tilde{\mathbf{k}}}_2\equiv {\mathbf{k}}_0+{\mathbf{G}}_1+{\tilde{\mathbf{G}}}_1.$$

$${\mathbf{k}}_0^{\prime }\equiv {\mathbf{k}}_2,$$

$${\mathbf{k}}_2$$

$${\tilde{\mathbf{k}}}_2$$

$${\mathbf{k}}_2(={\mathbf{k}}_0+{\tilde{\mathbf{G}}}_1-{\mathbf{G}}_2)$$

$${\tilde{\mathbf{k}}}_2(={\mathbf{k}}_0+{\mathbf{G}}_1+{\tilde{\mathbf{G}}}_1)$$

$${\mathbf{q}}^{\prime }\equiv {\mathbf{k}}_2-({\mathbf{k}}_0^{\prime }-{\tilde{\mathbf{k}}}_2)=\mathbf{q}.$$

Full text

I Introduction

When two graphene layers are overlapped on top of the other, the interlayer twist angle $\theta$ is an important physical quantity to determine the electronic structures. This twisted bilayer graphene (TBG) is essentially a quasi-periodic system, as the two lattice periods of individual graphene layers are generally irrational to each other. When $\theta$ is relatively small (less than about 10∘), however, the low-energy physics is governed by the long-range moiré interference pattern, and then the electronic properties are captured by the moiré effective theory that does not need an exact lattice matching. In brief, the effective theory approximately treats TBG as a translationally-symmetric system ruled by the moiré period. The exotic phenomena in the low-angle regime Mele (2010); Fu et al. (2018), such as the flat band formation Lopes Dos Santos et al. (2007); Trambly de Laissardière et al. (2010); Shallcross et al. (2010); Suárez Morell et al. (2010); Bistritzer and MacDonald (2011); Moon and Koshino (2012); Trambly de Laissardière et al. (2012) and the Hofstadter butterfly under magnetic field Moon and Koshino (2012); Dean et al. (2013); Hunt et al. (2013); Ponomarenko et al. (2013); Moon and Koshino (2014), can be understood in terms of the moiré effective theory.

In TBG of large $\theta$, on the other hand, the moiré period competes with the atomic length scale and the quasi-periodic nature emerges Koren and Duerig (2016). When $\theta=30^{\circ}$, in particular, the overlaid two hexagonal lattices is mapped onto a 12-fold rotationally symmetric quasicrystalline lattice without any translational symmetry [Fig. 1(a)], as first shown by Stampfli Stampfli (1986). Recently, the TBG with a precise rotation angle of $30^{\circ}$ was experimentally realized and its spectrum measured in epitaxially grown samples on top of SiC surface Ahn et al. (2018). In addition, similar TBGs have been realized on top of Ni surface Takesaki et al. (2016); Yao et al. (2018) and also by a transfer method Chen et al. (2016). Moreover, another $30^{\circ}$-rotated stack of atomic layers have also been realized in graphene on top of BN layer Wang et al. (2016) as well as MoSe2 bilayer system Choi et al. (2017). In such the quasicrystalline TBG (QC-TBG), the moiré effective approach sketched above breaks down because its main assumption that the moiré pattern governs the system is no longer valid.

In the literature, several theoretical approaches have been applied to understand the electronic structures of conventional quasi-periodic systems Roche et al. (1997), such as one-dimensional Fibonacci lattices Niu and Nori (1986); Kohmoto et al. (1987), two-dimensional non-periodic tiling including Penrose lattice Niizeki and Akamuatsu (1990); Gambaudo and Vignolo (2014), metal nanoparticles Dong et al. (2009), photonic quasicrystals Mnaymneh and Gauthier (2007) and three-dimensional alloys including Al-Mn, Al-Ni-Co, and Al-Cu-Co Smith and Ashcroft (1987); Fujiwara and Yokokawa (1991); Hafner and Krajčí (1992); Trambly de Laissardière and Fujiwara (1994); Roche and Fujiwara (1998); Rogalev et al. (2015). These systems can be viewed as intrinsic quasicrystals where the atomic sites are intrinsically arranged in the quasi-periodic order. In contrast, the QC-TBG is regarded as an extrinsic quasicrystal, in that it is composed of a pair of perfect crystals having independent periodicities, and the quasi-periodic nature appears only in the perturbational coupling between the two subsystems. Thus, the QC-TBG unites the quasicrystalline order and the relativistic nature of the massless Dirac particles of graphene, yet it is not obvious whether and in what form the essential features of quasicrystals emerge in the electronic properties. Since such a hybrid situation is out of the scope of the previous theories of intrinsic quasicrystals, we need an alternative theoretical framework to properly describe the quasicrystalline physics of QC-TBG.

In this paper, we develop a concise model Hamiltonian that fully respects both the dodecagonal rotational symmetry and the massless Dirac nature, to describe the quasicrystalline electronic states in the QC-TBG. We find that the electronic spectrum of QC-TBG is characterized by the 12-fold resonant states of relativistic Dirac fermions, and they can be well captured by a ring Hamiltonian composed of 12 Dirac cones. The resulting quasi-band structure comprises a series of the nearly flat bands corresponding to the resonant states, each of which is labeled by a 12-fold quantized angular momentum. The spatial pattern of wave functions exhibit the fractal inflations of the Stampfli tiling, which is a direct manifestation of the quasicrystalline nature Niizeki (1989). Since we can tune the twist angles in the model, the transition of electronic states from the approximants Goldman and Kelton (1993) of QC-TBG to a true dodecagonal rotational symmetry can be continuously described within a 12-fold ring model, and the emergence of quasicrystalline states and the validity of the approximant method are critically attested. We also show that the 12-fold resonant states appear as spatially-localized states in a finite-size geometry, which is another hallmark of the quasicrystalline nature Niu and Nori (1986); Kohmoto et al. (1987); Deguchi et al. (2012). The proposed theoretical approach is applicable to a broad class of extrinsic quasicrystals, and its simple structure of the closed Hamiltonian allows rigorous analysis on exotic quantum phenomena of quasicrystals.

The paper is organized as follows. In Sec. II, we present the tight-binding model for QC-TBG, and introduce the dual tight-binding approach in the momentum space. In Sec. III.1, we derive the approximate 12-wave ring Hamiltonian, and using this, we describe the quasi-band structure, the resonant states and the characteristic wave functions to respect the Stampfli tiling. In Sec. III.3, we calculate the electronic states of QC-TBG in an alternative method using the finite-size tight-binding model, and demonstrate the localization nature of the 12-fold resonant states. A brief conclusion is given in Sec. IV.

II Theoretical methods

II.1 Tight-binding Hamiltonian for QC-TBG

We define the atomic structure of QC-TBG by starting from AA-stacked bilayer graphene (i.e. perfectly overlapping honeycomb lattices) and rotating the layer 2 around the center of hexagon by $30^{\circ}$ [Fig. 1(a)]. We set $xy$ coordinates parallel to the graphene layers and $z$ axis perpendicular to the plane. The system belongs to the symmetry group $D_{6d}$, and it is invariant under an improper rotation $R(\pi/6)M_{z}$, where $R(\theta)$ is the rotation by an angle $\theta$ around $z$ axis, and $M_{z}$ is the mirror reflection with respect to $xy$ plane. The primitive lattice vectors of layer 1 are taken as ${\bf a}_{1}=a(1,0)$ and ${\bf a}_{2}=a(1/2,\sqrt{3}/2)$ with the lattice constant $a\approx 0.246\,\mathrm{nm}$, and those of the layer 2 as $\tilde{{\bf a}}_{i}=R(\pi/6)\,{\bf a}_{i}$. Accordingly, the reciprocal lattice vectors of layer 1 are given by ${\bf a}^{*}_{1}=(2\pi/a)(1,-1/\sqrt{3})$ and ${\bf a}^{*}_{2}=(2\pi/a)(0,2/\sqrt{3})$, and layer 2 by $\tilde{{\bf a}}^{*}_{i}=R(\pi/6)\,{\bf a}^{*}_{i}$. The atomic positions are given by

[TABLE]

where $n_{i}$ and $\tilde{n}_{i}$ are integers, $X=A,B$ ($\tilde{X}=\tilde{A},\tilde{B}$) denotes the sublattice site of layer 1(2), and \mbox{\boldmath\tau}_{X} and \mbox{\boldmath\tau}_{\tilde{X}} are the sublattice positions in the unit cell, defined by \mbox{\boldmath\tau}_{A}=-\mbox{\boldmath\tau}_{1}, \mbox{\boldmath\tau}_{B}=\mbox{\boldmath\tau}_{1}, \mbox{\boldmath\tau}_{\tilde{A}}=-R(\pi/6)\mbox{\boldmath\tau}_{1}+d{\bf e}_{z}, \mbox{\boldmath\tau}_{\tilde{B}}=R(\pi/6)\mbox{\boldmath\tau}_{1}+d{\bf e}_{z} with \mbox{\boldmath\tau}_{1}=(0,a/\sqrt{3}). Here $d\approx 0.335\,\mathrm{nm}$ is the interlayer spacing between graphene layers and ${\bf e}_{z}$ is the unit vector normal to the layer.

We model graphene by the tight-binding model of carbon $p_{z}$ orbitals. The Hamiltonian is spanned by the Bloch bases of $p_{z}$ orbitals at difference sublattices,

[TABLE]

where $|{\bf R}_{X}\rangle$ is the atomic $p_{z}$ orbital at the site ${\bf R}_{X}$, ${\bf k}$ and $\tilde{{\bf k}}$ are the two-dimensional Bloch wave vectors and $N=S/S_{\rm tot}$ is the number of graphene’s unit cells $S=(\sqrt{3}/2)a^{2}$ in the total system area $S_{\rm tot}$. We assume that the transfer integral between any two $p_{z}$ orbitals is expressed as Slater and Koster (1954)

[TABLE]

where ${\bf R}$ is the relative vector between two atoms, $V_{pp\pi}^{0}\approx-2.7\,\mathrm{eV}$, $V_{pp\sigma}^{0}\approx 0.48\,\mathrm{eV}$, and $r_{0}\approx 0.0453\,\mathrm{nm}$ Trambly de Laissardière et al. (2010); Moon and Koshino (2013a).

The total tight-binding Hamiltonian is expressed as $H=H_{1}+H_{2}+U$ where $H_{1}$ and $H_{2}$ are the Hamiltonian for the intrinsic monolayer graphenes of layer 1 and 2, respectively, and $U$ is for the interlayer coupling. The intralayer matrix elements of layer 1 are given by

[TABLE]

where ${\bf L}=n_{1}{\bf a}_{1}+n_{2}{\bf a}_{2}$ and \mbox{\boldmath\tau}_{X^{\prime}X}=\mbox{\boldmath\tau}_{X^{\prime}}-\mbox{\boldmath\tau}_{X}. Similarly, the matrix for $H_{2}$ is given by replacing ${\bf k}$ with $R(-\pi/6){\bf k}$.

The interlayer matrix element between layer 1 and 2 is written as Mele (2010); Bistritzer and MacDonald (2011); Koshino (2015)

[TABLE]

where ${\bf G}=m_{1}{\bf a}^{*}_{1}+m_{2}{\bf a}^{*}_{2}$ and $\tilde{{\bf G}}=\tilde{m}_{1}\tilde{{\bf a}}^{*}_{1}+\tilde{m}_{2}\tilde{{\bf a}}^{*}_{2}$ ($m_{1},m_{2},\tilde{m}_{1},\tilde{m}_{2}\in\mathbb{Z}$) run over all the reciprocal points of layer 1 and 2, respectively. We also defined

[TABLE]

where z_{\tilde{X}X}=(\mbox{\boldmath\tau}_{\tilde{X}}-\mbox{\boldmath\tau}_{X})\cdot{\bf e}_{z}.

II.2 Dual tight-binding lattice in momentum space

Equation (5) shows that the interlayer interaction occurs between the states satisfying the generalized Umklapp scattering condition ${\bf k}+{\bf G}=\tilde{{\bf k}}+\tilde{{\bf G}}$. When we start from the layer 1’s Bloch states at ${\bf k}_{0}$, for example, the interlayer Hamiltonian $U$ couples this state with layer 2’s Bloch states at $\tilde{{\bf k}}={\bf k}_{0}+{\bf G}-\tilde{{\bf G}}$. They are further coupled back to layer 1’s states at ${\bf k}={\bf k}_{0}+{\bf G}^{\prime}-\tilde{{\bf G}}^{\prime}$, and so forth. As a result, the space of the wave functions associated with ${\bf k}_{0}$ is spanned by $\{|{\bf k},X\rangle\,|\,{\bf k}={\bf k}_{0}+\tilde{{\bf G}}-{\bf G}\}$ and $\{|\tilde{{\bf k}},\tilde{X}\rangle\,|\,\tilde{{\bf k}}={\bf k}_{0}+{\bf G}-\tilde{{\bf G}}\}$ for $\forall{\bf G}$ and $\forall\tilde{{\bf G}}$. However, we actually need only a subset of these sets, since the BZ of each layer is translationally invariant with respect to the reciprocal lattice vectors of its own (i.e., ${\bf k}$ and ${\bf k}+{\bf G}$ stand for the same Bloch wavenumber of layer 1). Thus, without loss of generality, we can choose the subspace spanned by the QC-TBG Hamiltonian as $\{|{\bf k},X\rangle\,|\,{\bf k}={\bf k}_{0}+\tilde{{\bf G}},\forall\tilde{{\bf G}}\}$ and $\{|\tilde{{\bf k}},\tilde{X}\rangle\,|\,\tilde{{\bf k}}={\bf k}_{0}+{\bf G},\forall{\bf G}\}$. Here note that the $k$-points in each layer is regularly spaced with the reciprocal vectors of the other layer.

According to Eq. (5), the interaction strength between ${\bf k}={\bf k}_{0}+\tilde{{\bf G}}$ and $\tilde{{\bf k}}={\bf k}_{0}+{\bf G}$ is given by $t({\bf q})$ where ${\bf q}={\bf k}+{\bf G}=\tilde{{\bf k}}+\tilde{{\bf G}}={\bf k}+\tilde{{\bf k}}-{\bf k}_{0}$. Since $t({\bf q})$ decays in large ${\bf q}$, the relevant contribution occurs only when $|{\bf k}+\tilde{{\bf k}}-{\bf k}_{0}|$ is relatively small. The interaction strength can be visualized by the diagram Fig. 1(b), where all the layer 2’s wave points $\tilde{{\bf k}}$ are inverted to ${\bf k}_{0}-\tilde{{\bf k}}$, and overlapped with the layer 1’s wave points ${\bf k}$. In the map, the quantity $|{\bf k}+\tilde{{\bf k}}-{\bf k}_{0}|$ is the geometrical distance by given two points, so that the interaction takes place only between the points located in close distance. If the $k$-points are viewed as ‘sites’, the whole system can be recognized as a tight-binding lattice in $k$-space, which is the dual counterpart of the original tight-binding Hamiltonian in the real space. It should be noted that, unlike the real-space version, the intralayer Hamiltonians $H_{1}$ and $H_{2}$ now can be interpreted as $k$-dependent on-site potential in the $k$-space, which is nothing but the band energy of intrinsic graphene. Recently, the relationship between the real space and the momentum space was also noticed in the localized wave functions in moiré bilayer systems Carr et al. (2018).

In this $k$-space tight-binding model, the hopping between different $k$-space sites (the interlayer interaction $U$) is smaller by an order of magnitude than the potential landscape (the band energy), so that the eigen functions tend to be localized in the $k$-space lattice, in a similar manner to the Aubry-André model in one dimension Aubry and André (1980). In the practical calculation, therefore, we are allowed to take only a limited number of wave points around ${\bf k}_{0}$ inside a certain cut-off circle, and obtain the energy eigenvalues by diagonalizing the Hamiltonian matrix within the finite bases. If we plot the energy levels against ${\bf k}_{0}$, we obtain the quasi-band structures of the system. Here the wavenumber ${\bf k}_{0}$ works like the crystal momentum for the periodic system, so it can be called the quasicrystal momentum. The cut-off radius $k_{c}$ should be greater than the typical localization length in the $k$-space, but need not be too large, since the wave points discarded outside $k_{c}$ are properly considered by shifting ${\bf k}_{0}$. If we increase $k_{c}$, we will see more and more replicas of the identical quasi-energy band with different origins, because shifting ${\bf k}_{0}$ actually corresponds to taking a different origin in the $k$-space map of Fig. 1(b). The resonant band structure near ${\bf k}_{0}=0$ barely changes in this process because its wave function is very well localized to the 12-membered ring in the $k$-space. The replica bands are just duplication of the same states so they should be appropriately removed in calculating the physical quantities such as the density of states. The validity of the momentum cut-off is discussed in detail in Appendix A.

III Results and discussion

III.1 12-fold symmetric resonant states

At ${\bf k}_{0}={\bf 0}$, we see that the twelve symmetric points ${\bf Q}_{n}=a^{*}[\cos(n\pi/6),\sin(n\pi/6)](n=0,1,2,\cdots,11)$ form a circular chain in the dual-tight-binding lattice of which radius is $a^{*}\equiv|{\bf a}^{*}_{i}|=4\pi/(\sqrt{3}a)$, indicated by the dashed ring in Fig. 1(b). Noting that the layer 2’s wave points are inverted, these points are associated with layer 1’s Bloch wavenumbers ${\bf k}={\bf Q}_{n}$ for even $n$’s and layer 2’s $\tilde{{\bf k}}=-{\bf Q}_{n}$ for odd $n$’s. Figure 1(c) shows the original positions of ${\bf k}$ (layer 1) and $\tilde{{\bf k}}$ (layer 2) associated with ${\bf Q}_{n}$, in the first Brillouin zone. Due to the symmetry, the intrinsic graphene’s Bloch states at the twelve points are all degenerate in energy, and therefore the interlayer coupling hybridizes them to make resonant states. Here the coupling is only relevant between the neighboring sites of the ring, and it is given by $t_{0}=t(2a^{*}\sin 15^{\circ})\approx 157\,\mathrm{meV}$.

In the vicinity of ${\bf k}_{0}={\bf 0}$, the Hamiltonian of the ring can be expressed by a $24\times 24$ matrix,

[TABLE]

where $h_{X^{\prime}X}^{(n)}({\bf k}_{0})=h_{X^{\prime}X}[R(-7n\pi/6){\bf k}_{0}+{\bf Q}_{0}]$, $\omega=e^{2\pi i/3}$, and we neglect the ${\bf k}_{0}$ dependence of the interlayer matrix element $t({\bf q})$. The diagonal block $H^{(n)}$ represent monolayer’s Hamiltonian at ${\bf k}={\bf k}_{0}+{\bf Q}_{n}$ for even $n$ (layer 1) and $\tilde{{\bf k}}={\bf k}_{0}-{\bf Q}_{n}$ for odd $n$ (layer 2). In each $2\times 2$ block the sublattices are arranged in the order of $(A,B)$ or $(\tilde{A},\tilde{B})$ for $n\equiv 0,3$ in modulo of 4, and $(B,A)$ or $(\tilde{B},\tilde{A})$ for $n\equiv 1,2$. By doing this, the first base of a $2\times 2$ block is always mapped to the first base of other block under the operation of $R(\pi/6)M_{z}$. Note that the arrangement of $h_{AA}$, $h_{AB}$, etc. in the submatrix $H^{(n)}$ is fixed irrespective of $n$, and the dependence on $n$ solely comes from $R(-7n\pi/6){\bf k}_{0}$ in the argument of $h_{X^{\prime}X}$. Consequently, the total Hamiltonian ${\cal H}_{\rm ring}$ is obviously symmetric under rotation by a single span of the ring (i.e., moving ${\bf Q}_{n}$ to ${\bf Q}_{n+1}$), which actually corresponds to the operation $[R(\pi/6)M_{z}]^{7}$ (210∘ rotation and swapping layer 1 and 2) in the original system.

Figures 2 show (a) the density of states (DOS), (b) the band structures as a function of ${\bf k}_{\bf 0}$ in the negative energy region, and (c) its closer view near ${\bf k}_{0}={\bf 0}$. The twelve Dirac cones are arranged on a circle with a radius $\Delta k=4(2-\sqrt{3})\pi/(3a)$ and they are strongly mixed near ${\bf k}_{0}={\bf 0}$. As a result, the originally degenerate twelve states of graphene (in each of the electron side and the hole side of the Dirac cone) split into different energies, and exhibit the characteristic dispersion including flat band-bottoms and the Mexican-hat edges. This leads to a series of spiky peaks and dips (pseudogaps) in DOS. At ${\bf k}_{0}={\bf 0}$, the Hamiltonian can be analytically diagonalized to obtain a set of energies (neglecting the constant energy)

[TABLE]

where $h_{0}=h_{AB}({\bf Q}_{0})=h_{BA}({\bf Q}_{0})=$ 1.84 eV, $\pm$ corresponds to the conduction band and valence band, respectively, and $q_{m}=(7\pi/6)m$ with $m=-5,-4,\cdots,5,6$ is the wavenumber along the chain. The eigenvalue of $R(\pi/6)M_{z}$ is given by $e^{im\pi/6}$. Here the states with $m=\pm s$ $(s=1,2,3,4,5)$ form twofold doublets, and belong to two-dimensional $E_{s}$ irreducible representation of $D_{6d}$ point group. The $m=0$ and $6$ are non-degenerate, and belong to $A_{1}(A_{2})$ and $B_{2}(B_{1})$, respectively, for the conduction (valence) band. If we disregard the $z$-position difference, the index $m$ is regarded as quantized angular momentum in 12-fold rotational symmetry, and this is an essential characteristics of quasicrystal TBG.

We have similar resonant states also in the conduction band, while the energy scale of the band structures is much smaller than in the valence band. Equation (9) clearly explains such asymmetry; the dispersion of $E^{-}_{m}$ in $q_{m}$ is nearly three times wider than that of $E^{+}_{m}$, considering that $h_{0}\gg t_{0}$. Intuitively, the wave function of the conduction band of intrinsic graphene has the opposite phases between the sublattice $A$ and $B$, and then the interlayer coupling between incommensurate layers is tend to be suppressed by the phase cancellation.

III.2 Wave functions showing the quasicrystal tiling

The 12-wave resonant coupling also gives rise to a characteristic pattern in the wave function. Figure 2(d) shows the wave functions at ${\bf k}_{0}={\bf 0}$ where the hybridization is the most prominent, where we can see that the wave amplitude distribute selectively on a limited number of sites in a 12-fold rotationally symmetric pattern. The extent of the hybridization of different wave modes is characterized by the inverse participation ratio (IPR) on the dominant layer,

[TABLE]

Here $\psi_{i}$ is the amplitude at the site $i$ of the eigenstates $\psi$, and $\sum^{\prime}_{i}$ represents the sum over the sites on the dominant layer, which is defined as the layer having greater wave amplitude than the other. We have $P^{-1}=1$ for a pure single layer state, and $P^{-1}=1.5$ for a hybrid state of two plain wave modes. In Fig. 2(b), the blue dots represent IPR at several sample points in the band structures along $x$ direction, where the dot area is proportional to $P^{-1}-1$. We see that the IPR becomes large exclusively around ${\bf k}_{0}={\bf 0}$, where the 12 wave components are strongly hybridized. $P^{-1}$ remains almost 1 near the Dirac cones where the hybridization is almost negligible. We also have a region of $P^{-1}\sim 1.5$ along the arch-shaped gap below the Dirac cone, which is indicated by “2-wave” in Fig. 3(b). These “2-wave” states arise from the hybridization of the $K$ and $K^{\prime}$ of the same layer assisted by the second-order process of the interlayer coupling $U$.

We also show a large scale plot of $m=0,6$ states in Fig. 3(a). We see that the wave pattern perfectly follows the Stampfli tiling, where the red (left-half) and black (right-half) lines represent the third and fourth generations of the fractal inflation, respectively Stampfli (1986). Such the long-range structure of the quasicrystalline wave function is actually quite sensitive to a slight change of the twist angle. Figure 3(b) represents the wave pattern of the corresponding state in TBG with $\theta=29.84^{\circ}$, calculated by the same 12-wave method. The TBG of $29.84^{\circ}$ is a quasicrystal approximant, which is not quasi-periodic but has a translational symmetry with period of $3.31\,\mathrm{nm}$. We can see that the local wave pattern is quite similar to that of $30^{\circ}$, while the long-range quasi-periodic nature is completely lost and round to a periodic pattern. Here we confirmed that the quasi-band structures, DOS and IPR look almost the same as $30^{\circ}$, but the tiny change of the wave bases and the coupling matrix elements in the 12-ring Hamiltonian encodes the periodic / quasi-periodic transition.

The energy spectrum of the QC-TBG approximant can also be calculated by the original real-space tight-binding model since it has a finite superlattice unit cell. We can show that the DOS and the wave function of $29.84^{\circ}$ calculated by the original tight-binding model are virtually the same as the result of the 12-ring effective Hamiltonian, and this justifies the validity the effective approach. In the appendix B, we present the extensive study on the electronic structures of the quasi approximants in all the angle region from $0^{\circ}$ to $30^{\circ}$.

III.3 Localization in finite-sized QC-TBG

The emergence of quasicrystalline states in QC-TBG can also be confirmed by a finite-sized tight-binding lattice, while the computation is enormous. Here we consider a tight-binding lattice composed of two large disks of graphene with radius $R=39.4\,\mathrm{nm}$ stacked at exactly $30^{\circ}$, and calculate its electronic structures by diagonalizing the huge Hamiltonian matrix with the total number of atoms 371,532. As shown in Fig. 4(a), the DOS of the finite flakes (thick red line), which is obtained by broadening its discrete spectrum, is consistent with the DOS of the 12-wave effective model (thin black line) calculated by the effective Hamiltonian with a few wave bases [Fig. 2(b)]. In Fig. 4(b), we also present the wave functions at three energies, $\alpha$, $\beta$ and $\gamma$, which correspond to the band edges of the quasi-band structures in the effective Hamiltonian (Fig. 2).

The magnified plot of $\gamma$ is presented in the inset of Fig. 4(b) showing the characteristic pattern of 12-wave approximation. Interestingly, however, it is overlapped with an envelope function decaying in the radial direction. Such the localized feature is never seen in single layer graphene and it is the characteristics of the resonant states of QC-TBG. As a measure of the concentration to the center, we calculate the second momentum $\langle r^{2}\rangle/R^{2}$ for each eigenstate and plot it as green circles in Fig. 4(a). For a uniform state (i.e., the wave amplitude is constant throughout the system), $\langle r^{2}\rangle/R^{2}$ approaches $1/(2\pi)$, which is indicated by the dashed line. We can actually see that $\langle r^{2}\rangle/R^{2}$ lies around this line for most of the states, while it becomes exceptionally small at the energies of the quasi-band edges argued in the previous section. In terms of the quasi-band structures, these localized states actually correspond to the integral of the quasi-band states over the nearly flat region, and the length scale of the envelope function is related to the size of the flat area in the momentum space. For the state at $\gamma$ ($m=0$ state), for instance, the radius of the flat area is roughly given by $\delta k\sim 0.2/a$, and the corresponding real-space scale $r=2\pi/k\sim 7.7\,\mathrm{nm}$ matches the characteristic decaying and oscillating scale of the envelope function.

IV Conclusions

We revealed that the quasicrystalline nature emerges in the electronic properties of QC-TBG, or the twisted bilayer graphene stacked at 30∘. We developed a concise model Hamiltonian for this unique system, and demonstrated that the electronic structure is well described by the quasi-band picture despite of the lack of periodicity. The quasi-band states of the QC-TBG are characterized by the 12-fold resonant states of relativistic Dirac fermions, where the wave functions exhibit the spatial pattern fully respecting the dodecagonal quasicrystal tiling. Such a non-uniform distribution of electron may be observed by microscopy imaging techniques. The emergence of quasicrystalline states was attested by comparing the QC-TBG and a periodic approximant near $30^{\circ}$, and it was demonstrated that even a slight deviation from the QC configuration destroys the long-range quasicrystalline nature. Finally, we studied the electronic states of QC-TBG using the finite-size tight-binding model, where the 12-fold resonant states appear as spatially-localized states in a finite-size geometry.

While we considered the QC-TBG as a model example in this paper, the theoretical method based on the $k$-space tight-binding approach introduced here is applicable to any kind of extrinsic quasicrystals composed two-dimensional materials overlaid in incommensurate configurations, including heterostructures of two-dimensional materials having difference lattice symmetries (e.g., rectangle and hexagon).

Extrinsic quasicrystals also provide a unique opportunity to tune the quasicrystal bands by controlling the interlayer interaction strength $U$. As $U$ is an exponential function of the interlayer spacing $d$, Koshino et al. (2015) we can either increase $U$ by applying pressure, or decrease it through intercalation of ions or addition of barrier atomic layers Chittari et al. (2018). When $U$ becomes comparable to the width of the energy bands, we expect a transition from the weakly coupled regime to the strongly coupled regime where the quasicrystalline nature is even more pronounced. The detailed studies on exotic electronic natures in a broad class of extrinsic quasicrystals, such as the electronic transport, optical properties, the quantum Hall effect, and also the effects of $U$ modulation to these phenomena, are left for future research.

Acknowledgments

We thank L. A. Wray and A. Kent for fruitful discussions. P.M. was supported by NYU Shanghai (Start-Up Funds), NYU-ECNU Institute of Physics at NYU Shanghai, New York University Global Seed Grants for Collaborative Research. This research was carried out on the High Performance Computing resources at NYU Shanghai and CAC of KIAS. M.K. was supported by JSPS KAKENHI Grant Numbers JP25107005, JP15K21722, JP17K05496. Y.-W.S was supported by NRF of Korea (Grant No. 2017R1A5A1014862, SRC program: vdWMRC center).

Appendix A Validity of the momentum space cut-off

In this section, we argue about the validity of introducing the momentum space cut-off in calculating the quasi band structure. As we mentioned in Sec. II.2, the wave functions of the QC-TBG are localized in the $k$-space in a similar manner to the Aubry-André model in one-dimension Aubry and André (1980), because the hopping term in the $k$-space is much smaller than the potential landscape (the band energy). Figures 5(b)-(d) show some examples of the $k$-space amplitude map. The panel (b) shows a nearly decoupled state which is dominated by only a single state of monolayer graphene, and (c) is a state originating from the 2-wave mixing, where a pair of monolayer’s states on layer 2 are coupled though the mediation of a middle state on layer 1. The panel (d) is the 12-wave resonant state. Any eigenstates other than those examples are also localized within just a few reciprocal lattice constants in $k$-space. When we increase the number of total wave bases components (${|{\bf k},X\rangle}$ and ${|\tilde{{\bf k}},\tilde{X}\rangle}$) by increasing $k_{c}$, each eigenstate hardly changes as long as $k_{c}$ is greater than the typical localization length. In Fig. 5(a), we show the histogram of the number of dominant wave components in the eigenstates at a particular ${\bf k}_{0}$, calculated in the basis sets within $k_{c}$ of $150/a$, $75.1/a$ and $18.8/a$, respectively, where the total number of wave bases are $11630$, $2918$ and $182$, respectively. We actually see that each of eigenstates is composed only a few (mostly less than 10) bases. We note that, in large $k_{c}$, we often see a resonance between different localized states which are very distant in $k$-space. This does not much affect the calculation of the physical quantity because the overlap of the different localized wave functions are exponentially small.

We have infinitely many localized states far away from the first Brillouin zone, so one might think that it is necessary to take an infinite $k_{c}$ to properly include all the states. Note that, however, these localized states can be moved into the vicinity of the first Brillouin zone by shifting ${\bf k}_{0}$ with a proper amount, as we show in the following. Thus, instead of using a large $k_{c}$ requiring a large computational cost, we can obtain the full spectrum of the system by calculating the electronic structures as a function of ${\bf k}_{0}$ with a moderate $k_{c}$.

Let us consider two states $|{\bf k}_{1},X\rangle$ and $|\tilde{{\bf k}}_{1},\tilde{X}\rangle$ with

[TABLE]

for a given ${\bf k}_{0}$. Suppose ${\bf k}_{1}$ and $\tilde{{\bf k}}_{1}$ are outside the cut-off circle, i.e., $|{\bf k}_{1}|>k_{c}$ and $|\tilde{{\bf k}}_{1}|>k_{c}$, but they strongly interact with each other, i.e.,

[TABLE]

Now, for any such ${\bf k}_{1}$, we can always find ${\bf G}_{2}$ (${\bf G}_{2}\in{\bf G}$) which makes ${\bf k}_{1}$ move to the point ${\bf k}_{2}\equiv{\bf k}_{1}-{\bf G}_{2}$ in the first Brillouin zone, i.e.,

[TABLE]

And suppose $\tilde{{\bf k}}_{2}$, defined as

[TABLE]

Then, by shifting ${\bf k}_{0}$ to a new point ${\bf k}_{0}^{\prime}$ defined as

[TABLE]

we can see that

[TABLE]

are the member of the subspace spanned from ${\bf k}_{0}^{\prime}$. And by considering that $\tilde{{\bf k}}_{2}={\bf q}$, and from Eqs. 12 and 13, we can show that these two points are within the cut-off circle. Since

[TABLE]

$|{\bf k}_{2},X\rangle$ and $|\tilde{{\bf k}}_{2},\tilde{X}\rangle$ represent the Bloch states same to $|{\bf k}_{1},X\rangle$ and $|\tilde{{\bf k}}_{1},\tilde{X}\rangle$, respectively, interacting with the same interaction strength $t({\bf q})$, since

[TABLE]

Thus, by shifting ${\bf k}_{0}$ to ${\bf k}_{0}^{\prime}$, the points discarded outside $k_{c}$ with ${\bf k}_{0}$ are properly considered. And by calculating the electronic structures for every ${\bf k}_{0}$ in the first Brillouin zone, we can get every possible interaction pairs in this system.

Figure 6 shows the quasicrystal bands of QC-TBG calculated with 12-wave model ($k_{c}<3.76/a$ with ${\bf k}={\bf 0}$ removed) and 182-wave model ($k_{c}<18.8/a$). We can see that the band structure of 182-wave model fully includes the spectrum of 12-wave model, while also contains many other band lines. Actually, these extra lines are just replicas of the identical quasicrystal bands of 12-wave model with different origins (Sec. II.2). In other words, the physical properties can be well described by calculating the quasi-band structure with a relatively short $k_{c}$.

The minimum 12-wave models well reproduces the band structure near the 12-wave resonant states, while there are some small errors in the other energies. In Fig. 7, we see that the 182-wave model almost perfectly overlaps with the DOS of very large finite flakes, while 12-wave model slightly under-/overestimates the density of states far from the resonant-state peaks. We confirmed that further increase of $k_{c}$ does not change the DOS profile.

Appendix B Quasicrystal approximants

In Sec. III.1, we compared the quasicrystalline TBG stacked at $30^{\circ}$ (QC-TBG) and its periodic approximant at $29.84^{\circ}$. Actually there exist infintely many periodic TBGs in any finite region in $\theta$, just like rational numbers in the real number axis. As we will see the following, the peak structure in the density of states changes almost continuously in rotating the twist angle $\theta$, and tracing its evolution is useful to get insights on the connection between QC-TBG and the low-angle moiré TBGs, although the computation requires enormous number of atomic bases ($10^{3}-10^{5}$ atoms).

Here we calculated the electronic structures of periodic TBGs with various $\theta$ using a tight-binding model. Figure 8(a) shows the evolution of the DOS in a wide range of energy, and Fig. 8(b) is the magnified plot near the critical states of QC-TBG. The peaks marked with “MG” correspond to the van Hove singularity of monolayer graphene. Similarly, “2-wave-interlayer” represents the singularity originates from the two-wave mixing between the states in different layers Moon and Koshino (2013b), and “2-wave-intralayer” is the mixing between the states in the same layer. Koshino et al. (2015); Yao et al. (2018) We show the width of the band opening (pseudogap) from 2-wave interlayer/intralayer mixing by the green horizontal arrows.

The sharp peaks marked as $\alpha$, $\beta$, $\gamma$ correspond to the singularities coming from the 12-wave mixing of QC-TBG, which were described as the nearly flat bands in the quasi-band picture in Sec. III.1. We can see that the singular peaks rapidly grow as $\theta$ approaches $30^{\circ}$. The DOS of the TBG with $\theta=29.99^{\circ}$, the periodic TBG closest to 30∘ in this calculation, is consistent with that of QC-TBG [Fig. 2(a)]. It should be noted that, however, the wave functions of the approximants do not obey the quasicrystalline long-range structure with 12-fold rotational symmetry, as argued in Fig. 3(b). We also see that the peak-and-dip structure in the valence band is much wider than in the conduction band, and it is consistent with the analytic argument in the 12-wave ring model [Eq. (9)].

The peaks enclosed by the red box in Fig. 8(b) are associated with the resonant states other than $\alpha$, $\beta$, $\gamma$ [i.e., the solutions of Eq. (7) other than $\alpha$, $\beta$, $\gamma$]. In the 12-wave model, we can also show that corresponding states have flat dispersion in quasi-band structure at $\theta\sim 28^{\circ}$, and exhibit singularities in DOS. As $\theta$ approaches $30^{\circ}$, however, the quasi-band becomes dispersive [Fig. 2(b)] and the DOS singularities disappear.