Coexistence of ultraheavy and ultrarelativistic Dirac quasiparticles in sandwiched trilayer graphene

TL;DR

This paper demonstrates that twisted sandwiched trilayer graphene can stably host coexisting ultraheavy flat bands and ultrarelativistic Dirac quasiparticles near the Fermi energy, controllable by external fields, offering a new platform for exploring strongly correlated and topological phases.

Contribution

It introduces a stable twisted sandwiched graphene structure with coexisting flat and Dirac bands, controllable via external fields, advancing the understanding of strongly correlated quantum materials.

Findings

Stable flat bands coexist with Dirac cones at 1.5° twist.

Atomistic calculations predict relaxation to high-symmetry stacking.

External fields can tune the energy offset between bands.

Abstract

Electrons in quantum materials exhibiting coexistence of dispersionless (flat) bands piercing dispersive (steep) bands can give rise to strongly correlated phenomena, and are associated with unconventional superconductivity. It is known that in twisted trilayer graphene steep Dirac cones can coexist with band flattening, but the phenomenon is not stable under layer misalignments. Here we show that such a twisted sandwiched graphene (TSWG) -- a three-layer van der Waals heterostructure with a twisted middle layer -- can have very stable flat bands coexisting with Dirac cones near the Fermi energy when twisted to 1.5 degrees. These flat bands require a specific high-symmetry stacking order, and our atomistic calculations predict that TSWG always relaxes to it. Additionally, with external fields, we can control the relative energy offset between the Dirac cone vertex and the flat bands. Our work establishes twisted sandwiched graphene as a new platform for research into strongly interacting phases, and topological transport beyond Dirac and Weyl semimetals.

Equations

$$v(\theta )=v_0\frac{1-6{\alpha }^2(\theta )}{1+6{\alpha }^2(\theta )}.$$

$${\alpha }_{\ast }=\frac{1}{\sqrt{6}}\approx 0.408,{\theta }_{\ast }=\frac{w_{AB}}{{\alpha }_{\ast }{k}_D{v}_0}.$$

$$E_{\mathrm{intra}}[{\mathbf{u}}_i(\mathbf{b})]$$

$$\displaystyle=\int_{\Gamma}d\bm{b}\,\frac{1}{2}\Big{[}G(\partial_{x}u_{i,x}+\partial_{y}u_{i,y})^{2}$$

$$\displaystyle+K((\partial_{x}u_{i,x}-\partial_{y}u_{i,y})^{2}+(\partial_{x}u_{i,y}+\partial_{y}u_{i,x})^{2})\Big{]},$$

$$V_{\mathrm{GSFE}}(v,w)$$

$$+c_2(\cos (v+2w)+\cos (v-w)+\cos (2v+w))$$

$$+c_3(\cos 2v+\cos 2w+\cos (2v+2w)),$$

$$E_{\mathrm{inter}}[{\mathbf{u}}_1,{\mathbf{u}}_2]=2{\int }_{\mathrm{\Gamma }}d\mathbf{b}{V}_{\mathrm{GSFE}}(\mathbf{b}+{\mathbf{u}}_1(\mathbf{b})-{\mathbf{u}}_2(\mathbf{b})),$$

$$E_{\mathrm{tot}}(\mathbf{u}(\mathbf{r}))=\sum _{i=1}^3{E}_{\mathrm{intra}}({\mathbf{u}}_i(\mathbf{b}))+E_{\mathrm{inter}}({\mathbf{u}}_1(\mathbf{b}),{\mathbf{u}}_2(\mathbf{b})).$$

$$\mathbf{b}=(E_1^{-1}{E}_2-1)\mathbf{r},$$

$$\mathcal{H}=\left(\begin{array}{ccc}-i\hslash v_0({\mathbf{\sigma }}_{-\theta /2}\cdot \nabla )& T(\mathbf{r})& 0\\ T^{†}(\mathbf{r})& -i\hslash v_0({\mathbf{\sigma }}_{+\theta /2}\cdot \nabla )& T^{†}(\mathbf{r})\\ 0& T(\mathbf{r})& -i\hslash v_0({\mathbf{\sigma }}_{-\theta /2}\cdot \nabla )\end{array}\right),$$

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

