Symmetries as the guiding principle for flattening bands of Dirac fermions

Since the discovery of magic-angle twisted bilayer graphene (TBG), flat bands in Dirac materials have become a prominent platform for realizing strong correlation effects in electronic systems. Here we show that the symmetry group protecting the Dirac cone in such materials determines whether a Dirac band may be flattened by the tuning of a small number of parameters. We devise a criterion that, given a symmetry group, allows for the calculation of the number of parameters required to make the Dirac velocity vanish. This criterion is employed to study band flattening in twisted bilayer graphene and in surface states of 3D topological insulators. Following this discussion, we identify the symmetries under which the vanishing of the Dirac velocity implies the emergence of perfectly-flat bands. Our analysis allows us to construct additional model Hamiltonians that display perfectly-flat bands at certain points in the space of parameters: the first is a toy model of two coupled 3D TI surfaces, and the second is a quasi-crystalline generalization of the chiral model of TBG.


I. INTRODUCTION
Since the discovery of superconductivity and correlated insulating states in magic-angle twisted bilayer graphene (TBG) [1][2][3][4] moiré materials have drawn tremendous attention as a tunable platform for creating novel electronic effects.The main feature of TBG is that by tuning the twist angle between the graphene layers one can tune the Dirac velocity at the Dirac cones to vanish to a remarkable degree of precision.The vanishing of the Dirac velocity is accompanied by a large density of states (DOS) at charge neutrality, thereby enhancing correlation effects.Following the example of magic-angle TBG, similar fine-tuned systems were shown to exhibit band flattening, with some examples being twisted trilayer graphene [5][6][7][8], twisted superconductors [9,10], and moiré patterns on the surfaces of 3D topological insulators (TI) [11][12][13].
The emerging plethora of flat-band Hamiltonians in fine-tuned materials raises the question of how generic this phenomenon is.In other words, what characterizes the set of systems for which fine-tuning a small set of parameters leads to the formation of flat bands or almost flat bands?This question is interesting both from the theoretical and a practical point of view, as a criterion for band flatness should be a useful guide in searching for new materials where exotic correlated phenomena may be found.
In this work, we focus on flat bands in systems harboring Dirac fermions, and more specifically on the conditions for flattening a band by making the Dirac velocity vanish.This scenario is relatively convenient to analyze theoretically, as it requires the knowledge of the Bloch Hamiltonian at only a single k point.The choice to focus on the Dirac velocities can also be motivated by noting that generically an upper bound to the bandwidth may be estimated by a Debye-like approximation to the band dispersion.More rigorously, in many cases of interest (see the examples discussed below) the quadratic order of the band dispersion near the k point vanishes by symmetry.In these cases, the vanishing of the Dirac velocity guarantees that the DOS diverges at least as (δE) −1/3 near the Dirac point.Notice that a quadratic band-touching cannot be obtained when the Dirac points are fixed by the symmetries of the system (for example by a rotation symmetry) as a result of the π Berry phase of the Dirac cone [14].
We begin in Sec.II by defining an algebraic criterion: We show that the symmetry group G acting on the Dirac cones determines the number of parameters (e.g., twist angle, pressure, etc...) that should generically be tuned to make the Dirac velocities vanish.This criterion is used to analyze different symmetry groups which can protect a Dirac cone, and to find the classes which allow for the tuning of a small number of parameters to obtain vanishing Dirac velocities.
In Sec.III we apply our criterion to the analysis of two systems of interest: The first one is band flattening of TBG, where we show that the existence of an approximate particle-hole symmetry is necessary for the vanishing of the Dirac velocity.The second system is surfacestates of 3D TIs under a periodic potential.We show for such systems that the Dirac velocity at charge neutrality can be made to vanish entirely by varying a C 2 symmetric potential.The resulting system might enable a platform for realizing strongly-interacting phases on the surface of a TI.
The vanishing of the Dirac velocity may be the first step towards a further increase of the DOS that culminates in a perfectly flat band [15][16][17].The first example of a model that exhibits such a band was a toy model of TBG [15].In Sec.IV we extend the ideas raised by [17] and our discussion of the vanishing Dirac velocity to discuss the symmetry requirements that are needed to obtain exactly-flat bands in general settings.We show that such flat-band Hamiltonians naturally arise in Dirac Hamiltonians with an external SU (2) gauge field by tun-ing a small number of parameters.For Hamiltonians in class CI of the Altland-Zirnbauer classification [18] we show that the flat-bands condition is equivalent to the vanishing of the Dirac velocity.For Hamiltonians with more general symmetries, we show that the flat bands can be found by considering the vanishing of the velocity in a modified version of the original Hamiltonian, which is in class CI.
We employ our discussion of exactly flat bands in Sec.V where we discuss two new model Hamiltonians which realize such exactly-flat bands, along with an in-depth analysis of a recently-proposed model.The first example is a continuum model with a C 4 symmetry, which can be thought of as a toy model of two TI surfaces with spin-flipping tunneling and an in-plane positiondependent magnetic field.The second example is a model of a quadratic band-touching Hamiltonian first proposed by [19], on which our analysis can be used to prove analytically the existence of exactly-flat bands.The last model is a quasi-crystalline generalization of the chiral TBG Hamiltonian.While the latter model does not have well-defined bands, we show it to host "magic angles" with an extensive degeneracy at charge neutrality.
Sec. VI concludes with a discussion of possible future directions.The appendices contain a more rigorous definition of our algebraic criterion, reviews of known results published elsewhere, technical proofs, and a discussion of edge cases that are not treated in the main text.

