Twists and The Electronic Structure of Graphitic Materials

TL;DR

This paper investigates how twisting layers in graphitic materials affects their electronic properties, revealing complex phenomena like flat bands and resonances, with implications for material design.

Contribution

It introduces a generalized method to analyze electronic structures in twisted layered graphene configurations, applicable to various stacking sequences.

Findings

Presence of sharp resonances at small twist angles

Emergence of flat bands in twisted configurations

Method applicable to complex stacking arrangements

Abstract

We analyze the effect of twists on the electronic structure of configurations of infinite stacks of graphene layers. We focus on three different cases: an infinite stack where each layer is rotated with respect to the previous one by a fixed angle, two pieces of semi-infinite graphite rotated with respect to each other, and finally a single layer of graphene rotated with respect to a graphite surface. In all three cases we find a rich structure, with sharp resonances and flat bands for small twist angles. The method used can be easily generalized to more complex arrangements and stacking sequences.

Figures

Equations

$${\mathbf{G}}_1=\frac{4\pi }{\sqrt{3}L}(1/2,\sqrt{3}/2)\text{ and }{\mathbf{G}}_2=\frac{4\pi }{\sqrt{3}L}(-1,0).$$

$$H_{nn}$$

$${\hat{T}}_{\perp }(0)=\left(\begin{array}{ccc}{g}_1& & g_2\\ g_2& & g_1\end{array}\right),{\hat{T}}_{\perp }(-{\mathbf{G}}_1)=\left(\begin{array}{ccc}{g}_1{e}^{2i\pi /3}& & g_2\\ g_2{e}^{-2i\pi /3}& & g_1{e}^{2i\pi /3}\end{array}\right)={\left[{\hat{T}}_{\perp }(-{\mathbf{G}}_1-{\mathbf{G}}_2)\right]}^{\ast }.$$

$$H_{n,n+1}=\int d^{2}k\sum _{\mathbf{G}}{\Psi }_{n,3}^{†}(\mathbf{k}){\hat{T}}_{\perp }(\mathbf{G}){\Psi }_{n+1,1}(\mathbf{k}-\mathbf{G}-2{\mathbf{G}}_1-{\mathbf{G}}_2)+\text{h.c.},$$

$$H=\sum _{n=-\infty }^{\infty }\left(H_{n,n}+H_{n,n+1}\right),$$

$$H$$

$$H_{k_z}$$

$${\Phi }_{k_z}(\mathbf{k},\mathbf{r})=\sum _{nm}\left[\begin{array}{c}{\varphi }_{nm,k_z}^{A1}(\mathbf{k},\mathbf{r})\\ {\varphi }_{nm,k_z}^{B1}(\mathbf{k},\mathbf{r})\\ {\varphi }_{nm,k_z}^{A2}(\mathbf{k},\mathbf{r})\\ {\varphi }_{nm,k_z}^{B2}(\mathbf{k},\mathbf{r})\\ {\varphi }_{nm,k_z}^{A3}(\mathbf{k},\mathbf{r})\\ {\varphi }_{nm,k_z}^{B3}(\mathbf{k},\mathbf{r})\end{array}\right]\times e^{i\left(\mathbf{k}+n{\mathbf{G}}_1+m{\mathbf{G}}_2\right)\cdot \mathbf{r}},$$

$$S_{\text{eff}}$$

$${\mathcal{L}}_{\text{eff}}(\omega )$$

$$H_{LU}=\int d^{2}k\sum _{\mathbf{G}}{\Psi }_L^{†}(\mathbf{k},\omega ){\hat{T}}_{\perp }(\mathbf{G}){\Psi }_U(\mathbf{k}-\mathbf{G},\omega )+\text{h.c.}$$

$$\rho (\mathbf{k},\omega )=-\frac{1}{\pi }\sum _{\alpha }\mathrm{Im}\left\{{\xi }_{\mathbf{k},\alpha }(\omega +i\delta {)}^{-1}\right\},$$

$$S_{\text{stack}}=\int d\omega d^{2}k\sum _{n\ge 0}\left\{{\Xi }_n^{†}(\mathbf{k},\omega ){\mathcal{G}}_0^{-1}(\mathbf{k},\omega ){\Xi }_n(\mathbf{k},\omega )-\left[{\Xi }_n^{†}(\mathbf{k},\omega ){\hat{T}}_0(\mathbf{k}){\Xi }_{n+1}(\mathbf{k},\omega )+\text{h.c.}\right]\right\},$$

$$\hat{H}(\mathbf{k})=\left(\begin{array}{cccc}0& v_{F}k& -v_{4}k& v_3{k}^{\ast }\\ v_F{k}^{\ast }& 0& {\gamma }_1& -v_{4}k\\ -v_4{k}^{\ast }& {\gamma }_1& 0& v_{F}k\\ v_{3}k& -v_4{k}^{\ast }& v_F{k}^{\ast }& 0\end{array}\right),{\hat{T}}_0(\mathbf{k})=\left(\begin{array}{cccc}{\gamma }_2& 0& 0& 0\\ 0& -{\gamma }_2& 0& 0\\ -v_4{k}^{\ast }& {\gamma }_1& -{\gamma }_2& 0\\ v_{3}k& -v_4{k}^{\ast }& 0& {\gamma }_2\end{array}\right).$$