$$T_n$$

$${\hat{\Omega }}_{\varphi }$$

Full text

{affiliations}

Department of Physics, Harvard University, Cambridge, MA 02138, USA

Harvard J.A. Paulson School of Engineering and Applied Sciences, Cambridge, MA 02138, USA

Graphene, an atomically thin crystal of carbon, provides an experimentally favorable platform for two dimensional (2D) Dirac physics as it exhibits ultrarelativistic Dirac cones in its band structure, described with massless quasiparticles when weak spin-orbit coupling is neglected.[1] Bilayers of graphene in the energetically favorable Bernal ($AB$) stacking have quadratic dispersion and quasiparticles with well-defined effective mass.[1] Twisted bilayer graphene (TBG) – two rotationally mismatched graphene layers – can be fabricated at the so-called magic angle near $1.1^{\circ}$, where it hosts ultraheavy fermions with remarkably flat, almost dispersionless electronic bands[2, 3, 4, 5] of a topological origin.[5, 6, 7, 8] The twist angle serves as a precise control of the interlayer coupling between the graphene monolayers, revealing the band flattening phenomena as an ultimate manifestation of hybridization of Dirac cones. Flat bands and the corresponding large density of electronic states can lead to novel strongly correlated phenomena. Indeed, since the recent discovery of correlated insulators and unconventional superconductivity in TBG,[9, 10, 11, 12] van der Waals multilayer stacks have been further explored as a platform of exotic correlated physics. In particular, effectively 2D heterostructures consisting of flat sheets of graphene, transition metal dichalcogenides, and hexagonal boron nitride have been successful candidates for the moiré-induced correlated phenomena.[13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] Recent experimental progress in studying correlations in multilayer heterostructures with more than two twisted graphene layers[13, 24, 25, 26] has led to a search for novel multilayer platforms with a particular focus on the trilayer geometry.[27, 14, 28, 29, 30]

In this work, we provide a detailed ab initio study of a unique extension of the TBG system: the twisted graphene sandwich (Fig. 1b), which is a promising construct of a three-layer graphene heterostructure.[28] In general, different trilayer systems are also represented by the untwisted $ABC$ stack, the twisted monolayer on bilayer, and the doubly incommensurate twisted trilayer. The $ABC$ graphene stack has been well understood,[1] although it has recently been observed to host correlated states when placed on hexagonal boron nitride due to a lattice-mismatch induced moiré superlattice.[13] The twisted monolayer on bilayer system is of the same experimental complexity as the twisted bilayer, and can host both parabolic and Dirac cone bands near the Fermi energy,[31, 14] but has a less robust flat band at its magic angle regime due to its reduced symmetry (Fig. 1c). The doubly incommensurate twisted trilayer is challenging to model accurately due to a complicated umklapp scattering process mediated by two independent twist angles, and the existing studies do not show the same spectacular flat bands as in the twisted bilayer.[29]

In contrast, the graphene sandwich retains a high degree of symmetry.[28] In the meantime, it has an effectively stronger interlayer coupling between layers, promising flat bands at larger angles and thus smaller moiré length-scales, likely enhancing correlated effects. We show in this work that the trilayer system hosts a unique feature compared to the twisted bilayer: a symmetry-protected Dirac cone that pierces through the magic-angle flat band. However, the system poses an experimental challenge: to perfectly mimic the moiré pattern of TBG, the bottom and top layers of TSWG need to be aligned in $AA$ stacking. We show that this challenge is overcome by natural relaxation of the sandwiched heterostructure, leading to protected coexistence of ultraheavy and ultrarelativistic Dirac quasiparticles at the same energy scale.

Results