II. CONDITIONS FOR THE VANISHING OF THE DIRAC VELOCITY A. Zero-velocity co-dimension
Consider a two-dimensional Bloch Hamiltonian H(k) whose band structure has Dirac points for certain values of k.The velocity operators at the Dirac points are defined by The Dirac velocities (that is, the dispersion of the Dirac cone close to the Dirac point) are calculated using firstorder perturbation theory of v acting on the degenerate wavefunctions at the Dirac cone.They are the eigenvalues of the matrices where ψ m are the degenerate Bloch wavefunctions at the Dirac cone.We will use the notation ρ( Ô) to denote the projection of the operator Ô onto the subspace spanned by {ψ m }.Most commonly m = 1, 2, but we shall also consider the cases of n D degenerate Dirac cones, for which m = 1, ..., 2n D .We note that when the degenerate Dirac cones are protected by a local unitary symmetry (such as SU (2) spin rotation or a translation symmetry), we can consider each eigenspace of the symmetry separately as a single Dirac cone.
Our main interest is the condition for the Dirac velocities to vanish at some points in the space of parameters.We assume that H is controlled by a set of d parameters α 1 , ..., α d .We define the zero velocity codimension δ Z to be the codimension of the manifold in the space of α i for which ρ (v x ) = ρ(v y ) = 0. Roughly speaking, if the Hamiltonian H can have a vanishing Dirac velocity, δ Z is the number of parameters that should be tuned to make the velocity vanish.Our goal in this section is to show how δ Z can be calculated from the symmetries that preserve the Dirac cone.
Let G be the group of symmetries that preserve the Dirac cone.For such symmetries ρ(G) is a representation of G, for g ∈ G being a unitary operator.In the case where g is an antiunitary operator we obtain an antiunitary representation by multiplying the matrix ρ(g) mn by the complex-conjugation operator K. Since the elements g relate states |ψ i only to one another, they satisfy, for all g ∈ G.Note that lattice symmetries can relate v x and v y .The tuples (ρ(v x ), ρ(v y )) are therefore elements of the linear subspace V of tuples of 2n D -dimensional Hermitian matrices that satisfy (3).We have where f l (α 1 , ..., α d ) are real-valued and (M l,x , M l,y ) give a basis for V .Consequently, δ Z is the dimension of V .
Eq. ( 4) is then the statement of how the symmetries of the Dirac cones define δ Z .In the absence of symmetries protecting the Dirac cones, (4) means trivially that the matrices ρ(v i ) can be expanded in a basis of Hermitian matrices, and δ Z is the dimension of all possible tuples of such matrices, given by 8n 2 D (for example, with n D = 1, M i could be any of σ 0,x,y,z , giving a total of 8 tuples).
In the case of δ Z = 1, Eq. ( 4) reduces to ρ(v i ) are then fixed up to a real parameter and we obtain a vanishing Dirac velocity whenever f vanishes (see Fig. 1).In that case, we can make the Dirac velocity vanish exactly by tuning a single parameter.In Fig. 1 we show the trajectory for ρ(v x ) for δ Z = 2 and δ Z = 1 as we vary a single parameter.Noticeably, the Dirac velocity can vanish exactly as α is varied only when δ Z = 1.In general, given that we have d parameters and δ Z equations, the dimension of the zero-velocity solutions is d − δ Z .Notice that, since f in ( 5) is a continuous function, a change of sign in the Dirac velocity implies the existence of a point where it vanishes.Also, this sign is in agreement with the sign of v D obtained in perturbation-theory, e.g. in [15].A slightly more rigorous definition of δ Z is given in Appendix A. We show there that if the Dirac velocity is made to vanish at some point α 0 in the parameters space, there exists a manifold of dimension d − δ Z around α 0 in parameter space where the Dirac velocity remains zero.The proof relies on the assumption that the gap between the Dirac point and the rest of the bands does not close.Such closing of the gap results in the Dirac cone wavefunctions not being continuous and can create a boundary to the zero-velocity manifold.We treat an example of such gap closing in Appendix F.

