Extended magic phase in twisted graphene multilayers

Theoretical and experimental studies have verified the existence of ``magic angles'' in twisted bilayer graphene, where the twist between layers gives rise to flat bands and consequently highly correlated phases. Narrow bands can also exist in multilayers with alternating twist angles, and recent theoretical work suggests that they can also be found in trilayers with twist angles between neighboring layers in the same direction. We show here that flat bands exist in a variety of multilayers where the ratio between twist angles is close to coprime integers. We generalize previous analyses, and, using the chiral limit for interlayer coupling, give examples of many combinations of twist angles in stacks made up of three and four layers which lead to flat bands. The technique we use can be extended to systems with many layers. Our results suggest that flat bands can exist in graphene multilayers with angle disorder, that is, narrow samples of turbostatic graphite.

Introduction -Twisted bilayer graphene (tBG) has proven to be a fertile platform for condensed matter physics studies owing to the remarkable tunability of the Fermi velocity and electronic structure through the twist angle [1].In particular, at the so-called "magic angles" [2] for tBG, θ M tBG ≈ 1.1 • , the Fermi velocity vanishes, allowing for strongly correlated phases where interactions dominate over the kinetic energy [3,4].This has led to a resurgence of interest in moiré systems, including but not limited to Van der Waals heterostructures [5][6][7][8], transition metal dichalcogenides [9][10][11][12], and of course graphene multilayers [13][14][15][16].One particularly successful avenue of investigation has been the extension of the arguments demonstrating the existence of magic angles in tBG [17,18] to graphene multilayers [19], which found simple geometric relations between the bilayer and multilayer magic angles, predictions that were later verified experimentally [20][21][22].Later numerical studies have built upon these foundations, demonstrating for example that the existence of flat bands is robust against perpendicular magnetic fields [23] and asymmetry of the layer parameters [24], and that peaks in the density of states (DOS) appear, and are resistant to relaxation, for trilayer configurations beyond the canonical alternating-twist case [25].Recent work has predicted magic angles in the equal twist case [26], and a magic line, with magic angle a function of pressure, has been predicted in twisted WSe 2 [27].In this work we demonstrate the existence of flat bands in multilayers with generic, commensurate twist angles, indicating an extended N − 2 dimensional magic surfaces in the space of N − 1 twist angles characterising N layers, with particular emphasis on generic twisted trilayer graphene (tTG).
A prior study [19] noted that in the particular case of unshifted, alternating twist graphene multilayers, the Hamiltonian may be exactly decomposed into a direct sum of bilayer graphene Hamiltonians (with an additional monolayer for odd number of layers), with the moiré potential identical up to a scalar multiple.This allows one to directly read off the magic angles of multilayer stacks.In contrast, other trilayer configurations have more complex effective moiré potentials, though magic angle physics may still be possible [28].In particular, a strongly correlated superconducting state has been observed in a trilayer system with incommensurate twist angles [29].
Moiré lattices in twisted graphene trilayers and in other twisted stacks.-The electronic structure of a twisted bilayer can be reasonably approximated by a periodic moiré structure [1,2], irrespective of whether the two layers are commensurate at the atomic scale.The situation is more complex in twisted stacks with more than two layers [25,[30][31][32][33][34][35][36].In a trilayer, for instance, two twist angles between nearest neighbor layers can be defined, θ 12 and θ 23 , shown in Fig. 1(a).Except for the case mentioned above of multilayers where the twist angle is of equal magnitude and opposite sign between nearest neighbor pairs of layers [19], no simple global moiré structure can be defined.
An approximate moiré lattice can be defined if the relative angle between the Brillouin Zones defined by the twist angles between different pairs of layers is neglected [32,34,35].For instance, in a trilayer with two different twist angles, θ 12 , θ 23 , this approximation becomes exact when θ 12 = θ 23 in the limit θ 12 → 0, although corrections [36] are required for physically relevant angles, θ 12 , θ 23 ≈ 1 • .It is worth noting that this approximation, for a trilayer, has a dependence not only on the angles θ 12 , θ 23 , but also on the relative displacement of the moiré potentials, D illustrated in Fig. 1(b) [33].The dependence is periodic on D, with wavelength of the  are used for clarity, and less misalignment is expected for real magic angle samples.The lengthscale of this modulation is determined by the closest approach of second layer lattice sites formed from the two moiré patterns, indicated with a dotted line.The local expansion used in numerical calculation neglects this misalignment and, for commensurate angles, collapses this distance to zero.(d) Log-scale plot of the moiré-of-moiré length scale, lm2, in microns, as a function of twist angles.The left diagonal in white marks the line θ12 = −θ23, where there is no misalignment and lm2 diverges.lm2 also diverges when either angle is zero and there is no moiré-of-moiré pattern to speak of.order of the moiré length scales.
It can be shown [37] that the approximation mentioned above defined by a moiré lattice commensurate with the two twist angles and a displacement between the top and bottom layers describes locally the full problem of the trilayer.For the case of approximately commensurate twist angles, θ 12 ≈ mθ 0 , θ 23 ≈ nθ 0 , where m, n are coprime integers, the interlayer tunneling terms in the trilayer continuum Hamiltonian can be written as containing two periodicities: i) the periodicity, ℓ m , of the moiré problem obtained by neglecting the angle between the two Brillouin zones defined by θ 12 , θ 23 , of order d/θ 0 , where d is the lattice constant of graphene, and ii) a periodicity, ℓ m2 obtained by modulating the displacement D, over a larger unit cell whose dimensions scale as d/θ 2 0 .Each tBG pair in tTG defines a momentum-space moiré unit cell, purple and orange in Fig. 1(c), whose misalignment leads to misalignment of the moiré lattices, shown with the black lines and the finite δk.Neglecting this moiré misalignment forces the black lines to coincide, and for commensurate angles forces δk = 0. δk therefore is conjugate to ℓ m2 .For small values of θ 0 , shown in Fig. 1(d), this scaling behavior implies that ℓ m ≪ ℓ m2 , as ℓ m2 can reach values of microns for θ 0 ≲ 1 • .Further details may be found in [38].
It is worth noting that small and layer-dependent biaxial strains can exactly align the Brillouin zones associated to different twist angles, θ 12 , θ 23 , so that ℓ m2 → ∞ and the approximation which makes use of a single moiré lattice defined by ℓ m becomes exact, see [39].These strains imply a small relaxation of the lattice [39].This analysis has been extended to some commensurate cases of trilayers with unequal twist angles [40].
Results.Flat bands in the chiral limit.-We first analyze numerically the existence of flat bands in chiral model trilayers [17,18] for different combinations of twist angles, θ 12 , θ 23 and relative displacements between the top and bottom layers.The chiral limit for graphene trilayers with the same twist angle between nearest layers, θ 12 = θ 23 has been recently explored in [26,39,41].The parameters that we use [42] in the following are γ 0 = 2.6 eV (intralayer nearest neighbor hopping), and t AB = 0.0975 eV (interlayer hopping).The choice of γ 0 leads to a Fermi velocity v F ≈ 840 km/s.
The analytical arguments for the existence of infinitely flat central bands in twisted bilayer graphene [18] is based on the existence of zeros at the corners of the unit cell of a 2 component spinor which describes the periodic part of the wavefunction at the K point of the Brillouin Zone.As one of the two components, by symmetry, vanishes at the relevant points, the problem reduces to the tuning of the zeros of a single component.The argument can be generalized to a trilayer, where a three component spinor is involved.When the Hamiltonian has C 3 symmetry [18] two of the three components of the spinor vanish at the relevant point in the real space unit cell, so that the existence of infinitely flat bands reduces to the existence of a zero in a single component of the spinor.The C 3 symmetry exists at least for a displacement D = 0, and for the displacements which take the trilayer from the AAA to the ABA configuration.An alternative analysis, for the case of a trilayer with θ 12 = θ 23 , can be found in [26].We numerically find flat bands, whose width is continuously lowered as the angle is determined with better precision, for all combinations of integers, m, n and displacements, D, tried.Results are shown in Fig. 2. Interestingly, in some cases we find a quadruplet of bands nearly degenerate at zero energy.We show in Fig. 2 the local DOS at the magic angles corresponding to {m, n} = {1, 1} and {m, n} = {3, −4} obtained for a displacement between top and bottom layers D = 0.For each magic angle, we plot the DOS for the selected displacement, D = 0, and for two other displacements, shown in the inset in Fig. 2 (a).There is a significant dependence of the local DOS on D in the {1, 1} case, and a less pronounced dependence in the {3, −4} case.This trend towards a weaker dependence on D with increasing values of the integers {m, n} is observed for all combinations explored, see [38].
Alternatively, the 3×3 matrix equation which describes the spinor at the magic angles of the chiral bilayer [18,26] can be approximately reduced to an effective bilayer, making significantly easier the calculation of magic angles as function of commensuration and displacement, see [38].
The scheme discussed here for a trilayer [32,34,35,37], based on two twist angles, two integers, and one interlayer displacement can be extended to an N layer model, with N − 1 twist angles, and N − 2 independent displacements.For the case of N −layer stacks where the twist angles between neighboring layers are equal, an analytical expression for the first magic angle, which generalizes the results in [26], gives: ] where e is Euler's number.Figure 2(d) shows results for magic angle combinations in stacks of four layers, defined by three twist angles, and three integers, l, m, n.Calculations for the chiral helical case, θ 1 = θ 23 = θ 34 = θ 0 , can be found in [38].
The magic angles where flat bands appear depend on the integers m, n, and on the displacement D. To a first approximation, we can approximate the electronic density of states of the system by an average over the values of D. This approximation can be expected to become exact in the limit ℓ m2 /ℓ m ∼ 1/θ 0 → ∞, and also when the dependence of the band structure on D is small.
The flat bands at the magic angles reported here show a greater degree of plasmon screening compared to tBG and thus, a greater resilience against distortion from interaction effects.They also have a finite Berry curvature and, typically, non zero Chern numbers, see [38].
Results.Analysis of the trilayer studied in [29].We now analyze the experimental situation reported in [29]  with [29]).As shown in the figure, the dependence of the density of states on the displacement D is very weak, which suggests that the weighted average shown in Fig. 3 (d) gives a good approximation to the exact value.
Figure 4(a) is the tight-binding (TB) density of states for trilayer with θ 12 ≈ 1.42 • , θ 23 ≈ −1.88 • , which is obtained via the parameters γ 0 = 3.1 eV and γ 1 = 0.43 eV [43][44][45], the same as the ones in [29].The discrepancy between the TB and chiral limit results is due to the difference of the Fermi velocity.Moreover, if we reduce the Fermi velocity, a flat band appears in the charge neutral- ity point (not shown here).The states of the peak at the energy of 7 meV in Fig. 4(a) is highly localized with an independence of the displacement D, shown in Fig. 4(c), which is in agreement with the results in Fig. 3 in the chiral limit.
We extend our discussion to another configuration, the helical trilayer graphene with θ 12 = θ 23 = θ 0 .One magic angle is θ 0 = 2 • [25], which has peaks located at the charge neutrality point, shown in Fig. 4. If we consider the lattice relaxation [46][47][48], the magic angle reduces to 1.57 • [39].For AAA helical trilayer graphene with identical twist angles, the moiré-of-moiré structure locally relaxes into large areas that contains only one moiré pattern.That is, a large region of the periodic ABA could be reconstructed if a strong intralayer and interlayer potential is utilized.The states of the flat band in the relaxed system is highly dependent on the displacement D (see [38] Fig. 10), which is consistent with the results in Fig. 2(a).Furthermore, the helical trilayer graphene with θ 12 = θ 23 could be an ideal platform for realization of the chiral limit due to significant corrugation, or via heterostrain engineering [39].
Conclusions.We have considered the existence of flat bands in twisted stacks of graphene layers where the twist angles do not allow for the definition of a simple moiré structure, generalizing previous works [26,32,34,35,37,[39][40][41] which have focused on special configurations yielding exact or approximate moiré structures in trilayers.By analyzing the chiral limit [17,18] we conclusively demonstrate that flat bands appear for all combinations of twist angles where the angles are close to multiples of coprime integers.Our further extension to generic, incommensurate angles [38] connects the previously disparate magic points of specific trilayer configurations into an extended magic line.Building on the local approximation as developed in [37] for the specific case of helical trilayers with nearly equal twist angles, we accurately describe the electronic structure of generic twisted multilayers even away from the magic angle, with our results consistent with both large scale tight binding calculations [25] and contemporary experiments [29].The analysis presented here opens the way to the study of topological features in many types of twisted multilayers, and they suggest that the rich phase diagram found in twisted bilayer graphene may exist in many other graphene stacks.
Note added.We have noticed ref. ([54], posted a day before our paper.This reference discusses topics related to our work.Insofar the two manuscript overlap, the results are in agreement.from the "Severo Ochoa" Programme for Centres of Excellence in R&D (CEX2020-001039-S / AEI / 10.13039/501100011033). F.G. acknowledges funding from the European Commission, within the Graphene Flagship, Core 3, grant number 881603 and from grants NMAT2D (Comunidad de Madrid, Spain), SprQuMat (Ministerio de Ciencia e Innovación, Spain) and financial support through the (MAD2D-CM)-MRR MATE-RIALES AVANZADOS-IMDEA-NC.DCWF, MA, LP and SA acknowledge support from the Singapore National Research Foundation Investigator Award (NRF-NRFI06-2020-0003).Z.Z.acknowledges support funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 101034431.Numerical calculations presented in this paper have been performed on the supercomputing system in the Supercomputing Center of Wuhan University.
Supplemental Material for: Extended magic phase in generalized trilayer graphene DEFINITION OF PERIODICITIES.
For completeness, we specify the two relevant periodicities in a twisted graphene trilayer with almost commensurate twist angles, θ 12 and θ 23 following the steps described in [37].We write the two angles as: where m, n are coprime integers, and we assume that θ 0 ≪ 1 and δθ ≪ θ 0 .The angles θ 12 and θ 23 define moiré structures.The vectors connecting the points ⃗ K and ⃗ K ′ of the respective Brillouin zones are: where we have assumed that θ 0 ≪ 1.We obtain: Vectors ⃗ g 0 and δ⃗ g and three equivalent vectors rotated by ±120 • define the relevant periodicities of the trilayer.The two vectors satisfy |δ⃗ g| ≪ | ⃗ g 0 |, so that, in real space, the periodicity associated with δ⃗ g changes much more slowly than that asociated to ⃗ g 0 .Locally, ⃗ g 0 defines a moiré lattice with unit vector ℓ m ≈ d/θ 0 , where d is the lattice constant of graphene.
In the following, for simplicity, we assume that δθ ≪ θ 2 0 , although this assumption is not crucial.Then, δ⃗ g defines a new periodicity, with unit cell rotated 90 • with respect to the unit cell defined by ⃗ g 0 and a lattice unit of length ℓ m2 ≈ (4d)/(|m − n|θ 2 0 ) ≫ ℓ m .Note, finally, that our definition of δ⃗ g, for m ̸ = 1 or n ̸ = 1, differs from the one used in [37], although both definitions scale as |δ⃗ g| ∝ θ 2 0 .The analysis described above can be extended to an arbitrary number, N , of layers, and N 1 twist angles between nearest layers, provided that these angles are almost commensurate, defined by N − 1 coprime integers.Such a calculation leads to a periodicity, ⃗ g 0 , defined by an angle θ 0 , and N − 2 additional periodicities defined by vectors Insofar as the role of δ⃗ g can be reduced to an overall displacement between the layers (see below) the two definitions should give the same results.

SLOW PERIODICITY AS A DISPLACEMENT BETWEEN LAYERS.
We now elaborate on the argument that the slow periodicity, δ⃗ g, identified in the previous section can be interpreted as a local displacement of the top layer with respect to the bottom layer [37].The interlayer tunneling part of the continuum Hamiltonian [1,2], in real space, is: T j e im⃗ gj •⃗ r where T j define inter sublattice 2 × 2 matrices: and: The momenta of the Bloch waves used in the calculation are shown in Fig. S1.We can generalize the tunneling Hamiltonian to the case where the moiré pattern of layers 1 and 2 is displaced by a shift ⃗ D 1 and the moiré pattern between layers 2 and 3 is displaced by a shift ⃗ D 3 [33].In the initial Hamiltonian, the rotation is around a point with AA stacking between layers 1 and 2, and AA stacking between layers 2 and 3, that is, AAA stacking.The shifts lead to the Hamiltonian: T j e im⃗ gj •(⃗ r− ⃗ D1) Each Bloch wave in the calculation, see Fig. S1, can be modified by a gauge factor, e i ⃗ G• ⃗ D ′ , where ⃗ G is the momentum of the wave.This gauge transformation induces a phase in the interlayer tunnelings, which is equivalent to a redefinition of the shifts: This result is independent of the integers m and n, and of the angle θ 0 .An equal shift, ⃗ D 1 = ⃗ D 3 , of the two moiré patterns can be gauged away.This implies that the calculated electronic structure is only dependent on the value of The two periodicities defined by ⃗ g 0 and δ⃗ g in the previous section modify the continuum Hamiltonian, Eq. (S4), leading to: T j e im⃗ g0j •⃗ r e imδ⃗ gj •⃗ r T j e in⃗ g0j •⃗ r e −inδ⃗ gj •⃗ r (S9) The vectors δ⃗ g j define a periodicity with unit vectors in real space: where ⃗ R 1 and ⃗ R 2 are the lattice vectors of an untwisted graphene layer.We now assume that the arguments δ⃗ g j • ⃗ r vary slowly with respect the functions ⃗ g 0j ⃗ r, and assume that the value of ⃗ r = ⃗ D in δ⃗ g j • ⃗ r is constant.Then, using Eq.(S7) we obtain that the local Hamiltonian is equivalent to displacements of the top and bottom layers by ±(d/ℓ m2 )R π/2 ( ⃗ D), where R θ is a rotation operator.For m ̸ = 1 or n ̸ = 1 the loops in momentum space which connect a momentum in the central layer to itself through a path which visits both the top and bottom layers require at least m × n steps.The increased periodicity of these loops for m ̸ = 1 or n ̸ = 1 imply that the band structure in real space repeats itself for a displacement periodicity defined by the unit vectors: This result implies that the actual periodicity in real space is given by ℓ m2 /(mn).A sketch of the different regions in the unit cell, and the points where the band structure repeats itself is shown in Fig. S2.Densities of states (DOS) for two choices of the magic angle in the chiral limit, and {m, n} = {2, 1} are shown in Fig. S3.Similar results for {m, n} = {3, 2} are shown in Fig. S4.
The analysis discussed previously can be extended to multilayers with any number of layers greater than three.Results for the density of states in the chiral limit of an helical stack of four layers (θ 12 = θ 23 = θ 34 = θ 0 ) at different magic angles are shown in Fig. S5.
It is also instructive to generalize the alternative method discussed in [26].The analysis addresses the magic angles, in the chiral limit, of a stack of N layers, where each layer is twisted by the same angle, θ, and in the same direction, with respect to the previous layer.
The magic angles for N layers are determined by an N −component spinor which is the solution of N first order differential equations, similar to the case of a 2-component spinor, {ψ ⃗ k , ψ⃗ k }, which determines the magic angles in twisted bilayer graphene.The solution found in [26] describes the three component spinor for the trilayer in terms of three combinations of products of combinations of these wavefunctions: A similar equation can be written for an N −component spinor: We defined α 2 = [v F (4π)/(2d sin θ M /2)]/t AB , where θ M is the magical angle for a twisted bilayer.Then, Eq. (S13) has a solution for: We have checked that this equation gives the correct solution for α 3 = α 2 / √ 2 and α 4 = α 2 × 1/ 3 √ 6 agrees with our numerical results for a stack of four layers.For N → ∞ we find: where e is Euler's number.This result implies that a small shear force applied to the top and bottom layers of a graphite stack can lead to a peak at the density of states near the Fermi level.

VALIDITY OF THE LOCAL HAMILTONIAN EXPANSION
The moiré-of-moiré length scale emerges due to the misalignment of the moiré patterns between layers 1 and 2, and layers 2 and 3, (θ 12 + θ 23 )/2, as shown in the Fig. 1(c and we take the small angle approximation in the second line and commensurate angles (m ≥ 0 WLOG) in the third line.The moiré-of-moiré lengthscale is then l m2 = 2π/ √ 3k D min(P ab ), where we minimise P ab over all integers a ≥ 0 and b ≥ 0, except for the trivial solution P 00 = 0.For small commensurate angles with m and n small, this generally occurs at P |n|m = |mn(m+n)| 2 θ 2 0 , however, if m and n are large and of the same sign, P ab may be minimised at some other a ̸ = |n|, b ̸ = m.For example, at θ 12 = 1.0 • , θ 23 = 0.9 • , P 1,1 < P 9,10 .For completeness, we note that this description only considers lattice points lying on the two lines as shown in the sketch, when in fact there are others in the 2D plane.Provided the angles and thus the misalignment is small, the points on these lines remain the closest to each other.
Numerically calculated results for l m2 are shown in the Fig. 1(d) of the main text.As expected, l m2 diverges along the line θ 12 = −θ 23 , and is generally larger for θ 23 /θ 12 negative than for θ 23 /θ 12 positive due to the reduced misalignment of the moiré patterns.l m2 is also particularly pronounced at commensurate angles with simple ratios, where it is easy to minimise the first term of Eq. (S16) without the second term getting too large.For experimental samples smaller than l m2 , a particular local Hamiltonian expansion is expected to be a good description of the bands while for samples much larger than l m2 , it is necessary to average over expansion centers, or equivalently, over D 1 −D 3 .We note that the form for l m2 derived and calculated here agrees with that previously derived [37] for the special case of θ 23 /θ 12 ≈ 1.
Neglecting the slow periodicity in the local expansion, as is done in the numerical calculations of us and [37], neglects the misalignment and so forces the black lines of the upper panel sketch to coincide.In the case of commensurate angles, this leads to min(P ab ) = 0 and thus l m2 diverging, leading us to ignore the modulation of D 1 − D 3 .

FINDING THE MAGIC SURFACE IN TRILAYER CONFIGURATION SPACE
Our numerical calculations have shown that a magic angle pair may be found seemingly for any numerically accessible (commensurate) twist angle ratio, and so the natural question would be whether this may be extended to the continuum in angle-pair-space.We start by writing the trilayer Hamiltonian in layer space, with D l the Dirac cones of layer l and the H ij T as in Eq. (S7).The subscript L denotes layer space, and each "matrix element" here is itself a matrix in sublattice space.In the chirally symmetric limit, t AA = 0, a simple permutation allows us to rewrite the Hamiltonian in the dimensionless form where the subscript AB indicates sublattice space and the zeroes on the diagonals are the natural consequence of using the chirally symmetric model.The layer space submatrix is given by where ∂i = 1 2 (∂ x + i∂ y ) and their relative twists have been compensated for by corresponding countertwists of the spinors.U 12 (r) = j=1,2,3 e i 2π(j−1) D1) is the moiré potential between layers 1 and 2, and similarly for U 23 (r) replacing m with n and D 1 with D 3 .
The off-diagonal elements of the layer space submatrix have the form which may be collected into a single off-diagonal through the similarity transformation where the subscript L ′ reminds us of the nontrivial mixing of the layers.We are thus able to re-express the interfering double-moiré problem into a single effective moiré potential.We note that the transformation S is not unitary in general and so does not conserve the inner products of spinors, but does conserve the bandstructure.We further note that there is some freedom to choose the form of the single off-diagonal element; our choice most easily connects to previous work [19] where in the case of alternating twist angles, U 12 (r) = U 23 (−r) and so S is unitary and the form factor 1 + cd ab = √ 2, resulting in a constant scaling of the bilayer moiré potential and thus of the magic angle.Applying the transformation of Eq. (S21) to Eq. (S19) then allows us to obtain a single effective moiré potential, from which an estimate of the magic angle pair may be made by comparison with the bilayer moiré potential, a nontrivial extension of previous work [19] which focused on the special case of alternating twist angle multilayers.
Our estimate for the magic angle and the error of our estimate is calculated from the average and standard deviation of the form factor 1 + cd ab over space.While the form of the effective moiré potential we have chosen does indeed have singularities e.g. at the origin, these diverge as 1  x+iy and so the integral over the plane converges, which is unsurprising as it is possible to choose a form of effective moiré potential with no singularities at all, though in that case the analogy to twisted bilayer graphene is weaker.) for D = 0. Shaded regions denote the error in the magic angle pair prediction.The configurations studied in [26] and [19] are indicated with a black and purple diamond respectively, and those corresponding to our numerical calculations are marked with black diamonds.Dotted lines of constant θ12/θ23 are marked as guides to the eye.Bottom: Predicted range of magic angles over all D. The shaded region indicates pairs of angles for which there exists some D that would give flat bands.The configurations studied in [49][50][51] are marked in green, red and blue respectively.
Results are shown in the upper panel of Fig. S6 for zero shift, D = 0, where we indicate with a red line the predicted magic angle pairs according to the average-value analysis described above, with the uncertainty in the prediction shown by the shaded pink region.The alternating twist case previously studied [19] is marked with a blue diamond, and we note that this is the only point where the analysis is exact.We mark also our numerically calculated magic angle pairs with black diamonds, and note the excellent agreement with the more general analysis, giving us confidence that the results are reliable even for the numerically inaccessible case of incommensurate angles.Our analysis predicts two branches of magic angle pairs from which an infinite set descends, as in the bilayer case [18].At most two sets of numerically calculated magic angles per angle ratio are shown to connect with the two branches, though other values have been found.
We repeat the calculation over all moiré shifts to obtain a range of magic angle pairs, shown in the lower panel of Fig. S6.The magic angles vary continuously with shift and so the manifold of magic configurations is a threedimensional surface within the four-dimensional trilayer configuration space of two angles and a two-dimensional relative shift.This surface divides the configuration space into unconnected sectors, as seen in the lower panel of Fig. S6 where it is not possible to move from θ 12 = θ 23 = 0 to θ 12 = θ 23 = 2 • without passing through the magic surface.We remind the reader here that the plot is symmetric about the origin.This may have consequences on the topology of the bands as a function of trilayer configuration, as magic configurations may arise as the "midpoints" of band crossings, with the flat bands flattened by repulsion of the crossing bands.Previous work on alternating-twist magic angle trilayer graphene [49], with trilayer configuration (θ 12 = 1.57• , θ 23 = 1.57• , D = 0), a twisted monolayer-AA-bilayer system [50], configuration (θ 12 = 1.22 • , θ 23 = 0, D = 0), and a moiré quasiperiodic crystal [51], with configuration (θ 12 = 1.41 • , θ 23 = 1.88 • , D unknown), may be placed as single points within this extended magic phase.

TIGHT-BINDING RESULTS
We use a round disk method to construct the twisted graphene trilayers with arbitrary twist angles.The two independent twist angles θ 12 and θ 23 are chosen to be the rotation of the second layer 2 relative to the first layer 1 and the rotation of the third layer 3 relative to the second layer 2, respectively.The rotation origin is chosen at an atom site.We use a twist angel pair (θ 12 , θ 23 ) as the notation for different twist angle configurations.Positive (negative) values of the twist angle denotes counterclockwise (clockwise) rotations.The sample with (-θ, θ) has a mirror symmetry with the middle layer as the mirror plane.To calculate the property of these large scale systems with arbitrary twist angles, we construct the system in a large round disk.The radius of the disk should be set sufficiently large to rid the effects of edge states [44].In the actual calculation, the disk with radius of 172.2 nm (700d being d the lattice constant of graphene) and contains 10 million carbon atoms are utilized for the twist angles investigated in this work.text, there is a magic angle at θ 0 = 0.72 • .In the case of θ 12 = θ 0 , θ 23 = 2θ 0 , shown in Fig. S9, the first magic angle is θ 0 = 1.25 • , which is consistent with the result in the continuum limit in Ref. [52].
As discussed in the Ref. [25], in the case of AAA twisted trilayer graphene with θ 12 = θ 23 , the first magic angle decreases from 2 • to 1.57 • when the lattice relaxation is considered.A sharp peak is located at the charge neutrility point, shown in Fig. S10(a).The states of the peak in the real space are different for different displacement D are different, which is in agreement with the result in Fig. 2(a) of the main text.Such state feature can be detected by the scanning tunneling microscopy in experiment.However, For the relaxation of the system with the REBO and KC potentials, we could not find a strong reconstruction of the structure with the ABA region relax to large area and clear domain wall region, which may achieve via a stronger intralyer and interlayer potential [39].

DIELECTRIC SCREENING IN TRILAYERS
It has been widely shown that the dielectric response for multi-twisted materials increases drastically as the bands become flat [53].The plasmon screening on electron-electron interaction is believed to have a significant effect on the strong correlation phenomena.One potential advantage of magic multilayers over tBG is that the calculated noninteracting bandstructures are expected to be more resilient against Hartree effects due to the higher dielectric function from plasmon screening.
We may verify this by calculating the static component of the dielectric function.For simplicity, we ignore Umklapp processes, which would give a matrix structure to the polarization, with entries labelled by moiré reciprocal lattice vectors, ⃗ G, ⃗ G ′ .Within the framework of random phase approximation (RPA), the dynamic dielectric function can be calculated via where V (k) = 2πe 2 /κk is the Fourier component of the two-dimensional Coulomb interaction, with κ being the background dielectric constant.We choose κ = 4 in our calculation, in effect assuming the substrate is hexagonal boron nitride (hBN).Π(k, ω) is the dynamic polarization function given by Π(k, ω) = 2 q,ξ m,n f n q+k,ξ − f m q,ξ F nm q,q+k,ξ here m, n are band index, ξ is valley index, f m q is the Fermi-Dirac distribution, and F nm q,q+k,ξ = ψ † n,q+k,ξ ψ m,q,ξ 2 is the form factor of two different Bloch states.
Results are shown in Fig. S11, which shows ϵ(k, ω = 0) over the moiré Brillouin zone for alternating-twist magic angle trilayer graphene in the AAA and AAB configurations, and for tBG.We choose to compare with alternating twist trilayers to remove moiré-of-moiré effects, and to be able to make direct comparisons between the Brillouin zones of the two materials.The value of ϵ is seen to be D dependent, as expected from the D dependance of the bandstructure, and is noticeably higher than that of tBG everywhere except at Γ, indicating that there is a greater degree of plasmon screening in the trilayer, as expected, and that the calculated noninteracting bands are therefore expected to be more robust against Coulomb interaction effects than those of tBG.

ELECTRONIC BANDS AND BERRY CURVATURE.
The flat electronic bands at magic angles studied in the main text show finite Berry curvatures, and, generally, non zero Chern numbers, see also [39,41].Examples are given in given in

FIG. 1 .
FIG. 1. Schematic illustrating the trilayer parameters (a) twist angle pairs and (b) the moiré displacement D. The moiré unit cell of tBG with θ12 is outlined by a rhombus with alternating colors red and blue, and θ23 with colors red and green.The rhombus vertices lie at AA stacking regions of the respective layers.(c) Sketch of connected momentum states between layers with misaligned moiré patterns (orange and purple), resulting in a moiré-of-moiré modulation.Large angles θ12 = 4 • , θ23 = 3• are used for clarity, and less misalignment is expected for real magic angle samples.The lengthscale of this modulation is determined by the closest approach of second layer lattice sites formed from the two moiré patterns, indicated with a dotted line.The local expansion used in numerical calculation neglects this misalignment and, for commensurate angles, collapses this distance to zero.(d) Log-scale plot of the moiré-of-moiré length scale, lm2, in microns, as a function of twist angles.The left diagonal in white marks the line θ12 = −θ23, where there is no misalignment and lm2 diverges.lm2 also diverges when either angle is zero and there is no moiré-of-moiré pattern to speak of.

FIG. 2 .
FIG. 2. (a) Local DOS of a graphene trilayer in the chiral limit with θ12 = θ23 = θ0 = 0.72 • (magic angle for zero displacement) for three different displacements D between the top and bottom layers.The insert on the top panel are the three displacements, which are marked by a dot in the hexagon with the color corresponding to that of the local DOS.(b) As in (a), but for a combination of angles θ12 = 3θ0, θ23 = −4θ0, and θ0 = 0.404 • .Note the difference in the scales between the plots in (a) and (b).(c) Calculated magic angle pairs (θ12, θ23) for twisted trilayer graphene.The red and green colors are for magic angles with red and green displacements illustrated in insert of (a), respectively.The magic angle pair of (a) is marked by a blue star and that of (b) by a blue triangle.The blue square is the magic angle in twisted trilayer graphene with mirror symmetry.(d) Magic angle triplets (θ12, θ23, θ34) for stacks of four layers.Different colors are used for clarity.

FIG. 4 .
FIG. 4. The tight-binding results.(a) The density of states for a trilayer with angles θ12 = 1.41 • , θ23 = −1.88• .The gray line is the result in the chiral limit.(b) The density of states for a trilayer with angles θ12 = 2 • , θ23 = 2 • .(c) The local density of states at energy of 7 meV in the real space for trilayer in (a), with a zoom in of the states with displacements D marked by dashed square with colors.Red and gray correspond to the color in Fig. 2(a), and blue are in the middle region.
FIG. S1.Sketch of the momenta of the Bloch waves included in the continuum model, see eq.(S4), for m = 2, n = 1.The magenta arrows represent the vectors m⃗ gi, and the orange arrows represent the vectors n⃗ gi, see eq.(S6).

FIG. S2 .
FIG. S2.Sketch of the real space unit cell structure for (a) trilayer with {m, n} = {1, 2} and (c) trilayer with {m, n} = {3, −4}.The black hexagon defines the large unit cell, associated with moiré length ℓm2.The red hexagons show the moiré structure associated with angle θ0 and moiré length ℓm.Blue and green hexagons show moirés associated with angles θ12 and θ23.Note that the moiré lattices associated to θ0 and θ12 coincide.Special points in the momentum space unit cell with the AAA electronic structure for (b) {m, n} = {1, 2} and (d) {m, n} = {3, −4} cases.
FIG. S3. Results for the magic angles in the chiral limit and θ12 = 2θ0, θ23 = θ0.(a) -(d): The local DOS for a magic angle θ0 = 1.163 • with three displacements and their average, as in Fig. 3 of the main text.(e) -(h): As in the top panel, but for a magic angle θ0 = 1.14 • .

FIG. S6 .
FIG. S6.Top: Predicted magic angle pairs (θ12 , θ M 23 ) for D = 0. Shaded regions denote the error in the magic angle pair prediction.The configurations studied in[26] and[19] are indicated with a black and purple diamond respectively, and those corresponding to our numerical calculations are marked with black diamonds.Dotted lines of constant θ12/θ23 are marked as guides to the eye.Bottom: Predicted range of magic angles over all D. The shaded region indicates pairs of angles for which there exists some D that would give flat bands.The configurations studied in[49][50][51] are marked in green, red and blue respectively.
FIG. S7.Schematic of the AAA trilayer graphene with θ12 = θ23 = θ0 = 21.8 • .The θ12 and θ23 are the twist angles between layer 1 and layer 2, layer 2 and layer 3, respectively.The orange and purple rhombi are the moiré unit cells with θ12 and θ23, respectively, which has length ℓm = 0.65 nm.The black rhombus is the unit cell of the moiré-of-moiré with ℓm2 = 1.7 nm.The green square shows an ABA stacking region.
FIG. S10.The lattice relaxation effect.(a) The TB DOS of AAA trilayer graphene with magic twist pair θ12 = θ23 = θ0 = 1.57• .(b) The local DOS at zero energy in real space with a zoom-in of the states at different displacement D. The colors of the squares correspond to the ones in Fig. 2(a) of the main text.The hopping parameters are t0 = 2.8 eV and t1 = 0.44 eV.
FIG. S11.Static component of the dielectric function, ϵ(k, ω = 0) evaluated over the Brillouin zone for alternating-twist magic angle trilayer graphene at AAA and AAB (blue) stacking, compared to that of magic angle tBG (black).The value of ϵ for trilayers is greatly increased compared to that of tBG everywhere except the Γ point, indicating an overall greater degree of plasmon screening in trilayers and thus greater robustness of the calculated noninteracting bandstructures to Hartree effects.