$$S_{\text{stack}}^{(N)}$$

$${\tilde{\mathcal{G}}}_N^{-1}$$

$${\mathcal{G}}_N^{-1}$$

$${\hat{T}}_N$$

$${\tilde{\mathcal{G}}}_{\infty }^{-1}=\left(\begin{array}{ccc}{\mathcal{G}}_{11}^{-1}& & {\mathcal{G}}_{12}^{-1}\\ {\mathcal{G}}_{12}^{-1,†}& & {\mathcal{G}}_{22}^{-1}\end{array}\right),\text{then}{G}_U^{-1}={\mathcal{G}}_{11}^{-1}-{\mathcal{G}}_{12}^{-1}{\mathcal{G}}_{22}{\mathcal{G}}_{12}^{-1,†}.$$

$$G_L^{-1}(\mathbf{k},\omega )=\left(\begin{array}{ccc}\omega -\Sigma (\mathbf{k},\omega )& & -v_{F}k\\ -v_F{k}^{\ast }& & \omega \end{array}\right),G_U^{-1}(\mathbf{k},\omega )=\left(\begin{array}{ccc}\omega & & -v_{F}k\\ -v_F{k}^{\ast }& & \omega -\Sigma (\mathbf{k},\omega )\end{array}\right),$$

$$\Sigma (\mathbf{k},\omega )=\frac{\omega {\gamma }_1^2}{{\omega }^2-{\left(v_{F}k\right)}^2}\left[1+\sum _{n\ge 1}\frac{{\prod }_{l=1}^n{\left({\gamma }_1{a}_l\right)}^2}{1-\frac{\omega a_n{\gamma }_1^2}{{\omega }^2-{\left(v_{F}k\right)}^2}}\right],$$

$$a_n=\frac{\omega }{{\omega }^2-{\left(v_{F}k\right)}^2-\omega {\gamma }_1^2{a}_{n-1}},a_0=0.$$

$$G_U^{-1}(\mathbf{k},\omega )=G_0^{-1}(\mathbf{k},\omega )=\omega -v_F\sigma \cdot \mathbf{k}.$$

Full text

I Introduction

The discovery of both superconductivity and insulating behavior in twisted graphene bilayers[1, 2, 3], see also[4, 5], has lead to an extraordinary research effort on twisted graphene systems. The twist between neighboring layers leads to the formation of a Moiré superlattice. The size of the unit cell increases as the twist angle is reduced. For sufficiently large lattice units, we can ignore the atomistic nature of the carbon atoms, and describe the low-energy electronic bands by a continuum model [6, 7, 8, 9, 10]. These superlattice bands are very narrow, with a bandwidth of a few meV’s for certain twist angles[11, 12, 8]. This analysis has been generalized to twisted trilayers[13, 14], tetralayers[15, 16], and to some twisted infinite stacks[17].

The problem of twists is related to the occurrence of stacking defects in graphite and at graphite surfaces. Both of these have been extensively studied experimentally, see Refs. [18, 19, 20, 21], and Refs. [22, 23, 24, 25, 26]. The electronic properties of stacking defects in graphite have also been studied theoretically[27, 28]. It is worth noting that single graphene layers on graphite can be manipulated with atomic scale probes[29], and that unusual superconducting like features in bulk graphite have been ascribed to stacking misalignment[30, 31].

In the following, we study three representative cases of three dimensional systems with a twist. A sketch of these structures is shown in Fig. 1. As described below, the analysis presented here can easily be extended to other intermediate situations. The first case is a stack of a graphene layers where each layer is rotated by a constant amount with respect to the previous one. This arrangement can be achieved by applying a global twist to three-dimensional graphite, provided that neighboring layers slide and lose registry at the atomic scale in the same way. This structure is a natural generalization to three dimensions of a twisted graphene bilayer. The two other situations both describe rotations involving stacks of bulk graphite. We use the Bernal stacking sequence, $ABAB,\cdots$ for graphite, although extension to other sequences is straightforward. In the second case we consider two semi-infinite graphite crystals twisted with respect to each other, and the final one is a single graphene layer twisted with respect to a graphite surface. The twist between layers leads to narrow bands in many cases, which will enhance the effect of the electron-electron interactions and may lead to broken-symmetry phases. The study of such features lies beyond the scope of this work.

We now turn to a detailed description of the three models. In the continuum model, we take the parametrization of the interlayer hoppings between the twisted layers from Ref. [32], where the $AA$ and $AB$ hoppings are described by the two parameters $\{g_{1},g_{2}\}=\{79.7,97.5\}\,\mathrm{meV}$. We will use a twist angle $\theta=1.08^{\circ}$, which is the first magic angle for the parametrization given above.

II Infinite Stack of Twisted Layers