B. Calculation of δZ for different symmetry groups
We now follow the principles outlined above to calculate δ Z for different symmetry groups G which preserve the Dirac point.The symmetry groups we choose to focus on may contain two antiunitary symmetries Θ, Π that anticommute (Θ) and commute (Π), with the operators v i , as well as their unitary product Σ, In cases where the system has local time-reversal T and particle-hole P symmetries that map the Dirac cone onto itself, they may serve as Θ and Π respectively.Such is the case for a Dirac cone on the surface of a three-dimensional TI.When T, P map between different Dirac cones, such as in the case of TBG, we can combine them with other unitary symmetries to form Θ, Π.We will identify these combinations when we discuss examples of the latter case.
In general, we do not demand that the symmetries are local.Furthermore, while we assume that the symmetries either commute or anticommute with the Hamiltonian, we define them only by the commutation relations (6) and not by their commutation/anticommutation relations with the Hamiltonian.When the symmetry group G exists, the symmetries constrain the possible representations of the two components of the velocity operator (v x , v y ).However, as long as the relations (6) do not distinguish between v x and v y , they are not enough to constrain δ Z to 1, since for any M l = (M x , M y ) that satisfies them Ml = (M y , −M x ) will do as well, leading to δ Z ≥ 2. To find cases for which δ Z = 1 we need an additional symmetry that acts differently on v x and v y .We, therefore, consider also a (unitary) reflection symmetry R which takes x → −x and thus satisfies Our main result in this section is a calculation of δ Z for all symmetry groups constructed from Θ, Π, Σ, R. Since Θ, Π are antiunitary, different symmetry groups are given by different choices of ξ T,P = ±1 defined by Besides the sign choice of ξ Π,Σ , different symmetry groups are distinguished by allowing R to either commute or anticommute with Θ, Π, Σ. Namely, for R fixed by R 2 = +1 we can choose Again ζ Θ/Π/Σ = ±1 and, assuming that all are present, ζ Σ = ζ Θ ζ Π .In the case where we have both Θ, Π we use the notation R ζΘζΠ to denote the commutation/anticommutation relations, while a single subscript will be used in the classes where we have only one of Π, Θ, Σ.The resulting family of symmetry groups is similar to that considered in the Altland-Zirnbauer classification with reflection symmetry [18,[20][21][22][23].
We calculate δ Z for n D = 1, 2. Table I presents the results for n D = 1.Table II presents the same case with more details and table III presents the results for n D = 2 (Tables II and III are presented in the Appendices).The tables are constructed as follows: Assuming that the Dirac cone is n D -times degenerate, for each symmetry group we construct a representation ρ.We then look for the possible representations of v x , v y which satisfy the commutation relations with the symmetry operators obtained from ( 6), (7).Finally, we assume the presence of a crystalline symmetry relating v x , v y such that δ Z is determined only by the dimension of possible ρ(v x ).Examples of such symmetry are C 3 and C 4 symmetries, where C n is a rotation of the system by 2π/n.A similar analysis can be straightforwardly extended to n D -fold degenerate Dirac cones for higher n D , include additional symmetries, or extend to three-dimensional Weyl and Dirac nodes [24].

III. APPLICATIONS A. Magic angles in TBG
As a first application of our analysis, let us consider the band flattening in TBG.Since the two Dirac points are fixed in different k points, we calculate δ Z for a single Dirac cone.Time reversal T maps between the two valleys of the electronic band.However, when multiplied by C 2 , we obtain an antiunitary symmetry that preserves the Dirac point and commutes with the velocity operators, such that it may serve as Π.Furthermore, the C 3 symmetry preserves the Dirac point as well (see a review of the TBG Hamiltonian and symmetries in Appendix B).We fix the representation of these symmetries on the two Dirac point wavefunctions to be which is dictated by the requirements that ρ (C 2 T ) should be anti-unitary, should square to +1, and should satisfy ρ(C 3 )ρ(C 2 T )ρ(C 3 ) −1 = ρ(C 2 T ).For the Dirac cone to be C 3 symmetric we must have {v x , v y } = 0 such that (4) becomes Evidently, these two symmetries are not sufficient to ensure that the Dirac velocity may be made to vanish with a variation of a single parameter α.Luckily, TBG at small twist angle θ has an additional approximate unitary particle-hole symmetry (broken by a term of order O(θ)) given by [25] C : η y σ x K. (12) where η i are the Pauli matrices acting on the layer indices.This symmetry can be combined with the exact symmetry C 2,x to form an additional symmetry that preserves the Dirac cone (see Appendix B for a review of the symmetries of TBG).Thus, under the approximation of a small twist angle the operator CC 2,x maps x → −x and anticommutes with the Hamiltonian at low energies.Consequently, it commutes with v x and anticommutes with v y .Choosing its representation to be σ x , it fixes f 2 = 0. Consequently, the magnitude of the Dirac velocity is given by |f 1 (α)|, which may be made zero when f 1 changes sign.This calculation leaves us with an important lesson: in TBG, both the exact and the approximate symmetries are necessary for the Dirac velocity to vanish at the magic angle.Indeed, by diagonalizing the Bistritzer-MacDonald (BM) Hamiltonian [1], we find that when one does not impose the approximate symmetry to be exact the Dirac velocity does not reach zero when θ is varied.It instead has a minimum value of ≈ 4 × 10 −4 × v 0 where v 0 is the Dirac velocity at zero coupling between the layers.The value of v D as the twist angle is varied, with and without the approximate C symmetry, is depicted in Fig. 2.
In magic-angle TBG the C symmetry is weakly broken by a symmetry-breaking term proportional to the small twist angle θ.When the twist angle is not small, for example when the layers are slightly twisted away from an angle of commensuration [26], the symmetrybreaking term has a non-trivial dependence on the twist angle.In particular, it does not vanish at commensuration angles.In a recent work [26], the authors describe the case of twisted graphene bilayers when the layers are twisted slightly away from a commensurate twist angle.The Hamiltonian obtained is similar to the BM Hamiltonian, but the C-symmetry is broken by a parameter that is independent of the deviation from the commensurate angle.As observed there, the Dirac velocity does not reach zero.Our analysis shows this to be a result of the breaking of C symmetry.

B. δZ in 3D TI surface-states
Refs. [11,12] suggested that the velocity characterizing the Dirac cone of the surface of a 3D TI may be suppressed by the application of a periodic potential on the surface.Here we use our analysis of δ Z to show that there exist "magic parameters" leading to an exact vanishing of the Dirac velocity.We present two types of periodic potentials that lead to a vanishing Dirac velocity: the first possesses a C 4 symmetry and requires tuning a  13) with a C4 symmetric potential in the P-symmetric and P broken case.The P-symmetric case is defined by the potential ( 16) with u = ux = uy while the P broken plot is given for the potential ( 17) with u1 = u, u2 = u 4 .In the latter case, the velocity is defined as the velocity of the Dirac cone connected adiabatically to the Dirac cone at zero energy as u is increased.One can see that vx can reach 0 when the P antisymmetry is broken, but not when it is present.
single parameter.The second has only a C 2 symmetry, so that each of v x , v y can be made to vanish by tuning a single parameter.The velocity vector can then be made to vanish entirely by tuning two parameters.
The Dirac cone on the surface of a 3D TI is protected by time-reversal symmetry T = σ y K.A periodic potential is consistent with this symmetry, leading to the Hamiltonian of the form where v 0 is the Dirac velocity at zero external potential and u(r) is a periodic potential term.Note that T allows only for a potential term proportional to the identity in (13) and prohibits the opening of a gap.
In the case where, in addition to T , there exist additional C n (n = 4, 6) and M x = σ y (x → −x) crystalline symmetries [27], we find from Table I that δ Z = 1 (since the Dirac point maps to itself under time reversal we have When a C 2 symmetry is present in the system, ρ(v i ) have no diagonal terms in the basis defined by |ψ , T |ψ , provided that |ψ is a Dirac-cone wavefunction that is a C 2 eigenfunction.The Dirac velocity can then be calculated by To make the discussion concrete, we start by taking u(r) of the form u(r) = 2u x cos q 0 x + 2u y cos q 0 y.
Generally, the Hamiltonian ( 13) with ( 16) is symmetric under Furthermore, the Hamiltonian anticommutes with the operator We now analyze the system both in the case where there is an additional C 4 symmetry and where this symmetry is broken.

C4 symmetric case
Our numerical studies of the C 4 -symmetric case indicate that the anticommutation of H and P prevents the vanishing of the Dirac velocity.Indeed, by diagonalizing the Hamiltonian we do not find any magic values (see Fig. 3) as u = u x = u y is varied.Besides the calculation presented in the figure, we checked that the Dirac velocity does not vanish up to u = 10.We also find numerically that the velocity does not vanish even when one considers additional C 4 and P preserving terms in u(r) with higher wave vectors.We therefore conjecture, but cannot prove generally, that a Hamiltonian of the form (13) with C 4 , M x , and P symmetries cannot yield a vanishing Dirac velocity for the Dirac cone at charge neutrality.
One can break P by introducing higher wave vectors in the potential u.As an example we take u(r) = 2u 1 (cos q 0 x + cos q 0 y) + 2u 2 (cos(q 0 x + q 0 y) + cos(q 0 x − q 0 y)).
When P is broken, the Dirac cone, which for u(r) = 0 is at E = 0 ceases to be fixed in energy.We can nevertheless calculate the velocity in the Dirac cone connected adiabatically to the one at E = 0 as the amplitude of u(r) increases.As an example, in Figure 3 we plot the Dirac velocity for u 1 = u, u 2 = u/4 and find points along the line in which the velocity vanishes.Note that while the first-order dispersion around the Dirac cone vanishes,  13),( 16) (in logarithmic scale).One sees that vx vanishes on lines in the ux, uy space, giving rise to a low energy Hamiltonian of the form (23).
we still have quadratic terms in the dispersion, leading to a finite (but increased) density of states at the Dirac cone.
In this example, breaking the P-symmetry opened a way for the Dirac velocity to vanish, despite the insensitivity of δ Z to this breaking.This may indicate that symmetries might also have an impeding role in the tuning of systems parameters to make the Dirac velocity vanish.Indeed, our analysis of δ Z gives necessary conditions, but not sufficient conditions, for making the Dirac velocity vanish by tuning a given number of parameters.

