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 • , 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. 14,[27][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, 14,31 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 symmetryprotected 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.

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 × 10 6 m/s compared to 0.81 × 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 • (see Fig. 1) for a rigid system, and inclusion of realistic lattice relaxation effects sharpens it to 1.47 • (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 √ 2 stronger than that of the twisted bilayer graphene, moving the flat band regime from θ = 1.1 • to roughly θ = 1.5 • . 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 ∆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 ∆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 ∆E K = 2 meV. However, if using an older model with up to eight nearest-neighbours (but no strain corrections), 35 we obtain ∆E K = 10 meV. If we modify this model by truncating the range of the coupling, ∆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 ∆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 The inclusion of atomic relaxation changes the graphene sandwich's electronic structure (see Fig. 2d). The flattest bands occur at roughly 0.1 • 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 Γ 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 θ = 1.47 • 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 d 13 of layer 3 with respect to the layer 1, where 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 (d 13 = 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/nm 2 when the displacement between layers 1 and 3, labeled 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/nm 2 between relaxed and unrelaxed TBG at θ = 1 • , 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/nm 2 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. which for w AB = 110 meV corresponds to the magic angle θ * ≈ 1.55 • (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 (α * = 0.414..), and just above and below this angle. We observe the perfectly flat bands piercing the Dirac cone vertex at θ = 1.55 • (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 θ 1 = 1.55 • 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. α 2 1.57, and α 3 2.65, which corresponds to θ 2 ≈ 0.405 • and θ 3 ≈ 0.240 • ), 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 • vs 1.47 • in atomistic calculations with relaxations), but stretches the energy scales. The inclusion of a realistic w AA ≈ 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 α gives the following dependence of the flat band Fermi velocity on twist θ:

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

This estimate is also valid for a more realistic case of w AA = 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 √ 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 ∆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 ∆E K .

Applying an external vertical electric displacement field can control the magnitude of ∆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. [24][25][26] 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 • 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/nm 2 . 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 perspect-ive, 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.

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 Å 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. [34][35][36]45 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, u i (b), is defined in terms of the local configuration or the relative local stacking, b, and we obtain u i (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,

where Γ is the union of all configurations, E(∇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 meV/cell, K = 69518 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 GSFE (b), has been employed to explain relaxation in van der Waals heterostructures, 36 which depends on the rel-ative stacking between two adjacent layers. We obtain the V GSFE by applying a 9 × 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 b = v w T can then be expressed as follows,

where c i 's are coefficients found by fitting the ground state energy at each shift. c 0 = 6.832 meV/cell, c 1 = 4.064 meV/cell, c 2 = -0.374 meV/cell, c 3 = -0.095 meV/cell. In terms of the V GSFE , the total interlayer energy E inter can be then written as follows for a relaxed TSWG,

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 u 1 (b) = u 3 (b) due to the layer inversion symmetry in the TSWG system. Note that the V 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:

The relaxation u i (b) is computed by minimizing the total energy. The following linear transformation maps the relaxation from the local configuration to the real space positions r:

where E 1 and E 2 are the unit cell vectors of the first (unrotated) and the second layers (rotated counter clockwise by θ), respectively.

Tight binding calculations.-We take a 1 = α(1, 0) and

2 ) with α = 2.4768 Å as the unit cell vectors of graphene. Periodic supercells are constructed in terms of the integers n and m according to cos θ = 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 θ. 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 1 -na 2 )

and (-na 1 + (m + n)a 2 ) for the twisted layer (a i are a i rotated c.c.w by θ). 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 + α t r-r 0 r 0 . 34 The vertical electric displacement field effect is treated in the leading order by adding on-site energies δ = 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

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

T n e -iqnr , with the three-fold star of q i , |q j | = 2k D sin(θ/2), each equirotated by φ = 2π/3, and

where

θ = q 3 -q 1 are the moiré reciprocal cell vectors and d is the relative displacements of one layer with respect to another one, and

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 √ 2. See Ref. 28 for further details.