In twisted bilayer graphene we find a Moiré superlattice, which does not require perfect alignment. A a commensurate superlattice at the atomic scale can only be defined for certain angles where reciprocal vectors of the two layers coincide. The continuum approximation developed in Refs. [6, 8] describe an approximate Moiré superlattice for any twist angle, with a lattice unit $L=d/(2\sin(\theta/2)$, where $d$ is the lattice unit of graphene and $\theta$ is the twist angle. This continuum description has been shown to provide an excellent approximation to the electronic bands at commensurate angles at the atomic scale, and is expected to describe the properties of the bilayer at arbitrary angles.

The general analysis of the electronic properties for more than two misaligned layers is more challenging, as the superposition of different Moiré patterns leads to complex structures with no obvious periodicity. A general analysis for the case of three rotated layers is given in Ref. [13]. The complexity of that method increases exponentially with the number of layers. The situation simplifies when the twist angle between successive layers is of equal magnitude and opposite sign between neighboring pairs of layers [17]. In the case of constant twist angle, even in the simple case of a trilayer, the Moiré structures defined by each pair of layers, which are of equal size, do not coincide, and a full description requires the consideration of many additional reciprocal vectors. As the angle is reduced, the misalignment between the two Moiré patterns decreases. By assuming that the two interlayer rotations define exactly the same Moiré pattern, a continuum approximation to the bands of the twisted trilayer can still be obtained[14], where the required number of reciprocal lattice vectors is of the same order of magnitude as that for a twisted bilayer.

We describe the infinite rotated stack as a succession of trilayers embedded in an environment to be defined self consistently, in an analogous manner as the Coherent Phase Approximation[33]. We use three sets of in-plane reciprocal lattice momenta, one for each layer, labeled by a two-dimensional vector, $\boldsymbol{\mathrm{k}}_{\parallel}^{i},i=1,2,3$. These are then repeated periodically in the $z$ direction. In a similar way as for the bilayer [8], each set of momenta is defined in a triangular lattice with reciprocal lattice spacing $|\boldsymbol{\mathrm{G}}|=(4\pi)/(\sqrt{3}L)$. We introduce a wavevector $k_{z}$ along the $z$ direction. The relative phase between successive trilayers is $k_{z}\bar{d}_{z}$, where $\bar{d}_{z}$ is three times the interlayer distance, $c$. This allows us to classify the eigenstates as function of $\boldsymbol{\mathrm{k}}_{\parallel}$ and $k_{z}$, where $\boldsymbol{\mathrm{k}}_{\parallel}$ is the same for the three layers. We consider only couplings between nearest-neighbor layers. Couplings between second-nearest-neighbor layers, given in the Bernal stacking by the hopping parameter $\gamma_{2}$, (see below) give rise to couplings between $AA$ sites modulated over the Moiré supercell. We have not included these terms, since they are only a small perturbation to the hamiltonian used here.

The assumption of a well defined Moiré structure is only valid for small angles. As the angle increases, the misalignment between layers at a given distance also increases. This misalignment tends to decouple the layers, and thus the model used here overestimates the interactions between distant layers, so that our results may overestimate the dispersion along the $k_{z}$ direction. Our framework assumes that each block of three consecutive twisted layers defines an approximate Moiré pattern. The same approach applies for different stacking arrangements: $AAA$, $ABA$ or $ABC$. If the twist angle is small enough, the approximate Moiré pattern identified by each block varies slowly along the stack. As we assume that we need to consider only nearest-neighbor interlayer couplings, we can locally describe the stack as a single Moiré pattern, disregarding the type of stacking arrangement.

Typical bands, as function of $\boldsymbol{\mathrm{k}}_{\parallel}$ and $k_{z}$ are shown in Fig. 2a,b, along with the corresponding density of states.

The results show the coexistence of very narrow and wider bands over a broad range of energies. The bands at $k_{z}=0$ and $k_{z}=\pi$ are equivalent under the simultaneous transformations $\boldsymbol{\mathrm{k}}_{\parallel}\leftrightarrow-\boldsymbol{\mathrm{k}}_{\parallel}$ and $E\leftrightarrow-E$. The narrow bands give rise to sharp peaks in the density of states. The overlapping of narrow and broad bands is reminiscent of a twisted trilayer[14, 34]. We find two sharp minima in the width of the lowest energy bands for $\theta\simeq 1^{\circ}$ (see Fig. 2d)), suggesting the possibility of defining two “magic angles”. The charge density distribution within the unit cell as function of wavevector is shown in [34]. The charge densities for the narrowest band near the neutrality point are very similar for the $K,K^{\prime}$ and $\Gamma$ points. This is rather different from the significant momentum dependence of the charge distributions in twisted graphene bilayers[14, 36].

III Twisted Stacking Fault in Bulk Graphite

We now consider two semi-infinite graphite stacks rotated with respect to each other by an angle $\theta$. The misalignment at the interface between the two stacks induces a Moiré pattern, in an analogous case to a twisted graphene bilayer. The standard calculation of the band structure of a twisted graphene layer[6, 8] requires the knowledge of the Green’s function of the states with momenta mixed by the Moiré superstructure. For the twisted bilayer this Green’s function can be obtained from the inversion of the hamiltonian of each individual layer, $G_{i}(\boldsymbol{\mathrm{k}}_{\parallel},\omega)=(\omega{\cal I}_{\sigma_{i}}-{\cal H}_{i,\boldsymbol{\mathrm{k}}_{\parallel}})^{-1}$, where $\omega$ is the frequency, $i=1,2$ is a layer index, ${\cal I}_{\sigma}$ is a $2\times 2$ matrix defined in the sublattice space of layer $i$, and ${\cal H}_{i,\boldsymbol{\mathrm{k}}_{\parallel}}$ is a $2\times 2$ matrix which defines the hamiltonian in layer $i$. In the case of a defect embedded in a three dimensional stack, the Green’s function at each side of the defect can be expressed in terms iof the slef energy, $G_{i}(\boldsymbol{\mathrm{k}}_{\parallel},\omega)=[\omega{\cal I}_{\sigma_{i}}-\Sigma_{i}(\omega,\boldsymbol{\mathrm{k}}_{\parallel})]^{-1}$, where $\Sigma_{i}(\omega,\boldsymbol{\mathrm{k}}_{\parallel})$ has to be chosen so that $G_{i}(\boldsymbol{\mathrm{k}}_{\parallel},\omega)$ gives the Green’s function at the surface of a semi-infinite stack. These self energies can be obtained by recursive methods[37, 38, 34]. The required interlayer hopping parameters that describe the bands of bulk graphite are well known in terms of the SWM parameters, and have been extensively studied [39, 40, 41, 42].