Electronic structure.— The electronic states of the twisted graphene sandwich consist of two main features near the Fermi energy: a set of four nearly flat bands, similar to those found in twisted bilayer graphene, and a Dirac cone reminiscent of monolayer graphene. Much like twisted bilayer graphene, the flatness of the first feature depends sensitively on the twist angle, crystal relaxation, and external perturbations.[32, 33] The sandwich’s Dirac cone is nearly identical to that of the monolayer cone, with a Fermi velocity of $0.75\times 10^{6}$ m/s compared to $0.81\times 10^{6}$ m/s in the monolayer case. We note that the sandwiched trilayer graphene has an advantage compared to the twisted monolayer on bilayer graphene (BG/MG), where the effects of band flattening are obscured by non-symmetric band hybridization due to the absence of a layer-inversion symmetry between the monolayer and bilayer materials (see Fig. 1). In contrast, the magic angle graphene sandwich shows a Dirac cone piercing nearly flat bands, which in the non-interacting picture already classifies it as an unconventional semimetal. In our ab initio calculations, the principal magic angle is found at $1.61^{\circ}$ (see Fig. 1) for a rigid system, and inclusion of realistic lattice relaxation effects sharpens it to $1.47^{\circ}$ (see Fig. 2).

The electronic structure of the TSWG can be explained by considering the top and bottom layers as one effective layer. We assume that the top and bottom layers are aligned to ensure that they have the identical electronic coupling to the middle layer. Then, the effective states are odd or even combinations between the $p_{z}$ orbitals of the top and bottom carbon atoms. The even combinations can couple with the middle layer, with an interlayer coupling a factor of $\sqrt{2}$ stronger than that of the twisted bilayer graphene, moving the flat band regime from $\theta=1.1^{\circ}$ to roughly $\theta=1.5^{\circ}$. The odd combination is decoupled from the middle layer, leaving a copy of the pristine monolayer Dirac cone being pierced by the flat bands. In fact, this decoupling is not exact, and we report a weak coupling between the “bread” of the sandwich (top and bottom graphene layers) of $6\,$meV in density functional theory (DFT) calculations, which is ignored in our tight-binding and continuum models in this work. The weak coupling could be important for the strongly-correlated electron physics. If included, the coupling can cause weak hybridization between the flat bands and the Dirac cone, similar to the effect of vertical displacement fields (see Fig. 4c).

In realistic TSWG, the vertex of the Dirac cone is slightly off-set from the flat bands at the $K$ point (see Fig. 2). The off-set energy $\Delta E_{K}$ (the energy difference between the Dirac cone vertex and the flat bands) is always small and positive, meaning the flat electronic bands are always piercing the Dirac cone slightly above the vertex in our study. We report that the exact value of $\Delta E_{K}$ appears to be very sensitive to the parameterization of the model. Throughout this work, we use a model for monolayer graphene that includes up to the third nearest-neighbour and accounts for strain effects.[34] This model gives $\Delta E_{K}=2$ meV. However, if using an older model with up to eight nearest-neighbours (but no strain corrections),[35] we obtain $\Delta E_{K}=10$ meV. If we modify this model by truncating the range of the coupling, $\Delta E_{K}$ reduces smoothly to the $2\,$meV result at third nearest-neighbour. To explain this strong dependence on the monolayer model, it is important to understand the origin of $\Delta E_{K}$. The Dirac cone is effectively decoupled from the twisted system, and so its vertex lies at the same energy as it would in the monolayer case, at the monolayer Fermi level. The flat bands of TSWG, however, have a modified Fermi energy $\Delta E$ due to the interlayer coupling over the moiré cell. The shift $\Delta E$ in the flat bands Fermi energy is not well documented in the existing TBG literature, as it can always be safely ignored by fixing the Fermi energy of the bilayer system to zero after a band structure calculation is performed. Yet the sandwiched graphene is different. The “monolayer” energy reference is preserved in the decoupled Dirac cone, causing a relative offset between the flat bands and the cone’s vertex. Comparison to experiments and fully self-consistent modeling is necessary to accurately assess $\Delta E_{K}$, yet the coexistence of flat bands and a weakly off-set Dirac cone at the magic angle is robust.