C4-broken case: vanishing of the velocity in a single direction
When the periodic potential breaks C 4 symmetry, the difficulties of making the Dirac velocity vanish are alleviated.In this case, each velocity component v x,y vanishes on a codimension-one manifold, which results in δ Z = 2. Indeed we find lines of "magic parameters" u x , u y which give one vanishing component of the velocity at charge neutrality (see Fig. 4).By tuning both parameters, we find points at which both components of the velocity vanish.
A simple, analytically-solvable example of this scenario can be found in the case where u y = 0 and the potential is one-dimensional.In that case, we can find zero-energy states of the form The Dirac velocities can be obtained from the represen-tations of the velocity operators in this basis.We obtain where J 0 is the Bessel function of the first kind.We notice the somewhat surprising result that it is v x , and not v y , which is independent of the potential [28], even though the potential varies along the x-direction.The matrix ρ(v y ) is proportional to σ y as a result of the inversion symmetry of u(x) and vanishes on the zeros of J 0 .We then find the "magic parameters" at u x /q 0 v 0 = .60,1.38, 2.16, ....A velocity that vanishes only in the y direction gives rise to a low-energy Hamiltonian of the form for some parameters ṽx , dx , dy .The DOS resulting from ( 23) vanishes as g(E) ∝ E 1/3 at low energies.An interesting question, which will not be elaborated on here, is the behavior of the Hamiltonian ( 23) when interactions are also considered.Renormalization-group analysis [29,30] suggests that in the presence of strong enough interactions, the dynamics of the system become quasione-dimensional, forming a Luttinger liquid phase in the x direction [31], possibly with spontaneous breaking of translation symmetry in the y-direction.
Since each velocity component vanishes on lines in (u x , u y ) space, we expect the entire velocity to vanish on points in that space.Fig. 5a shows that this is indeed the case.In Fig. 5b,c, we plot the band structure and DOS at one of these points.One can observe a peak in the DOS at charge neutrality, with two additional peaks corresponding to Van Hove singularities (vHs) at nonzero energies.Note that the functional dependency of the two peaks is different: while the DOS at the vHs diverges as − log(|δE|) (with δE being the deviation from the singularity), the DOS diverges as |δE| −1/3 around charge neutrality.Thus, the points in parameter space where the Dirac velocity vanishes might provide good candidates for strongly-interacting states at charge neutrality.Notice that the divergence is functionally similar to the higher vHs discussed in [12], but the low-energy Hamiltonian around the critical point is different [32,33].
The authors of [11] propose two methods for realizing a Hamiltonian of the form (16) on the surface of a TI, either by creating a moiré pattern on the surface or by posing a dielectric pattern on it [34].Here we note that the degree of tunability required to achieve the "magic coupling values" obtained in the C 2 -symmetric model can be achieved either by a C 2 symmetric dielectric pattern or by a potential generated by acoustic waves on the surface [35,36].Notice that the dimensionless parameter controlling the coupling strength is u/q 0 v 0 .Ideally, one would keep both u, q 0 high to mitigate disorder effects and to have a large range of momenta affected by the modulation.
Also, we note that in [13] the authors show that for two 3D TI surface states a spin-flipping tunneling term allows for the velocity to vanish.When only spinindependent tunneling between two such surfaces is allowed, the surface Hamiltonians can be written in fourcomponent spinors as when η i are the "surface" indices.The Hamiltonian H twisted-TI , therefore, splits into layer symmetric and antisymmetric sections, each described by the Hamiltonian (16).Our findings then provide two additional mechanisms for obtaining a Dirac cone with vanishing velocity in such systems.