Results for the density of states of the two layers at the stacking defect are shown in Fig. 3. The nature of the bands, or resonances, at very low energies depends sensitively on the value of the hopping between second-nearest-neighbors layers, $\gamma_{2}$, and results for the graphite value $\gamma_{2}$ and $\gamma_{2}=0$ are shown. This parameter is the smallest and the least precisely determined of the set which describes the electronic structure of bulk graphite.

The assumption $\gamma_{2}=0$ in modelling graphite leads to a nodal line of Fermi points instead of the observed two-dimensional Fermi surface. A finite value of $\gamma_{2}$ is also needed in order to explain the finite bandwidth of the $n=0$ Landau level of graphite[43, 44, 45]. Experiments in bulk graphite are consistent with a value $\gamma_{2}\approx-0.02\,\text{eV}$, see, for instance[46, 47]. On the other hand, recent measurements on 4-8 Bernal-stacked multilayers [48, 49, 50, 51] suggest that the second-nearest-neighbor interlayer hopping is much smaller, $\gamma_{2}\sim 0$.

The broadening of the features shown in Fig. 3 is due to the overlap of states localized at the defect with continuum states in bulk graphite. A simple analysis of the bulk bands projected onto a two-dimensional surface shows that the edges of the continuum lie approximately at $\epsilon(\boldsymbol{\mathrm{k}}_{\parallel})=\pm\gamma_{2}\pm(3\gamma_{0}|\boldsymbol{\mathrm{k}}_{\parallel}|d)|^{2})/(8\gamma_{1})\pm(\sqrt{3}\gamma_{3}|\boldsymbol{\mathrm{k}}_{\parallel}|d)/2$, where we use the conventional notation for the band structure of graphite [41], $\gamma_{4}$ and $\gamma_{5}$ have been neglected, and we assume that $\gamma_{2}\ll\gamma_{0},\gamma_{1},\gamma_{3}$. The above expression for $\epsilon(\boldsymbol{\mathrm{k}}_{\parallel})$ shows that states at $\epsilon=0$ do not overlap with the bulk continuum except for $|\boldsymbol{\mathrm{k}}_{\parallel}|=0$. For $\gamma_{2}\neq 0$, and neglecting the effect of $\gamma_{3}$, the states at $\epsilon=0$ overlap with the continuum for $|\boldsymbol{\mathrm{k}}_{\parallel}|=k_{0}\lesssim\sqrt{8\gamma_{1}\gamma_{2}}/(3\gamma_{0}d)$. The distance between the $\Gamma$ and $K$ points in the Moiré Brillouin Zone is $k_{\Gamma-K}=[3\sin(\theta/2)]/(2\pi d)$. For $\theta\lesssim 0.4^{\circ}$, we find $k_{\Gamma-K}\lesssim k_{0}$, and the entire Brillouin Zone overlaps with the continuum[52]. Note that the role of $\gamma_{2}$ in the case of the stacking fault considered here is very different from the role of $\gamma_{2}$ in the continuously twisted stack analyzed previously, as, in the present case, the series of second nearest neighbors at both sides of the defect are in registry.

Information about the real-space origin of the spectral features of the twisted stack shown in Fig. 1(b) can be obtained from Fig. 4. Here we show the local DOS on the boundary surface, computed in the $AA$ (red line), $AB$ (blue line) and $BA$ (brown line) regions of the Moiré unit cell. As can clearly be seen, the DOS in the $AA$ region strongly exceeds those in the $AB/BA$ regions at the Van Hove singularity, which is in good agreement with the scanning-tunneling-microscopy maps reported in twisted graphene bilayer [53, 54, 55].

IV Twisted Layer on a Graphite Surface

The analysis presented above can also be applied when a graphene layer is rotated with respect to the surface of a semi-infinite stack of graphene layers. In this case, the problem can be reduced to a twisted bilayer where the self-energy term influences only one layer. Results are shown in Fig. 5.

The results in Fig. 5 are rather similar to those in Fig. 3, but the features of the flat bands at low energies are sharper. The reduction in the value of $\gamma_{2}$ observed in not too large multilayers[48, 49, 50, 51], might be related to an increase in the interlayer distance. If that is the case, a similar effect can take place at surfaces, enhancing the sharpness of the low energy resonances.

V Lattice Relaxation