Lattice relaxation effects.— Lattice relaxation effects are indispensable for understanding TSWG electronic structure at small twist angles, inducing renormalization of the quasiparticle spectrum near the Fermi energy and providing robust energetic stability to the flat bands. An example of the relaxation patterns in TSWG is given in Fig. 2a. The top and bottom layers (layers 1 and 3) are in $AA$ stacking and the middle layer (layer 2) is twisted $1.47^{\circ}$ counter clockwise relative to them. The relaxation pattern is similar to that of twisted bilayer graphene: the relaxation fields form spirals around the $AA$ and $AB$ stacking regions, causing the effective twist angle between layers to change locally. The spirals around $AA$ enhance the local twist angle, while the spirals around $AB$ reduce the local twist angle. Overall, this maximizes the $AB$/$BA$ low energy stacking configuration and minimizes the area of the high-energy $AA$ stacking. The displacements in layer 2 are the opposite sign and roughly twice the magnitude compared to layer 1 or 3. This is expected as layer 2 experiences twice the interlayer potential of the other layers (see Methods).

The inclusion of atomic relaxation changes the graphene sandwich’s electronic structure (see Fig. 2d). The flattest bands occur at roughly $0.1^{\circ}$ smaller of a twist angle when compared to the unrelaxed case as shown in Fig. 1c, and the gaps on both the electron and hole side of the flat bands at the $\Gamma$ point are significantly increased. This is similar to the effects of relaxation in TBG, and has to do with changes in the relative interlayer coupling strength between $AA$ and $AB$ stacked domains as well as the pseudo-gauge fields caused by in-plane strains.[36]

Energetical stability, experimental viability— To address a technical challenge of aligning external layers towards an $AA$ stacking, we investigate the effect of graphene layers misalignment. In particular, we focus on the $\theta=1.47^{\circ}$ graphene sandwich as a representative moiré supercell with the flat band regime (Fig. 2). For a trilayer graphene, one of the interlayer shifts can be safely eliminated by a corresponding shift of the reference frame; however there is still a remaining degree of freedom which alters the electronic band structure. To be concrete, we consider the relative displacement $\bm{d}_{13}$ of layer 3 with respect to the layer 1, where $\bm{d}_{13}=0$ is a desirable condition for coexistence of flat bands with Dirac cones. We find that after shifting layer 3 with respect to layer 1 ($\bm{d}_{13}\neq 0$), the relaxation pattern does change, while the overall optimized energy does not. In fact, the TSWG system after relaxations behaves in the way that layers 1 and 3 translate to remove the initial displacement away from $AA$ stacking (Fig. 2a). We remark that this phenomenon may be general for multilayered van der Waals structures, meaning stacking misalignment will not occur easily in fabricated devices.

To confirm that the $AA$ stacking of the bread layers is robust against interlayer shifts, we fix the displacement of each layer to prevent translation back to $AA$ stacking (Fig. 2b) and calculate the total energy as a function relative displacement between layers 1 and 3 (Fig. 2c) The $AA$ stacking order has an energy barrier of 20 meV/nm2 when the displacement between layers $1$ and $3$, labeled $\mathbf{d}_{13}$, is fixed. To compare with important energy scales, we use our model for twisted bilayer graphene which gives an energy difference of 30 meV/nm2 between relaxed and unrelaxed TBG at $\theta=1^{\circ}$, and produces relaxation patterns in a good agreement with those observed in diffraction experiments.[37] We thus conclude that the TSWG relaxation barrier of 20 meV/nm2 is experimentally vital, and we expect that relaxation will cause graphene sandwich devices to naturally align the bread layers as $AA$, circumventing the experimental challenge of aligning the top and bottom layers manually. We note that the $AA$ stacked alignment is the highest symmetry configuration of the system, yielding a 3-fold rotation ($D_{3}$) center. This is a natural result, as most crystals minimize their energy by maximizing internal symmetry, except in some more exotic situations – such as Peierls distortions and charge density waves.