IV. EXACTLY-FLAT BANDS
An intriguing aspect of the BM Hamiltonian is the presence of a limit [15] in which the Hamiltonian has an additional chiral symmetry and in which the bands at charge neutrality become perfectly flat at the magic angle.These bands may then be chosen to have non-zero Chern numbers and a well-defined sublattice polarization [37].The Hamiltonian at that limit, referred to as cTBG, allows for exactly flat bands to be reached by tuning a single parameter.In this section, we provide symmetry requirements under which the vanishing of the Dirac velocity implies that the bands are exactly flat.Our analysis provides conditions for small-codimension exactly-flat bands.
We begin by considering a generalized form of the cTBG Hamiltonian, described by a Dirac electron in a background SU (2) gauge field [38].That is, it is of the form where ∂ = 1 2 (∂ x + i∂ y ), Ā = A x − iA y with A being a non-abelian traceless gauge potential (here we focus mostly on the SU (2) gauge group).We assume that Ā is periodic on some lattice.The Hamiltonian H comes naturally with a chiral symmetry, which we choose to be Hermitian, S = σ z .We follow our notation for TBG and use σ i as the isospin indices and η i as the gauge indices.For any two solutions ψ a,b of Dψ = 0 (that is, Dirac cone wavefunctions in the same S indices) we can define the Wronskian for the functions ψ a,b as where the second index is a spinor index.In [17] the authors show that I(r) is position independent and further show that the condition for exactly flat bands is the existence of two orthogonal ψ a,b for which I(r) = 0 (see a review in Appendix D).This condition may be expressed as where W is the "Wronskian operator" defined by The Wronskian is antiunitary, commutes with S, and satisfies W 2 = −1 (the choice W = iη y σ z K satisfies the same requirements and gives an equivalent condition).Since the Wronskian inverts the direction of the spinor that it operates on, Eq. ( 27) implies that with ν(r) being a scalar.Let us consider the symmetry groups which allow a Dirac Hamiltonian as in (25).We first assume that the symmetries of the system keep ψ a,b orthogonal (this assumption is broken, e.g., in C 3 -broken cTBG, which we treat in Appendix E).Since we want the momentum operator to be diagonal in D, it must satisfy the condition If the Hamiltonian has additional time-reversal symmetry T , then T and P = ST preserve the space of the two degenerate Dirac cones and therefore satisfy the definition (6).Using ( 6) and (30) we then find that This requirement already restricts the possible AZ symmetry classes which can support a continuum Hamiltonian with exactly flat bands to AIII (S only), DIII (where , and CI (where T 2 = −P 2 = 1).Notice that the Hamiltonian ( 25) can be thought of as a surface Hamiltonian for class AIII, CI, or DIII topological superconductors, all of which can have protected Dirac cones on the surface [39,40].This proves that the Hamiltonian ( 25) cannot open a gap at zero energy.
We now treat each one of the above symmetry classes.We first consider the two time-reversal symmetric classes, DIII and CI.We find that class DIII Hamiltonians are too constrained to allow for the condition (26), while for class CI the condition can be fulfilled and is, in fact, equivalent to the vanishing of the Dirac velocity.We then show that the analysis of class AIII Hamiltonians can be mapped onto the analysis of class CI.Therefore, understanding the criteria for obtaining a flat band in the CI case is sufficient for the more general case of Hamiltonians of the form (25) with no time-reversal symmetry.

A. Class DIII
Here we prove that a 4 × 4 Hamiltonian of the form (25) of class DIII cannot support exactly flat bands.We begin by fixing v x = σ x , v y = σ y , S = σ z and T = σ y K (the choices T = η x σ y K and T = η z σ y K, where η are the gauge field indices, are equivalent).This restricts Ā to be of the form for some real a x (r), a y (r).Equivalently, we can write where a(r) = a x (r) + ia y (r).Since Ā commutes with itself at different positions the equation Dψ = 0 can be straightforwardly solved by integrating both sides.To do so we decompose a(r) as where the sum over G is a sum over reciprocal lattice vectors and G = G x + iG y .The zero modes of D can then be calculated explicitly as where f (z) is holomorphic.Since e ±u(r) is periodic (up to a phase) and therefore bounded, ψ ± is normalizable only for f (z) = const.We can conclude that for any a(r) there are only two zero modes for D, which are given by (35).Since these solutions do not satisfy (29), there are no exactly-flat bands for any choice of a(r).One can explicitly check that the Wronskian I(r) of ψ + and ψ − is constant and nowhere vanishes, since the spinors are never parallel.

B. Class CI
While class DIII symmetries limit Ā to the form (33), for class CI Ā does not necessarily commutes with itself at different points.As a result, the zero modes cannot be obtained by an integration procedure similar to (35).This allows for a richer structure of the zero modes and, most interestingly to us, allows for the vanishing of I(r).
Further notice that, for class CI Hamiltonians, the combination v x T satisfies and thus must be proportional to either η y K or σ z η y K.
We therefore have Since the RHS is an element of ρ(v x ) the vanishing of the Dirac velocity implies (27).This argument shows that in class CI Hamiltonians, as long as the solutions ψ a and ψ b are kept orthogonal, the vanishing of the Dirac velocity implies the existence of exactly-flat bands.From Table III we see that for class CI (that is, with Θ 2 = +1, Π 2 = −1) we have δ Z ≤ 2.

C. Class AIII
For class AIII we can write Ā in the general form Here W represents the U (1) part of the gauge potential (physically, ∇ × W is a magnetic field), while X, Y, Z are the three components of the SU (2) part.When W = 0 this is a system of class CI with T = η y σ y K. Let us first assume for simplicity that W (r) is periodic with mean zero (this represents a staggered magnetic field).We notice that, for a zero mode ψ of D, W (r) can be absorbed to ψ by defining where the operator ∂−1 is formally defined by with q = q x + iq y .Under this definition, we find that Dψ = 0 if and only if D ψ = 0.That is, we can reduce the problem of finding a zero mode for D to the case in which Ā is traceless.We conclude that for finding exactly-flat bands of ( 25) it is sufficient to solve for the time-reversal symmetric (CI) case, where, as we showed above, the vanishing of the Dirac velocity implies an exactly-flat band.
The flat bands created by this procedure have an exact correspondence with lowest-Landau-level (LLL) wavefunctions [17,41,42], and therefore have a nonzero Chern number for each S polarization (they can, in fact, be written in a form that resembles the LLL wavefunctions on the plane, see Appendix D).By considering the "squared Hamiltonian" H = H 2 , for which S acts as a local unitary symmetry, we see that the nonzero Chern number on each S index gives the middle bands a Z × Z topological index when T is absent.This index collapses to a (nonzero) Z index in the presence of T [43].Finally, note that cTBG is in class CI, as a result of an emergent intra-valley T symmetry [44].
It is important to distinguish between the exactly-flatbands models discussed here and the flat bands in tight binding models, e.g. in bipartite lattices [45][46][47][48] or in line-graph lattices [49,50].The models we discuss allow for the creation of exactly-flat bands by the tuning of a small number of parameters, assuming that a given set of symmetries are preserved.The small value of the codimension δ Z implies that even when symmetry-allowed terms give the flat band a dispersion, this dispersion can always be compensated by the lowest-momentum tunneling, and the flatness is recovered.This property is not there in the tight-binding examples.The bipartite lattice models have a flat band in all possible parameters, provided that the lattice remains bipartite.For the linegraph lattices, on the other hand, there is an infinite set of parameters that may be varied to destroy the band flatness while preserving the lattice symmetries, and can not be compensated by other parameters.A further difference is that the models discussed here are continuum, rather than tight binding, models.This property allows for the separation of the exactly-flat bands to bands of opposite Chern numbers by a symmetry-breaking perturbation, such as a sublattice potential in the case of cTBG.In lattice models, on the other hand, the existence of exactly-flat bands with nonzero Chern number is prohibited [51,52] (but such bands may carry fragile topology [48]).

V. EXAMPLES OF CLASS CI FLAT-BAND MODELS A. Chiral C4-symmetric model
The insights gained in the previous section can be used to construct a continuum Hamiltonian of class CI that can be tuned to have exactly flat bands at zero energy.Our model is manifestly distinct from cTBG in that it has a C 4 , instead of C 3 symmetry.We will therefore call it the C 4 symmetric flat-band model (C4FB).The presence of the C 4 symmetry has additional interesting implications, which will be discussed shortly.
The C4FB model consists of two Dirac cones on two topological insulator (TI) surfaces on the x − y plane, connected by z-reflection symmetry, and coupled via a spin-dependent tunneling term modulated by an in-plane magnetic field (a system with slightly similar features was analyzed in [53]).The two 3D TI surface-states are described by [54][55][56][57]: where σ = (σ x , σ y ).Each of these Dirac cones has a chiral symmetry S = σ z and a spinful time-reversal symmetry T = Kσ y , with T 2 = −1.Since we want the system to be in class CI, we need to preserve S and replace T by a symmetry T satisfying T 2 = +1 and {T , S} = 0. To that end, we introduce spin-flipping tunneling between the layers, whose phase is modulated by an in-plane magnetic field: where A z (r) is the vector potential associated with the magnetic field.The resulting symmetry T = Kη x σ x is a combination of T and the z-reflection R z = σ z η x .The Hamiltonian is then We additionally require the symmetries C 4 = e i π 4 σz (r → R 4 r), M x = σ y (x → −x) and a translation invariance by the unit cell.While any form of A z and u satisfying these requirements will give similar results, we choose for concreteness Notice that H can be made of the form (25) by a gauge transformation.
By diagonalizing H we find that, as expected from our calculation of δ Z (we have Θ = T , Σ = S, R = M x , see table III), the Dirac velocity vanishes on a codimensionone manifold in the u, φ space (see Fig. 6).Our discussion above shows that when the Dirac velocity vanishes the two degenerate σ z = +1 wavefunctions ψ a,b at k = 0 satisfy I(r) = 0. We can therefore write where ∂ν = 0 since both ψ a,b satisfy Dψ = 0.The function ν(z) is periodic on the lattice and C 4 symmetric, inherited from ψ a,b .Therefore ν (z) must have at least four poles per unit cell, located at four C 4 -related points.At these four points ψ 1 must be zero [58].Using the fact that ψ 1 has four zeros per unit cell we can construct four σ z = +1 linearly-independent zero-energy wavefunctions at each k, in the form [59] Λ k,n (z) = e ikxz i=1,...,4 where ϑ 1 (z | τ ) is the Jacobi theta function [60] (we use the convention defined in (D6)), z i = x i + iy i are the zeros of ψ 1 and w i satisfy where k = k x + ik y .That is, the positions of w i determine the momentum of ψ k , and the possible configurations of w i satisfying (50) give the four degenerate wavefunctions.This construction is similar to that of lowest-Landau-level wavefunctions on the torus [61].Note that our arguments did not rule out the possibility of having more than four wavefunctions per unit cell, but we expect four to be the general case.We give an example of ψ k=0 in Fig. 7.

B. Magic parameters in Hamiltonians with symmetry-protected quadratic band-touching
Here we consider a system, previously analyzed in [19], which displays symmetry-protected quadratic bandtouching (QBT).We show how the symmetry analysis can provide the condition in which a perfectly-flat band can be created.
The Hamiltonian we consider will be of the form: where u = u x + iu y .In [19] the authors suggest realizing this Hamiltonian by stacking two twisted layers of a material hosting a QBT point at the K-points.The Hamiltonian shown there realizes two copies of H with Below we provide a proof that H with u(r) given by ( 53) has flat bands with codimension 1 in α (this was shown numerically in [19]).

From a QBT to a Dirac Hamiltonian
As a first step, we show the relation between H and the Hamiltonians discussed in Section IV.We first note that H has a time-reversal symmetry T = σ x K and a chiral symmetry S = σ z , and is therefore in class CI, just as the flat-bands Hamiltonian discussed there.Furthermore, we can relate the 2-band Hamiltonian (51), which has second-derivative operators, to a 4-band class CI Hamiltonian of the form (25) having only first derivative operators, and having the same number of zero-energy states.
To do so, we notice that for any scalar wavefunction ψ satisfying Dψ = 0, we have where and v 0 = (2m 0 ) As a result, the four-component spinor 0, 0, ψ, v 0 ∂ψ is a zero-energy state of the Hamiltonian Notably, H inherits the class CI symmetries of H, which are given by T = η y σ y K, S = σ z (σ i , η i are the Pauli matrices in the spinor and gauge components, respectively).Similarly, it inherits any crystalline symmetry that H has.This shows that finding conditions for a flat band in Hamiltonians of the form ( 51) is equivalent to finding the conditions for a flat band in class CI Dirac Hamiltonians, which were analyzed in Section IV.There we showed that a flat band will be created as a result of a vanishing Dirac velocity, provided that D has two orthogonal zero modes.We notice that here, for u = 0, H has a QBT, which amounts to having a vanishing Dirac velocity (i.e. a vanishing expectation of the operator (1)) with only a single zero-mode of D. To conclude our argument, we show that when the quadratic term in the dispersion is made to vanish by the application of u, the QBT can be separated into two Dirac cones with vanishing velocity, resulting in perfectly-flat bands.

From a Dirac Hamiltonian to a perfectly-flat band
We now apply our analysis to provide conditions for the emergence of exactly-flat bands in H.We first notice that the QBT is stable when u is modified if and only if H is C 4 symmetric.When the QBT is unstable we get a similar situation to the one discussed in Appendix E: the model has two zero-energy Dirac points, whose momenta change as u is modified.As a result, the vanishing of the Dirac velocity will result in the Dirac points fusing to a QBT with nonzero quadratic dispersion, and exactlyflat bands will not generally form.To go on further we, therefore, need to assume that H is C 4 symmetric with C 4 = σ z (x → y, y → −x) (which is also the case for the model discussed in [19]).
With the assumption of a C 4 symmetry, we investigate which additional symmetries can help us make the bands perfectly flat.The form of H ensures that the dispersion remains quadratic around the band-touching point, so the projected Hamiltonian near the K point is restricted to the form Consequently, the quadratic dispersion of H vanishes with codimension 2. To reduce the codimension to 1 we require an additional symmetry that guarantees the vanishing of f 2 , which is a reflection symmetry.Acting on H, this additional symmetry is equivalent to the requirement that u is either purely real or purely imaginary.
Having shown that the quadratic term in the dispersion can be made to vanish, we now show that this vanishing leads to perfectly-flat bands.To do so, we map the problem to the problem studied previously, where we had two separate zero modes ψ a,b of D. This can be done by adding a small C 4 -breaking perturbation to u, of the form where v r (α) = f 1 (α) + if 2 (α) is the square-root of the renormalized quadratic dispersion around the K point and k ε is small.The resulting Hamiltonian has, therefore, two Dirac points at k = (±k ε , 0), which remain separate for any α.In addition, the Dirac velocity of the Dirac cones vanishes whenever α is tuned to make v r vanish.From the analysis of section IV, we conclude that H with the additional perturbation has exactly-flat bands whenever v r vanishes and, taking k ε to zero, we conclude that H and therefore H, has exactly-flat bands whenever v r vanishes.The analysis in this section shows, therefore, that Hamiltonians of the form (51) with a C 4 symmetry have exactly-flat bands with a small codimension (1 if there is a reflection symmetry present and 2 otherwise).We note that, while an example of such a Hamiltonian was first provided in [19] using two twisted layers of a 2D material with QBT points, the u-term in (51) can be obtained by a long-wavelength modification of the hopping in a single layer hosting a QBT point (for example by periodic strain from surface-acoustic waves).

C. Quasi-crystalline generalization of cTBG
Here we discuss a quasi-crystalline generalization of the cTBG Hamiltonian, namely the generalization of the C 3 -symmetric model to a C n symmetric model for odd n ≥ 3. We focus on the chiral case as it is the easiest to analyze theoretically.The family of Hamiltonians is given by the form (25) with with q j = (cos(2πj/n − π/2)), sin(2πj/n − π/2)).Clearly (59) reduces to the cTBG Hamiltonian [15] for n = 3.For n > 3 the model is not crystalline anymore, but nevertheless the formal analysis of magic angles in cTBG continues to hold.That is, we can calculate by perturbation theory in α the correction to the zero-energy wavefunctions.We get a zero-energy wavefunction of the form where the operator ∂−1 is defined in (42).Note that ∂−1 is undefined for q = 0, but the C n -symmetry prevents zero-momentum terms from appearing in the perturbation series.Another zero-energy wavefunction is given by acting on ψ K with the intra-valley C 2 symmetry While the Dirac velocity is no longer well-defined (as there are no Bloch wavefunctions), the Wronskian operator W given by (28) still is.Since the "Dirac cone" wavefunctions are still reflection symmetric, we have δ Z = 1 for the vanishing of the formal Dirac velocity given by When v D vanishes we find from equation (D4) that ψ K must have extensively-many zeros (that is, the number of zeros in a given area being proportional to the area.By repeating the analysis in Appendix D we, therefore, find a "band" with an extensive degeneracy of zero-energy wavefunctions. In Appendix G we describe the results of a perturbative calculation of v D similar to the one detailed in [15].For n = 5 we find the first "magic angle" at α 0 = .32.We plot the resulting wavefunction in Fig. 8.Each zero of the wavefunction can be used to construct a zero-energy state.If the system is confined to a finite size L, the number of zero energy states (up to corrections of order 1/L), will equal the total number of zeros of the magicangle wavefunction, and will be extensive with the system size.

VI. DISCUSSION
In this work, we discussed the symmetry structure that is required for the flattening of bands of Dirac fermions.We started with the requirement for the Dirac velocity to vanish and gave examples in TBG and a Dirac cone on the surface of a 3D topological insulator.Afterward, we discussed the vanishing of any dispersion, namely the symmetry requirements for the formation of exactly flat bands.We showed that, for a certain set of symmetries, the vanishing of the Dirac velocity implies that the band is exactly flat.
The symmetry considerations which allowed us to calculate δ Z do not provide us with a recipe for writing a Hamiltonian with δ Z parameters which can have a vanishing velocity, but generically suggest natural candidates for flat-band Hamiltonians.In one of the cases, which we studied in section III, the existence of an extra symmetry seemed to impede such vanishing by the most natural candidate Hamiltonian.This observation may indicate that there may be further symmetry considerations that may guide the search for such Hamiltonians.These are left here as a subject for future research.
While in this work we focused mainly on 2D moiré materials, much of our discussion can be straightforwardly extended to other systems and different tuning parameters.For example, one can consider 3D nodal line materials, where each k z slice can be viewed as a 2D subsystem, and k z can serve as an adiabatic parameter.Another interesting future question is the generalization of our results in section IV to the case of SU (N ) gauge fields.Studies of specific models, such as alternating-twist nlayer graphene [5] and chiral twisted graphene multilayers [62,63] show the richness of states that might arise in such cases.In the former, one can encounter exactly-flat bands coexisting with dispersive bands, while in the latter we see exactly flat bands with Chren numbers C > 2. It would be interesting to see whether it is possible to give a classification of the possible states in that case, in terms of the underlying symmetries.
Finally, another interesting direction is an experimental realization of the models we presented here.The flatband models discussed in section V are theoretically intriguing, but more work is needed if one wishes to find candidates for experimental systems which host them.On the other hand, we believe that the TI models discussed in section III can be realized using currently available experimental capabilities.Such an increase in the density of states on the surface of a TI could give rise to intrinsic superconductivity or correlated insulators [64] on the surface of a TI, or, more exotically, a gapped state that is symmetric to both time-reversal and charge conservation.Such a state must be topologically ordered with quasiparticles satisfying non-abelian statistics [65][66][67][68][69][70].
2. Here δ Z is an upper bound to the codimension of the zero-velocity manifold.Cases where δ Z is strictly larger than the codimension should arise in the case where there are additional low-energy emergent symmetries at the Dirac cones.An example can be given in the C 3 -broken cTBG Hamiltonian (eq.(E1) below): In this case, the exact Hamiltonian does not have a rotational symmetry relating v x and v y .On the other hand, the velocity operators satisfy (30), giving rise to an additional constraint on the codimension.
3. When the gap with the upper bands closes the Dirac velocity representations are no longer required to be continuous since ψ i are no longer continuous.A gap closing can therefore create a boundary (of dimension < d − δ Z ) to the zerovelocity manifold.We give an example of this scenario in Appendix F.
Proof.Since the gap between the degenerate point and the other bands does not close, we can calculate the correction to ρ(ô) for any operator ô via first-order perturbation theory, that is where for the velocity operator with the sum running only on ψ i outside the degenerate space.In the case of unitary operators g ∈ G which preserve the degenerate subspace we have so we find ∂ α ρ(g) = 0. We then get for v i that We also have where γ g i,j are real constants that depend on whether g commutes or anticommutes with H, as well as the transformation that g induces on k.Combining (A7) and (A8) gives For a given set of matrix representations ρ(g) for g ∈ G (A9) gives a set of linear equation on the tuples (ρ(v x ), ρ(v y )).The tuples satisfying (A9) then form a vector space whose dimension is δ Z (see the definition of δ Z in (4))).There must therefore be at least d − δ Z directions in α space in which the velocity doesn't change, giving a (local) zero-velocity manifold around α 0 whose dimension is at least d − δ Z .
Appendix B: Review of the Bistritzer-Macdonald model and symmetries Here we review the continuum model of twisted bilayer graphene (TBG) proposed by Bistritzer and Macdonald [1].