Up to now we have assumed that there is no reconstruction of the lattice at the interface; it has been shown that this is an important effect for bilayers, see, e.g., Refs. [56, 57, 58], and can have important consequences for the electronic structure as well. We focus on in plane relaxation, which is the largest contribution to the total relaxation [59, 60, 61, 56, 57], especially for the (semi-)infinite structures studied here. We show results for a graphene layer on a graphite surface, as, in this case, the top layer is weakly coupled to the rest of the system. We optimize the positions in the atoms in the twisted layer and in the topmost layer at the graphite surface. The intra-layer potential used is the AIREBO-M potential [62], and the Kolomogorov-Crespi inter-layer potential [63]. We approximate the graphite as three layers of graphene, with the bottom two layers fixed, with one twisted graphene layer on top. Results are shown in Fig. 6. Additional results for the relaxation of a twisted bilayer on a graphite surface are shown in [34]. The plots in Fig. 6 show a significant relaxation of the topmost layer, in line with calculations for a twisted graphene bilayer[56, 57]. The first layer of the graphite stack is more weakly perturbed. The maximal displacement of atoms in the graphene layer is $0.36\,\text{\AA}$, and in the top layer of the graphite $0.24\,\text{\AA}$. Comparing Fig. 6b to results for bilayer graphene in [56, 57], we see that the deformation is less than that for a bilayer, but not by much. Since this would largely impact the hopping parameters at the interface, we expect that the electronic structure could be quite sensitive to this change.

VI Conclusions

We have analyzed the changes of the electronic properties due to small angle rotations between neighboring layers in three dimensional stacks of graphene layers.

We first looked at a uniformly rotated stack, which corresponds to the application of a shear strain to bulk graphite. In this case we find both flat and dispersive bands, which overlap at low energies. We find, nevertheless, sharp resonances near the neutrality point at well defined twist angles.

In the case where one semi-infinite stack is rotated with respect to another, i.e., a stacking fault in bulk graphite, sharp features in the density of states appear in the defect region. Localized states with low dispersion exist within bulk gaps. The regions in energy and momenta where these gaps can be defined depend sensitively on the hopping between second-nearest-neighbor layers.

A single layer of graphene on a graphite surface shows similar low energy resonances and states within bulk gaps as a rotated stacking fault, although the resonances become sharper due to the disappearance of one half of the bulk ”reservoir”, and the coupling to this continuum becomes weaker.

The momentum dependence of the charge density distributions within the Moiré unit cell varies in the three cases discussed here. In twisted graphite the charge is peaked at the $AA$ regions, independently of momentum. This is different from the case of twisted bilayer graphene, and suggests the existence of simple Wannier functions. The momentum dependence is significant for a graphene layer on a graphite surface. This dependence implies that the shape of the bands will be modified when the system is doped away from the neutrality point[36], and an additional interaction, electron assisted hopping, can be defined. The strength of this interaction will be discussed elsewhere.

The results presented here will be sensitive to relaxation of the atomic positions at the interface, which we have also calculated. It is smaller than in bilayers, but not by a significant amount.

The problems discussed here are just one of many one could envisage in stacks of graphene layers. The method used for the last two problems relies on the substitution of the periodic stack of layers near the defect by an effective self energy. This method is very general, and it can be extended to other combinations of Bernal, or rombohedral, graphite with rotated stacking faults and/or rotated layers at surfaces. Inhomogeneous electrostatic potentials induced by external gates can also be included, see [34].

Acknowledgements.

Acknowledgements. We would like to thank Pablo San José for useful conversations. This work was supported by funding from the European Commission under the Graphene Flagship, contract CNECTICT-604391.

S7 Infinite Stack of Twisted Layers. Continuum Model

We consider the continuum model for an infinite stack of layers of graphene with a constant small twist $\theta$ between each layer and the consecutive one, as depicted in the Fig. 1(a) of the main text.

As mentioned in the main text, a commensurate Moiré pattern cannot be defined in rotated graphene multilayers. However, if the angle between successive layers is small and constant, a relatively small number of layers still generates an aprroximate periodicity similar to that found in twisted bilayer for a general twist angle where the layers are not commensurate. An example is shown in Fig. S1, where we show the case of the approximate Moiré patterns formed by rotating two (left), four (center) and six (right) stacked layers. If the out-of-plane hopping between layers is restricted to nearest neighbor ones, one can attempt to describe the infinite structure in terms of the Moiré pattern identified by each pair of consecutive layers. This argument can be made rigorous by defining a succession of trilayers embedded in an environment to be defined self-consistently, in manner analogous to that applied the Coherent Phase Approximation[1]. However, we will not focus on this issue in the following and we just assume that the infinite stack identifies a single Moiré pattern, with characteristic length $L=\frac{d}{2\sin\theta/2}\simeq d/\theta$, where $d=2.46\,\text{\AA}$ is the lattice constant of monolayer graphene. This allows to define a reciprocal lattice generated by the two basis vectors

[TABLE]