Perfectly flat bands piercing Dirac cones.— The behaviour of flat bands piercing the Dirac cone can be most clearly understood in the chirally-symmetric limit of the effective continuum model (see Methods). At the magic angle twist 1.5*∘*, the moiré superlattice is approximately 40 graphene unit cells, and the effective behaviour of electrons in TSWG is governed by a large-period moiré field built on the three symmetry-related wave vectors $|\bm{q}_{i}|=2k_{D}\sin({\theta}/{2})$. For the given input monolayer Fermi velocity $v_{0}$, the continuum model for TSWG is captured by three key parameters: interlayer couplings between $AA$ and $AB$ sites ($w_{AA}$ and $w_{AB}$) as well as by the twist degree of freedom ($\theta$). Out-of-plane lattice relaxation affects the relative ratio of $w_{AA}/w_{AB}$, which is strongly suppressed at small angles justifying the use of a chiral-symmetric model with $w_{AA}=0$. Importantly, this model is ruled entirely by the only dimensionless twisting parameter $\alpha(\theta)={w_{AB}}/{2k_{D}v_{0}\sin\frac{\theta}{2}}$, and shows the perfectly flat bands piercing the Dirac cone at $\alpha_{*}\simeq 0.414$ which for $w_{AB}=110$ meV corresponds to the magic angle $\theta_{*}\approx 1.55^{\circ}$ (see Fig. 3). The principle magic angle of TSWG, being exactly defined in the continuum model, is nearly $40\%$ larger than the reported magic angle in TBG, which makes TSWG experimentally attractive.

Fig.3 shows the band structures for three instances of TSWG with twist angles exactly at the principal magic angle ($\alpha_{*}=0.414..$), and just above and below this angle. We observe the perfectly flat bands piercing the Dirac cone vertex at $\theta=1.55^{\circ}$ (Fig. 3b), while for slightly different angles (Fig. 3a,c) the bands are dispersive. We further track the renormalized Fermi velocities for TSWG, one from the Dirac cone itself and one from the flat bands. We define the first Fermi velocity as the slope corresponding to the flattened bands at the $K$ point and the second Fermi velocity corresponding to the Dirac cones. The first Fermi velocity vanishes at $\theta_{1}=1.55^{\circ}$ and then reappears (see Fig 3b), thus showing the behaviour similar to the case of twisted bilayer graphene.[3, 5] On the contrary, the second Fermi velocity is constant in this model and equals to the monolayer value $v_{0}$, confirming that the Dirac cone is a robust feature of TSWG, and the second Fermi velocity is very weakly dependent on twist as is reflected in our rigorous atomistic calculations. The continuum model for the twisted graphene sandwich also predicts higher-order magic angles (e.g. $\alpha_{2}\simeq 1.57$, and $\alpha_{3}\simeq 2.65$, which corresponds to $\theta_{2}\approx 0.405^{\circ}$ and $\theta_{3}\approx 0.240^{\circ}$), but these have not been confirmed by ab initio calculations and are likely suppressed by lattice relaxations as in the bilayer case.[33] We remark that the flawlessly flat band TSWG model ($w_{AA}=0$) captures the principle magic angle accurately ($1.55^{\circ}$ vs $1.47^{\circ}$ in atomistic calculations with relaxations), but stretches the energy scales. The inclusion of a realistic $w_{AA}\approx 90$ meV produces the energy scales similar to ab initio band structures, preserving the magic angle value.

The continuum model is also interesting because the flat band condition can be derived analytically up to an arbitrary precision. The leading order perturbation theory in $\alpha$ gives the following dependence of the flat band Fermi velocity on twist $\theta$:

[TABLE]

The principal magic angle is thus precisely defined when the two Dirac cones hybridize into the flat band with $v(\theta_{*})=0$, which has a solution at

[TABLE]

This estimate is also valid for a more realistic case of $w_{AA}\neq 0$. The chiral symmetric model with $w_{AA}=0$ also has an exact mapping to the twisted bilayer graphene[28], and thus, in continuum theory, the graphene sandwich can have an infinite set of magic angles related to the TBG magic angles by the factor $\sqrt{2}$. For our discussion, only the principal magic angle is relevant.