The TBG Hamiltonian
The Bistritzer-Macdonald (BM) Hamiltonian describes twisted bilayer graphene at small angles and low energies, at a single valley of the graphene layers.It is given by [1,71,72] where σ θ = e −iθσz/2 (σ x , σ y ) e iθσz/2 .The single-layer h are the Hamiltonians for a single Dirac cone in each graphene layer, twisted by a small angle.The tunneling matrices T i are with φ = 2π/3 and q 1 = k θ (0, −1) , q 2,3 = k θ ± √ 3, 1 /2.We have k θ = 2 sin (θ/2) k D ≈ θk D where k D = 4π 3 √ 3a0 and a 0 ≈ 1.4 Å is the distance between atoms in graphene.The scale w ≈ 110meV is the energy scale associated with the tunneling between the layers and the factor 0 ≤ κ ≤ 1 determines the ratio between AA and AB tunneling between the sublattices of the graphene layers.Real-world TBG has κ ≈ 0.7 as a result of lattice relaxation [73].
Important to some of our discussion is the chiral limit of TBG (cTBG) obtained by setting κ = 0.Under this assumption, we can remove the θ dependence in h(θ/2) by a gauge transformation.To write the resulting Hamiltonian in a form compatible with (25) we further re-scale the Hamiltonian by defining H = H/E 0 where E 0 = k θ w; define the dimensionless parameter α = w/k θ v; and rearrange the rows so that the Hamiltonian acts on the spinor (ψ 1 , ψ 2 , χ 1 , χ 2 ) (here the indices are layer indices, and ψ, χ live on the A, B sublattices, respectively).We obtain the chiral Hamiltonian [15,38] where z = x + iy, ∂ = 1 2 (∂ x + i∂ y ) and U (r) = e iq1•r + e iφ e −iq2•r + e −iφ e −iq2•r .