The corresponding Brillouin zone (BZ) is shown in Fig. S2, left, where we have $\mathbf{K}\equiv 4\pi(0,1)/(3L)$ and $\mathbf{M}\equiv\pi(1,\sqrt{3})/\left(\sqrt{3}L\right)$. It is interesting to note that the model has a non-symmorphic symmetry, as a translation by one layer, followed by a twist, leaves the Hamiltonian invariant. The calculations are done using a unit cell which contains three layers, see below. This implies that a new, smaller, Brillouin Zone can be defined, see Fig. S2, right. We will follow the convention used in the study of twisted bilayers, and use the reciprocal lattice vectors $\mathbf{G}_{1},\mathbf{G}_{2}$ of Eq. (S1). When studying its spectrum, we will use the high-symmetry points in the small triple-layer Brillouin Zone, shown on the right in Fig. S2. It is worth noting that the same situation occurs in rhombohedral graphite which has a $ABCABC\cdots$ stacking.

The basic building block for the model presented here is twisted trilayer graphene, and we built an infinite system by stacking these on top of each other. The top layer of each block is coupled to the bottom layer of the next one, so that neighboring layers of consecutive blocks form a twisted bilayer with the same twist angle. Each block can be modeled within the continuum framework first introduced in [2, 3, 4] for the twisted bilayer and recently extended to the twisted trilayer in Ref. [5]. Within this continuum model, the two valleys of each graphene layer at $K$ and $K^{\prime}$ can be treated separately. This approximation is still valid for the case of bulk three-dimensional systems considered here, where inter-valley effects occur only at very high energies [6, 7, 8], and hence they do no not affect the band structure close to the Dirac points. Consequently, in the following we will focus only on the $K$ valley of each layer. Using the notation of [2, 4], we write the Hamiltonian of the $n$-th block as

[TABLE]

where $n,l$ are the block and layer indices, respectively, $v_{F}=d\gamma_{0}\sqrt{3}/2$ is the Fermi velocity of monolayer graphene, and $\gamma_{0}=2.47\,\text{eV}$ [9, 10] is the nearest-neighbor hopping amplitude. The wavefunction $\Psi^{T}_{n,l}\equiv(\psi^{A}_{n,l},\psi^{B}_{n,l})$, with $A,B$ the sub-lattice indices, $\boldsymbol{\mathrm{\sigma}}\equiv(\sigma_{x},\sigma_{y})$ are the Pauli matrices, $\mathbf{G}=n\mathbf{G}_{1}+m\mathbf{G}_{2}$, $n,m$ integers. Finally $\hat{T}_{\perp}(\mathbf{G})$ is the $2\times 2$ matrix which parametrizes the hopping between consecutive layers. To a first approximation we consider only the contribution of the first three components of $\hat{T}(\mathbf{G})$: $\hat{T}_{\perp}(0),\hat{T}_{\perp}(-\mathbf{G}_{1}),\hat{T}_{\perp}(-\mathbf{G}_{1}-\mathbf{G}_{2})$, with the corresponding matrix elements given by (see [2, 10])

[TABLE]

We shall employ the values recently obtained in [10]: $g_{1}=0.0797$eV, $g_{2}=0.0975$eV. The difference between $g_{1}$ and $g_{2}$, not considered in the former models of twisted bilayer graphene, introduces energy gaps between the lowest bands and the excited bands, in qualitative agreement with the experimental results [11, 12, 13].

Notice that, in the absence of the interlayer interaction, the Hamiltonian $H_{nn}$ gives rise to Dirac cones at the point $\pm\mathbf{K},\Gamma$ of the BZ. It is also worth noting that in Eq. (S7) we are neglecting the small rotation of the pseudo-spin $\boldsymbol{\mathrm{\sigma}}$, which arises due to the twist between consecutive layers. This approximation has been shown to be valid for the case of trilayer graphene in Ref. [5]; we assume that it still holds in the case of infinite layers, as long as we can identify a Moiré Broillouin zone, i.e., for small twist angles.

The Hamiltonian coupling consecutive blocks can been written as

[TABLE]

and describes the hopping between the third layer of the $n$-th block and the first layer of the $(n+1)$-th block.

We note that the Dirac points of the $n$-th block are reciprocally connected to the ones of the $n+1$-th block by the reciprocal lattice vector: $2\mathbf{G}_{1}+\mathbf{G}_{2}$, which is why the same vector appears in the argument of $\Psi_{n+1,1}$ in the rhs of the Eq. (S4). Such a situation occurs only if the building blocks of the stack are arranged in multiples of three layers, which therefore justifies our choice of partitioning the stack in such blocks.

The full Hamiltonian is then

[TABLE]

and can be easily diagonalized by exploiting its translational invariance along the $z$ axis perpendicular to the layers by Fourier transformation. In term of the momentum $k_{z}$ we find

[TABLE]

where $c=3.35$Å is the interlayer spacing and $\Psi_{k_{z},l}({\bf k})\equiv\sqrt{\frac{3c}{2\pi}}\sum_{n}\Psi_{n,l}({\bf k})e^{-ik_{z}3cn}$.

The interlayer interaction of the Hamiltonian (S7) then couples states of layer 1 (2) with momentum ${\bf k}$ to states of layer 2 (3) with momentum ${\bf k}$, ${\bf k}+\mathbf{G}_{1}$ and ${\bf k}+\mathbf{G}_{1}+\mathbf{G}_{2}$, and states of layer 3 with momentum ${\bf k}$ to states of layer 1 with momentum ${\bf k}-2\mathbf{G}_{1}-\mathbf{G}_{2}$, ${\bf k}-\mathbf{G}_{1}-\mathbf{G}_{2}$ and ${\bf k}-\mathbf{G}_{1}$.