Tunability under external electric fields.— An advantage of the graphene sandwich is the tunability under external fields, with which the energy off-set between the Dirac cone and the Fermi level can be precisely controlled. We show this with application of external fields and uniform strains to the unrelaxed structure to understand how the position of the Dirac vertex could be tuned in experimental devices. We find that uniform planar strain can smoothly tune $\Delta E_{K}$, with a larger (smaller) lattice giving a smaller (larger) gap, which happens essentially due to rescaling of the system’s characteristic energy. As the lattice expands, the nearest-neighbour bonding distances increase and the electronic couplings become weaker. The Fermi velocity of the Dirac cone is directly proportional to this nearest-neighbour coupling energy, and as previously discussed, the Fermi velocity sets an overall energy scale for the twisted graphene systems. Thus, as the bonds are made weaker or stronger, the energy scale becomes smaller or larger, which directly relates to the value of $\Delta E_{K}$.

Applying an external vertical electric displacement field can control the magnitude of $\Delta E_{K}$ and tune the amount of hybridization between the flat bands and the Dirac cone (Fig. 4). As mentioned earlier, this hybridization is set to be zero in our simplified tight-binding model in the absence of displacement electric field, but we expect some non-zero hybridization at zero field because of the weak electronic coupling between the bread layers. When the displacement field is stronger than 0.4 V/nm, the low-energy band structure is no longer easily comparable to the zero field case, as the flat bands become very dispersive and higher-energy bands begin intersecting the original flat bands near the Fermi energy. To compare the displacement field modeled here to experimental devices, the dielectric screening of the graphene layers and any encapsulating substrate should be considered. Regardless, sufficiently weak external electric fields provide tunable control of the relative position of the Dirac cone vertex and the intensity of the van Hove singularity associated with the flat bands.

Discussion.— We first discuss the implication of our TSWG results for experiments. Most importantly, the magic angle in TSWG is 1.5*∘, which is $40\%$ larger than in the parent TBG heterostructure. In general, this is advantageous for two reasons: first, it is easier to fabricate a multilayer heterostructure at larger twist angles; second, a qualitative trend is that larger angles (smaller moiré patterns) generally correspond to a higher superconducting $T_{C}$. As an example, in TBG the magic angle is 1.1∘* and $T_{C}=1.7$ K, but when applying appropriate hydrostatic pressure the magic angle superconductivity can be observed at larger angle (1.27∘) and larger $T_{C}$ (3.5 K).[11] Another example is the twisted double bilayer graphene (TDBG) which has the flat band region around 1.25*∘*, and superconducting $T_{C}=3.5$ K.[26, 24, 25] Additionally, the twisted sandwich graphene can be fabricated from the single flake of graphene. First, one makes the twisted bilayer graphene with a standard approach.[38] At small angles this heterostructure relaxes on the moiré length scale, which will assist the deposition of a third layer at the correct $0^{\circ}$ alignment with the first layer. As we have shown in Fig. 2(c), even if the third layer is not positioned perfectly, it will tend to relax towards the energy minimum with depth of order 20 eV/nm2. We note that this energy barrier is relatively large, as it is the same order of magnitude as between unrelaxed and fully relaxed TBG.

The flat bands in the graphene sandwich exhibit a strong van Hove singularity in the density of electron states (see Fig. 2d). If realized experimentally, this will promote strongly correlated electron states as the kinetic energy is suppressed. Therefore, we may expect the emergence of correlated insulation and unconventional superconductivity under fine tuning, similar to TBG and TDBG. One of the main result of this study is that although at first glance the TSWG seems challenging to fabricate, our ab initio results predict that the system will always relax towards the beneficial $AA$ stacking between the first and third layers, providing the proper atomic geometry for the remarkably flat bands.

Last but not least, TSWG hosts both the stable flat bands and a Dirac cone in close proximity to one another in the band structure. There are only a few systems with this property, such as the exotic Kondo Weyl Semimetals and some Kagomé lattice systems with relatively flat bands offset from Dirac cones. All these systems are difficult to realize experimentally. From this perspective, TSWG may open a feasible experimental path for realizing the coexistence of both strongly localized and ultramobile quasiparticles simultaneously, important for the “steep band/flat band” scenario of superconductivity. As the flat bands can pierce the Dirac cone extremely close to its vertex, TSWG can also be viewed as a potential platform for reaching the so-called triple point states, a fine tuning for which the low-energy physics is effectively described by the Dirac equation with pseudospin-1 (see also Refs.[39, 40]). To conclude, the twisted sandwiched graphene represents a novel, experimentally feasible platform for a broad range of exotic electronic phenomena.