Symmetries
Let us discuss the symmetries of the BM Hamiltonian (B1).We define the Pauli matrices σ i , η i in sublattice and layer space, respectively.The point symmetries acting within the valley are given by [74] C 2 T : where K is the complex conjugation operator and R 3 is the rotation matrix by 2π/3.Of the three symmetries described above, only the first two preserve the Dirac points.
As a result of the small angle between the layers, the BM Hamiltonian has an additional approximate particlehole symmetry.If we take the approximation of setting θ = 0 in h(θ) the resulting Hamiltonian has a particlehole (PH) symmetry given by [25] C : η y σ x K. (B7) In the real-world model of TBG, the symmetry is broken in order O(θ).The combination CC 2,x gives an additional antisymmetry that preserves the Dirac cone, that is Finally, the Chiral model has, besides C, the additional chiral symmetry S : σ z . (B9) Since in the Chiral model the θ dependence in h is removed by a gauge transformation, the unitary PH symmetry C is exact here.We can therefore combine S and C to obtain an intra-valley unitary rotation symmetry that sends r → −r [44].By combining the intra-valley rotation with C 2 T we obtain the intra-valley time-reversal symmetry T = σ y η y K which satisfies (T ) 2 = +1.This shows that the cTBG model is indeed in class CI.Here we elaborate on our discussion of the Dirac cone on the surface of a 3D TI.In particular, we study additional parameters (besides the potential amplitude) which can be tuned to obtain a vanishing velocity for a C 2 -symmetric Dirac cone in a potential.The Hamiltonian of the form ( 13) and ( 16) is defined to be consistent with the T, M x , and C 2 symmetries, but is not the most general form consistent with these symmetries.More generally we can write an anisotropic form for the Dirac cone H = v x σ x p x + v y σ y p y + 2u x cos q x x + 2u y cos q y y. (C1) Here v x /v y can be controlled by applying strain on the TI while q x /q y can be controlled (for example) by an asymmetry in the dielectric pattern.By rescaling the y axis we can make v y = v x = v 0 .We, therefore, write the Hamiltonian H = v 0 p • σ + 2u(cos q 0 x + β u cos q 0 β q y).
(C2) which is controlled by the dimensionless parameters β u , β q , u/q 0 v 0 (the first two define the C 4 -symmetry breaking).In Fig. 9 we plot the velocity of the Dirac cone of (C2) at charge neutrality as a function of u, β q .The results show similar "magic parameters" to the case discussed in the main text (see Fig. 5a).
Appendix D: The Wronskian operator and requirements for exactly flat bands In this section, we restate some results from [17] that are useful to our discussion of the condition of exactly flat bands in chiral-symmetric continuum models.We begin with a Hamiltonian of the form where ∂ = 1 2 (∂ x + i∂ y ), Ā = A x − iA y with A being an SU (2) gauge potential.We assume that H is symmetric under translations by the lattice vectors a 1 , a 2 .Given two solutions ψ a , ψ b of the zero mode equation one can write the Wronskian Importantly, we find that I(r) = const.This is because When there is no external magnetic field we have tr Ā = 0, so I(r) = I(z).However, since I(z) has no singularity it must be constant.Next, we find that when I = 0 there must be an exactly flat band at zero energy.This is because, in that case, we must have We have, however, ∂ν(r) = 0, as can be seen by applying D on both sides of (D4).We can therefore write ν(r) = ν(z).Assuming that ψ a,b are orthogonal (which must be the case if they are positioned on different points in the BZ) ν(z) is a non-constant meromorphic function.It therefore must have a pole at some point z 0 in the unit cell.At this point, ψ a must have a zero.We can then construct additional wavefunctions in the flat band by writing where a i = a i,x +a i,y , with a i being lattice vectors.Here ϑ 1 (z | τ ) is the Jacobi theta function, defined by (−1) n e iπ(n+1/2) 2 τ sin (2π(n + 1/2)z).
(D6) Importantly, the pole of the ϑ 1 cancels the zero at ψ a making ψ k as defined above normalizable.
One can follow an alternative approach for the construction of the flat-band wavefunctions, which makes the similarity between the flat-band wavefunctions and the lowest Landau levels manifest.We can choose a basis of such functions in the form [42] where f (z) is any holomorphic function, A is the unit cell area, and G(r) is a structure function that captures the lattice dependency of the wavefunction.It is given by Interestingly, it can be checked that |G(r)| is periodic with the lattice.similar argument shows that in the presence of a weaker C 2 symmetry we must have at least 4 flat bands (2 per S eigenvalue).Interestingly, we find (Fig. 12a) that the magic angle obtained by tuning α at β 1,2 = 0 is unstable when the C 4 symmetry is broken by a nonzero β 1,2 .That is, taking β 2 = 0, β 1 = 0, for example, the codimension of the zero Dirac velocity manifold in the (α, β 1 ) space is not δ Z = 1 as can naively be expected from Table III.Rather, the velocity vanishes on a point in (α, β 1 ) space (where This apparent contradiction is resolved by noting that once either β 1 or β 2 are non-zero, the conditions of the theorem we proved in Appendix A are not satisfied, and Table III cannot be used to infer the codimension.Namely, the theorem requires that there is a gap between the Dirac cone and the higher bands.In the case discussed here, for β 1,2 = 0 the C 4 symmetry requires that there are 8 degenerate zero-energy bands at the magic angle.Since the Dirac cone is only 4-fold degenerate, We therefore find that there are 8 bands connected to E = 0, all of which become exactly flat when the velocity vanishes. there must be 4 additional states closing the band gap (see Figure 13).In the presence of the C 4 symmetry the additional bands do not hybridize with the Dirac point wavefunctions as they have different C 4 -eigenvalues, and we can still use the results of Appendix A. On the other hand, C 4 breaking couples the Dirac cones and higherband wavefunctions, breaking the assumptions made in the theorem.
A different way of understanding this phenomenon is that the definition of f as in ( 5) is by adiabatic continuation: we always require that the wavefunctions are changed continuously as we vary the control parameters.Since the Dirac points wavefunctions are changed discontinuously at the magic angle f cannot be defined to be both continuous and single-valued.As an example, in Fig. 12b we draw f (α, β 1 ) going in two trajectories, one with β 1 = 0 and the other with nonzero β 1 .The sign obtained for f (α, β 1 ) is opposite as a result of the discontinuity.

Appendix G: Perturbation theory for the
Cn-symmetric quasi-crystalline models Here we calculate in perturbation theory the formal Dirac velocity for the quasi-crystalline models (59).We can numerically calculate higher-orders in the perturbation series for n = 5.We find The codimension δZ of the zero-Dirac velocity manifold for a single non-degenerate Dirac cone, according to the symmetry group, with the symmetries Θ, Π, Σ, R defined in (6).In the column of each symmetry a zero denotes the absence of the symmetry, while the sign denotes the square of the symmetry.The signs of the reflection operator R ζ Θ ,ζ Π reflect the commutation relations of R with Θ, Π.We omitted the rows that could not give rise to a single Dirac cone with the listed symmetries.
For each symmetry group we specify a representation for the Θ, Π, Σ, R operators and write matrices spanning the linear space of possible vx representations that satisfy ( 6) and (7).δZ is then the dimension of this linear space.

FIG. 2 .C 4 FIG. 3 .
FIG.2.The Dirac velocity vD of the BM model (B1), with and without the additional C symmetry, as the twist angle θ is varied near the magic angle.The C symmetry is imposed by setting θ = 0 in the h(θ) terms.

FIG. 4 .
FIG.4.The x component vx of the velocity for the Dirac cone at charge neutrality for the Hamiltonian (13),(16) (in logarithmic scale).One sees that vx vanishes on lines in the ux, uy space, giving rise to a low energy Hamiltonian of the form(23).

FIG. 5 .
FIG. 5. (a) The "absolute" Dirac velocity v 2 x + v 2 y of the Dirac cone at charge neutrality for the Hamiltonian (13),(16) as a function of ux, uy, in logarithmic scale.The dark spots are codimension-two manifolds on which the Dirac velocity vanishes entirely.(b,c) Band structure and DOS of the Hamiltonian (13),(16) at the "magic angle" obtained with ux = 0.95, uy = 0.25 (the point is marked by a circle in (a)).We find a divergent DOS at charge neutrality with additional Van Hove singularities at E = ±0.21.

FIG. 6 .FIG. 7 .
FIG. 6. Dirac velocity for the Dirac cone of the Hamiltonian (45) as a function of the parameters u, φ, in logarithmic scale.The Dirac velocity vanishes exactly on the dark line, leading to an exact flattening of the bands.

FIG. 8 .
FIG. 8.The wavefunction ψK in the chiral C5 symmetric model (59) at the magic angle.The orange points signify zeros of the wavefunction.
Appendix C: Additional parameters for tuning a C2-symmetric vanishing-velocity Dirac cone

FIG. 9 .
FIG. 9.The "absolute Dirac velocity" of the Dirac cone at charge neutrality of the Hamiltonian (C2) (similar to Fig 5a).The dark valleys are points of vanishing velocity.

FIG. 10 . 4 FIG. 11 .
FIG. 10.(a) The trajectory of the Dirac points in a C3 broken cTBG near the magic angle when α is varied.For a small symmetry-breaking parameter β the Dirac cones remain close to the K, K points away from the magic angle.Near the magic angle, the displacement δkD diverges and the Dirac cones travel around the BZ, meeting twice to form QBT points at the Γ and M -points.(b) The normalized Dirac velocity at the Dirac cones for C3-symmetric (β = 0) and weakly C3 broken (β = 0.05) cTBG Hamiltonian (E1).In the C3 broken case, the Dirac velocity vanishes twice, at the two QBT points.

FIG. 12 . 4 EFIG. 13 .
FIG. 12. (a) The value of the Dirac velocity at the Dirac cone for the Hamiltonian (F2) with a broken C4 symmetry.The codimension of the zero Dirac velocity manifold is 2 (the dark point).(b) The values of the "signed" Dirac velocity f (as defined in (5)) going through two different trajectories.Since the wavefunctions are discontinuous at the magic angle the velocities acquire an opposite sign.

TABLE III .
δZ for two degenerate Dirac cones: Same as Table II but for doubly-degenerate Dirac cones.The dashed rows signify symmetry groups which cannot support the algebra of a Dirac cone.