The Hamiltonian (S7) can be diagonalized in the reciprocal space by looking for the eigenfunctions as a superposition of Bloch waves,

[TABLE]

where ${\bf k}$ belongs to the Moiré BZ shown in Fig. S2. We perform the numerical diagonalization by truncating the Bloch waves expansion at momenta large enough to achieve the convergence of the lowest energy bands. The results shown in the Fig. 2 of the main text and in the following have been obtained by considering a set of 169 Bloch states.

In Fig. S3 we show the bands and the density of the states (DOS) for a few additional twist angles, $\theta=0.82^{\circ},0.55^{\circ}$. The red line superimposed to the DOS represents the case of three uncoupled layers.

In Fig. S4 we show the charge-density distribution for the three lowest energy bands obtained for the angle $\theta=1.08^{\circ}$, at $k_{z}=0$ and at the points $\Gamma$ (left), $K$ (central), $K^{\prime}$ (right) of the small Moiré BZ. The middle row, which corresponds to the flat band, shows that the charge density is zero in the $AA$ region (see the schematic representation of the layer alignment in the unit cell in the bottom panel). The central band shows similar charge distributions for the three high symmetry points, while the neighboring bands show inequivalent distributions for each of the three points. This situation is different from that found in a twisted bilayer, where the $K$ and $K^{\prime}$ points are equivalent, and different from the $\Gamma$ point[35]. We do not explore here how these symmetries influence the Wannier functions which describe these bands[10].

S8 Twisted Stacking Fault in Bulk Graphite

In contrast to the previous section, where we considered an infinite stack of graphene layers with a constant twist between a layer and the consecutive one, here we study the case of one single twist between two semi-infinite stacks placed one above the other, as shown in the Fig. 1(b) of the main text. Every stack is assumed to be of Bernal type. This situation is indeed particularly meaningful for describing defects in graphite.

We neglect hopping between next nearest neighbor layers in different layers across the stacking defect, so that the coupling between the two stacks is restricted to the layers at the interface. Our approach consists in studying the twisted bilayer represented by the bottom layer of the upper ($U$) stack and the top layer of the lower ($L$) stack, once the interior (core) layers have been integrated out. We describe this bilayer within the path integral formulation, by means of the following effective action

[TABLE]

where $\Psi^{T}_{i}=(\psi^{A}_{i},\psi^{B}_{i})$ is the spinor for the layer placed at the surface of the stack $i$, $G_{i}$ is the corresponding effective Green’s function, which encodes the information concerning the spectrum of graphite, and $H_{LU}$ is the interlayer Hamiltonian of a twisted bilayer[2, 4],

[TABLE]

We investigate the spectrum of the system by computing the DOS per unit of momentum $\rho({\bf k},\omega)$, which is defined by the eigenvalues $\xi_{{\bf k},\alpha}(\omega)$ of $\mathcal{L}_{\text{eff}}(\omega)$ according to

[TABLE]

where $\alpha$ is the band index and $\delta>0$ is a finite spectral broadening.

In order to describe the procedure for computing the Green’s function $G_{i}$, in the following we focus on the upper stack $i=U$, without loss of generality. We consider all its sub-stacks $\left(BA\right)_{0},\left(BA\right)_{1}\dots$, where [math] denotes the bottom block. The action of the stack is

[TABLE]

where $\Xi^{T}_{n}=(\psi_{A1,n},\psi_{B1,n},\psi_{A2,n},\psi_{B2,n})$ is the spinor for the bilayer in the $n$-th sub-stack, $\mathcal{G}^{-1}_{0}({\bf k},\omega)=\omega\mathcal{I}-\hat{H}({\bf k})$ is the inverse of the non-interacting Green’s function. The operators $\hat{H}({\bf k})$, $\hat{T}_{0}({\bf k})$ are the $4\times 4$ matrices of the couplings in each sub-stack and between consecutive sub-stacks, respectively [7, 8, 15]:

[TABLE]

Here $k\equiv(k_{x}-K^{\prime}_{x})-i(k_{y}-K^{\prime}_{y})$, the Dirac point lying at $\mathbf{K^{\prime}}$, $\gamma_{1}$ is the coupling between orbitals in atoms that are nearest neighbors in consecutive layers, $v_{3,4}\equiv v_{F}\gamma_{3,4}/\gamma_{0}$, where $\gamma_{3},\gamma_{4}$ are the hopping amplitudes between orbitals at next nearest neighbors in consecutive layers, and $\gamma_{2}$ is the coupling between orbitals at next nearest neighbors layers. Notice that we are neglecting the parameter $\Delta$ and we are assuming $\gamma_{5}=\gamma_{2}$. The values of the parameters $\gamma$ in graphite have been intensively studied (see[6, 16, 17, 18, 19, 20]). In the following we refer to those reported in [21]: $\gamma_{0}=3.16$eV, $\gamma_{1}=0.39$eV, $\gamma_{2}=-0.02$eV, $\gamma_{3}=0.315$eV, $\gamma_{4}=0.044$eV.

Integrating the sub-stacks $1,\dots N$ out of the action (S12), the effective action of the remaining sub-stacks can be written as:

[TABLE]

where $\tilde{\mathcal{G}}_{N}$ is the effective Green’s function of the bottom sub-stack, $\mathcal{G}_{N}$ that of the $(N+1)$-th sub-stack and $\hat{T}_{N}$ is the effective coupling between them. It can be easily shown that $\tilde{\mathcal{G}}_{N},\mathcal{G}_{N}$ and $\hat{T}_{N}$ are obtained from the $(N-1)$-th step according to the recursive relations:

[TABLE]

with $N>0$ and $\tilde{\mathcal{G}}_{0}=\mathcal{G}_{0}$. Because the inter-layer hopping amplitudes $\gamma_{1}\dots\gamma_{4}$ are smaller than the intra-layer one $\gamma_{0}$, then $\hat{T}_{N}\to 0$ as $N\to\infty$. Thus $\tilde{\mathcal{G}}_{N}$ converges to a fixed point $\tilde{\mathcal{G}}_{\infty}$, which describes the dynamics of the bottom sub-stack accounting for all the interactions with the upper ones. This method can of course also be used to study a finite stack with $2N$ layers.

The $2\times 2$ inverse Green’s function $G^{-1}_{U}$ introduced in Eq. (S9) is the projection of $\tilde{\mathcal{G}}^{-1}_{\infty}$ on the first $2\times 2$ block, meaning that, if

[TABLE]

The procedure for obtaining $G^{-1}_{L}$ is equivalent to that described above. The only changes are the swap of the indices $1\leftrightarrow 2$ and $3\leftrightarrow 4$ and $k\leftrightarrow k^{*}$ in the matrices of the Eq. (S13), as the stack is now arranged as $(AB)_{0},(AB)_{1},\dots$.

The spectra in Fig. 3 of the main text show the DOS, $\rho({\bf k},\omega)$, given by the Eq. (S11) and obtained by truncating the iteration (S15) at $N=50$, where convergence can be considered reached. We use the spectral broadening $\delta=10^{-3}\gamma_{0}$.

Finally we note that in the simplified case where $\gamma_{2}=\gamma_{3}=\gamma_{4}=0$ while $\gamma_{1}$ is finite, $G^{-1}_{i}$ takes the following analytic form:

[TABLE]

where

[TABLE]

and the coefficients $a_{n}$ satisfy the recursive equation:

[TABLE]

We show the DOS resulting from this analytic Green’s functions (S16) for the twist angle $\theta=1.08^{\circ}$ in Fig. S5. The isolated flat bands appearing at zero energy are due to the absence of the $\gamma_{2},\gamma_{3},\gamma_{4}$ couplings.

S9 Twisted Layer on a Graphite Surface

The formalism developed in the previous section can be easily adapted to describe the situation in which a single layer of graphene is rotated with respect to bulk graphite. We only need to compute $G^{-1}_{L}$ as above, and can use the graphene result for the other:

[TABLE]

This is then substituted into Eq. (S9). The resulting spectra are shown in the Fig. 4 of the main text. In Fig. S6 we show the local DOS at the surface between graphite and the monolayer, integrated over a narrow energy window of $40\,\text{meV}$ centered at the Fermi level, at the points $\Gamma$ (left) and $K$ (right) of the Brillouin zone.

These distributions resemble the case for twisted bilayer graphene[35]. They suggest that the symmetries present in twisted bilayer graphene are at least approximate symmetries in this system.

We can also apply a voltage at the interface between graphite and the twisted layer, by putting a gate above the top layer. We include the effect of the gate by means of introducing two different chemical potentials, $\mu_{L}\neq\mu_{U}$, in $G^{-1}_{L}$ and $G^{-1}_{U}$, respectively. Figure S7 shows the spectra obtained for $\{\mu_{L},\mu_{U}\}=\{5,10\},\{5,20\},\{5,30\}$meV.

The changes due to the gate potential are rather small, although noticeable. The small effect of the gate is probably due to level repulsion between the low energy band and the continuum coming from the semi-infinite stack.

S10 Lattice Relaxation

The idea of the lattice relaxation is discussed in great detail in Refs. [22, 23], and the results presented here are a based on the work there. We use classical potential models as implemented in the LAMMPS code [24] to relax the atomic positions. The potential models used here are the AIREBO-M[25] for inlayer interactions, overlaid by the Kolmogorov-Crespi potential [26] for the long-range interactions between the layers.

We have modelled the atomic positions for the last two problems only using this approach, since the first one is not easily attacked by such a simulations due to the lack of periodicity. We show the results for the single graphene layer in Fig. 5 of the main text. We model the graphite as an $ABA$ stack with the lowest two $AB$ layers completely fixed, at the graphite spacing. The top $A$ layer and the rotated graphene layer are left completely free. We have performed a similar calculation for graphite on graphite, where each stack is modelled as a movable $A$ layer on a fixed $B$ layer, and we assume all layers are spaced at the graphite layer spacing, which fortunately is very close to the equilibrium graphene-graphene spacing. As we can see in Fig. S8 the surface reconstruction is very similar to that for graphene on graphite: The alignment between the two movable layers seems very similar. This means that the effect of relaxation should be rather universal as well, since the hopping parameters at the twist interfaces are expected to be similar.