1 Methods

Density functional method.– We use the VASP [41, 42, 43] implementation of density functional theory (DFT) to calculate electronic structure for untwisted bilayer and trilayer graphene systems. The semi-local meta-GGA functional SCAN+rVV10 [44] is used for its good performance in van der Waals materials and low computational cost. Multiple calculations of these untwisted systems are performed, with the graphene layers shifted in-plane to accurately capture electronic and mechanical effects of different stacking orders. The graphene sandwich (bilayer) systems have 6 (4) carbon atoms, and we include a vertical vacuum space of $20\,\textrm{\AA}$ to prevent interactions between periodic images in the $z$ direction of the heterostructures. Then ab initio tight-binding parameters are extracted by the method of maximally localized Wannier functions with the wannier90 package. From this, one obtains a fully parameterized model for in-plane and inter-plane $p_{z}$ orbital interactions, which can also accurately account for corrugations and in-plane strain.[35, 34, 45, 36] We find that the largest effective tight-binding coupling between $p_{z}$ orbitals of the top and bottom layers of the sandwich is $6$ meV, roughly $2\%$ of the maximum coupling between adjacent layers. For this reason, we ignore couplings between the top and bottom layers in tight-binding and continuum simulations.

Lattice relaxations modelling.–The relaxation of TSWG is obtained using a continuum model to account for in-plane distortions due to relaxation and the Generalized Stacking Fault Energy (GSFE) to account for the interlayer coupling. The relaxation of the layer $i$, $\bm{u}_{i}(\bm{b})$, is defined in terms of the local configuration or the relative local stacking, $\bm{b}$, and we obtain $\bm{u}_{i}(\bm{b})$’s by minimizing the total energy.[46] The energy has two contributions. The first is the intralayer energy of the $i$-th layer, which is calculated based on linear elasticity theory,

[TABLE]

where $\Gamma$ is the union of all configurations, $\mathcal{E}(\gradient\bm{u}_{i})$ is the strain tensor, $C_{i}$ is the linear elasticity tensor of the $i$-th layer (which is identical for all $i$’s in this case), $G$ and $K$ are shear and bulk modulus of a monolayer graphene, which we take to be $G=47352\,\mathrm{meV/cell}$, $K=69518\,\mathrm{meV/cell}$ and the graphene unit cell size is $5.3128\,$Å2. The values of $G$ and $K$ are obtained with DFT calculations by isotropically straining and compressing the monolayer and performing a linear fitting of the ground-state energy as a function of the applied strain. The second energy contribution is the interlayer energy, which is described by the GSFE.[36, 47, 48] The GSFE, denoted as $V_{\mathrm{GSFE}}(\bm{b})$, has been employed to explain relaxation in van der Waals heterostructures,[36] which depends on the relative stacking between two adjacent layers. We obtain the $V_{\mathrm{GSFE}}$ by applying a $9\times 9$ grid sampling of rigid shifts to layer 1 in the unit cell with respect to layer 2 and extract the relaxed ground state energy at each shift from DFT. The GSFE of graphene at a given configuration $\bm{b}=\begin{pmatrix}v&w\end{pmatrix}^{T}$ can then be expressed as follows,

[TABLE]

where $c_{i}$’s are coefficients found by fitting the ground state energy at each shift. $c_{0}=6.832\,\mathrm{meV/cell}$, $c_{1}=4.064\,\mathrm{meV/cell}$, $c_{2}=-0.374\,\mathrm{meV/cell}$, $c_{3}=-0.095\,\mathrm{meV/cell}$. In terms of the $V_{\mathrm{GSFE}}$, the total interlayer energy $E_{\mathrm{inter}}$ can be then written as follows for a relaxed TSWG,

[TABLE]

where the factor of 2 comes from the sum of couplings between layers 1, 2 and layers 2, 3, and we use the fact that $\bm{u}_{1}(\bm{b})=\bm{u}_{3}(\bm{b})$ due to the layer inversion symmetry in the TSWG system. Note that the $V_{\mathrm{GSFE}}$ is a function of the sum of unrelaxed configuration and the relaxation displacement vectors in order to describe the inter-layer stacking energy after relaxation. The total energy is the sum of the interlayer and the intralayer energies:

[TABLE]

The relaxation $u_{i}(\bm{b})$ is computed by minimizing the total energy. The following linear transformation maps the relaxation from the local configuration to the real space positions $\bm{r}$:

[TABLE]

where $E_{1}$ and $E_{2}$ are the unit cell vectors of the first (unrotated) and the second layers (rotated counter clockwise by $\theta$), respectively.

Tight binding calculations.–We take $a_{1}=\alpha(1,0)$ and $a_{2}=\alpha(\frac{1}{2},\frac{\sqrt{3}}{2})$ with $\alpha=2.4768\textrm{\AA}$ as the unit cell vectors of graphene. Periodic supercells are constructed in terms of the integers $n$ and $m$ according to $\cos\theta=\frac{n^{2}+4nm+m^{2}}{2(n^{2}+nm+m^{2})}$. For the trilayer sandwich, the bottom and top layers are untwisted while the middle layer is twisted counterclockwise by $\theta$. This gives supercell vectors of $(na_{1}-ma_{2})$ and $(-ma_{1}+(m+n)a_{2})$ for the untwisted layers and supercell vectors of $(ma^{\prime}_{1}-na^{\prime}_{2})$ and $(-na^{\prime}_{1}+(m+n)a^{\prime}_{2})$ for the twisted layer ($a^{\prime}_{i}$ are $a_{i}$ rotated c.c.w by $\theta$). Band structures are calculated by diagonalization of these supercell Hamiltonians. The interlayer electronic couplings deal with relaxation and strain effects easily, as they depend directly on the $p_{z}$ orbital positions by design. Strain is included for in-plane coupling with a simple bond-length approximation, where the coupling $t$ is dependent on the bond-length $r$ by $t=t_{0}+\alpha_{t}\frac{r-r_{0}}{r_{0}}$.[34] The vertical electric displacement field effect is treated in the leading order by adding on-site energies $\delta=E_{z}r_{z}$ in the existing tight-binding model.

Continuum model. We use an effective continuum model based on approach of Refs[5, 28]. In the pristine setting, the effective continuum model for the twisted graphene sandwich can be described by

[TABLE]

where the moiré-induced interlayer coupling is taken up to the first shell in momentum space,

[TABLE]

with the three-fold star of $\bm{q}_{i}$, $|\bm{q}_{j}|=2k_{D}\sin(\theta/2)$, each equirotated by $\phi=2\pi/3$, and

[TABLE]

where $\bm{\mathcal{G}}^{(0)}_{\theta}=0$, $\bm{\mathcal{G}}^{(1)}_{\theta}=\bm{q}_{2}-\bm{q}_{1}$, $\bm{\mathcal{G}}^{(2)}_{\theta}=\bm{q}_{3}-\bm{q}_{1}$ are the moiré reciprocal cell vectors and $\bm{d}$ is the relative displacements of one layer with respect to another one, and

[TABLE]

Note that in the continuum model, one can eliminate one of the displacements after redefining the reference frame. The Hamiltonian of Eqn. (8) acquires additional chiral symmetry and perfectly flat bands piercing Dirac cones at neutrality for $w_{AA}=0$. In this idealistic setting, the twisted trilayer graphene have a family of well-defined magic angles, enlarged compared to the similar sequence in the TBG by $\sqrt{2}$. See Ref.[28] for further details.

{addendum}

We thank Pablo Jarillo-Herrero, Bertrand Halperin and Oleg Yazyev for useful discussions. This work was performed under Grant No. P2ELP2_175278 from the Swiss National Science Foundation (SNSF), Grants No. DMR-1664842 and No. DMR-1231319 from the National Science Foundation (NSF), and ARO MURI Award W911NF-14-0247. The computations were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.

The authors declare no competing financial interests.