Higher-order van Hove singularity in magic-angle twisted trilayer graphene

We study the presence of higher-order van Hove singularities in mirror-symmetric twisted trilayer graphene. This geometry has recently emerged experimentally as a fascinating playground for studying correlated and exotic superconducting phases. We find that the trilayer hosts a zero-energy higher-order van Hove singularity with an exponent -1/3. The singularity is protected by the threefold rotation symmetry and a combined mirror-particle-hole symmetry and it can be tuned with only the twist angle and a perpendicular electric field. It arises from the combined merging of van Hove singularities and Dirac cones at zero energy, beyond the recent classifications of van Hove singularities. Moreover, we find that varying a third parameter such as corrugation brings the system to a topological Lifshitz transition, with anomalous exponent -2/5, separating regions of locally open and closed semiclassical orbits.

Introduction.-Strong correlations generally result from a quenching of electronic motion, comparatively magnifying the strength of electron-electron interactions. This is the case for engineered flat bands, such as for fractional quantum Hall states, but also more generally in proximity to a singularity in the density of states. Tuning the chemical potential at a van Hove singularity (VHS) [1,2] introduces a large number of singleparticle states with negligible energy likely to form a correlated state. Conventional van Hove singularities entail a logarithmic singularity but there are also higherorder types [3,4] with more diverging power law scaling which have been recently classified in single-band electron models [5,6]. Such strong divergence amplifies correlation effects and plays a key role in determining the ordering instabilities in various materials, such as twisted bilayer graphene [4,[7][8][9][10][11], biased Bernal stacked bilayer graphene [3], twisted bilayer transition metal dichalcogenides [12], Sr3Ru2O7 [13], heavy fermions materials [14], and high-T c superconductors [15]. Moreover, it has been recently shown that the higher-order van Hove singularity gives rise to a novel non-Fermi liquid critical state dubbed supermetal [16].
In this letter, we argue that a strong higher-order van Hove singularity emerges in the single-particle spectrum of TTG upon tuning the displacement field and rotation angle. It results from the symmetric merging at zero energy of two standard VHS with opposite energies. Located at the band touching K point of the Moiré Brillouin zone, it falls outside the single-band classification of VHS performed in Refs. [5,6], and also differs from the higher-order VHS identified [4] in TBG. It exhibits the power law scaling ω −1/3 , stronger than ω −1/4 predicted for TBG. An even stronger exponent −2/5 is found by tuning the corrugation, indicating a Lifshitz transition between two topologically incompatible energy contours.
Model and mapping to TBG.-The starting point is the continuum model [17,18] for trilayer graphene, where the three layers are stacked with alternating twist angles ±θ [64], coupling the three Dirac cones in each valley. It is convenient to take advantage of the mirror symmetry M z with respect to the middle layer in the absence of displacement field D/ε 0 , and write the Hamiltonian in an already layer-rotated basis [64,67,69,70] wherek = −i∇ r , v 0 is the electrons velocity in graphene, T (r) = 3 j=1 e iqj ·r T j are the interlayer hoppings with sublattice structure, q 1 = k θ (0, 1), q 2/3 = k θ (∓ √ 3/2, −1/2), where k θ = |K|θ is the distance between the K points of consecutive layers. ±U are the gate potentials applied on the top and bottom layers, 3 nm is the interlayer distance and 0 is the bare dielectric constant. For simplicity we neglect screening of the displacement and take ε = 1. The interlayer tunneling matrices take the form T j+1 = w AA σ 0 + w AB σ + e −2iπj/3 + σ − e 2iπj/3 , where j = 0, 1, 2, w AB = w and w AA = w r with r < 1 due to lattice relaxation effects [83,84], and the σ 0,± matrices act in sublattice space. Hereafter we take r = 0.8 unless stated otherwise. As Eq. (1) governs the valley K, the valley K is simply obtained by time-reversal symmetry, i.e. complex conjugation.
At zero displacement U = 0, Eq. (1) decouples a highvelocity Dirac cone in the odd mirror sector, located at the K point of the Moiré Brillouin zone, from two coupled Dirac cones in the even mirror sector, one at K and the other one at K . The even sector maps exactly onto the TBG Hamiltonian with a rescaled interlayer tunneling [64], and a corresponding √ 2 enhancement of the magic angle θ m 1.541 • at which the Dirac velocity vanishes. Computing the density of states, we recover the symmetric VHS of TBG in the two active bands, see Figs.1 (a) and (b), observed in STS experiments [85][86][87][88][89]. The VHS is in fact composed by three equivalent saddle points related by C 3z symmetry, i.e. 2π/3 rotation in the graphene plane. At a critical angle θ v 1.571 • , a higher-order VHS [4] occurs below which each saddle point splits in two. The DOS singularity at θ v is a power law with exponent −1/4 and asymmetry ratio ∼ √ 2, which convincingly matches the experimentally observed VHS peak [88]. The C 2x symmetry denotes a π rotation around x exchanging layer and sublattice indices. It enforces the saddle points to be on the ΓM lines of the Moiré Brillouin zone for θ larger than θ v . Throughout the rest of this paper, we will focus on angles θ such that the higher-order VHS inherited from TBG does not play a role. Aside from the symmetries M z , C 3z , C 2x already introduced, the model in Eq. (1) is also invariant under the combination C 2z T . It anticommutes with the particle-hole symmetry P in the even mirror sector and with C 2x P in both odd and even subspaces.
VHS merging and effective model.-We discuss the case of non-zero displacement field U = 0 and explore the evolution of the density of states at the twist-angle θ = 1.59 • . U breaks the mirror symmetry M z , the rotation symmetry C 2x and the particle-hole symmetry P . However, the product M z C 2x P is preserved [70] together with the remaining symmetries. M z C 2x P acts as anticommuting particle-hole symmetry. It will be very important in stabilizing the new van Hove singularity, see below. Fig. 1 (a) shows the density of states as a function of the number of electrons per Moiré unit cell and electric displacement field D/ε 0 . The two particle-hole symmetric VHS, represented as red dots in Fig. 1 (a), move towards charge neutrality as D increases and merge at zero energy at a critical D c /ε 0 0.37 V/nm. The VHS peaks become concomitantly more pronounced with increasing D as shown in Fig. 1 (b), moving towards what seems to be a zero-energy divergence. The corresponding evolution of the C 3z symmetric saddle points with the displacement field is shown in Fig. 1 (c). As a consequence of C 2x symmetry breaking, the saddle points leave the ΓM lines and converge towards the K points in the Moiré Brillouin zone. By further increasing D above D c the VHS split again and move away from charge neutrality together with a substantial reduction of the VHS peaks.
In order to gain more analytical insight into the zeroenergy merging of VHS and the marked singularity in the density of states, we derive a low-energy approach close to the K point in the Moiré Brillouin zone. Being a high-symmetry point, K retains some of the symmetries of the model that leaves it invariant: C 3z , C 2z T and M z C 2x P . C 2z T protects a Dirac cone at K even at non-zero displacement field. The two degenerate states at K, denoted u ω and u ω * , can be classified by their C 3z eigenvalues ω = e 2iπ/3 and ω * , respectively. M z C 2x P pins these two states at zero energy and more generally enforces a fully particle-hole symmetric spectrum at K. The form of the low-energy Hamiltonian in the vicinity of K is constrained by the symmetries. In the basis of the two degenerate states (u ω , u ω * ) defining the Pauli matrices τ 0,x,y,z , it takes the C 3z -symmetric form to third order in q = k − K, where q 2 = q 2 x + q 2 y and the couplings {v, η, γ, ξ} are calculated by employing a k · p approach with degenerate perturbation theory (see below). The symmetry M z C 2x P , acting as , forbids other C 3zsymmetric terms, such as q 2 τ 0 , (q 2 x − q 2 y )τ x − 2q x q y τ y or (q 3 y − 3q y q 2 x )τ 0 , to appear in Eq. (2). The crucial reduction to four coupling constants {v, η, γ, ξ} implies that the higher-order VHS, taking place when v and η both vanish, requires the fine-tuning of only two parameters of the original model, such as the twist angle θ and the displacement field D.
Beyond the zero-energy subspace (u ω , u ω * ), it is convenient to introduce the other eigenstates of Eq. (1) at K, | u n =| u nK with energies n , and the operator P = − n | u n u n | / n . The values of the coupling constants are calculated from a k·p approach at K [90] and depend on the details of the spectrum of the trilayer. They are shown in Fig. 1 (d) as a function of the displacement field for the angle θ = 1.59 • measured in [80]. Remarkably, v and η vanish for values of D in very close proximity, suggesting the vicinity to the higher-order VHS. It explains the strong feature observed Fig. 1 in the density of states although, rigorously speaking, reaching the higher-order VHS requires the fine-tuning of an additional parameter. This can be done by changing the twist angle θ. In order to get a closer look at the exact position of the higher-order VHS, we trace out the values of D and θ for which the velocity and the curvature of the band dispersion both vanish at the K point. The result as function of the corrugation parameter r is shown in Fig. 2 (a). As anticipated, we find a line of higher-order VHS in this tridimensional parameter space. We stress again that, under the reasonable assumption of a corrugation parameter r = 0.8 and for the twist angle θ = 1.59 • , one gets very close to the higher-order VHS by simply tuning the electric displacement. The trajectory of higher-order VHS in Fig. 2 (a) originates from a critical angle θ CH 1.536 • where both the displacement field D and the corrugation r are vanishing. This point corresponds in fact to the chiral limit of TBG [91][92][93][94][95] where the whole active band is rigorously flat. At r = 0, θ = θ CH is also the first magic angle. As D increases, the line of higher-order VHS extends till a critical value of the atomic corrugation r 0.842 above which it is not possible to have v and η both vanishing. At smaller twist angles, we find that each magic angle θ (n>1) CH in the chiral limit r = 0 is the starting point of a similar line of higher-order VHS.
The higher-order van Hove singularity.-The higherorder VHS appears for a vanishing velocity v = 0 and curvature η = 0 in Eq. (2). In this case, the dimensional scaling H K (λ 1/3 q) = λ H K (q) indicates a powerlaw scaling close to charge neutrality, with exponent −1/3. The Hamiltonian describes a cubic band touching and therefore extends the one-band classification of VHS [5,6]. We would recover effectively one band by setting ξ = 0 and having H K (q) = 2γ(q 3 x − 3q x q 2 y )τ 0 . In the vicinity of the higher-order VHS, there are three C 3z -symmetric saddle points at positive energy, three at negative, all merging at K as v and η are tuned to zero, as shown in Fig. 1 (c), reminiscent of the saddle point merging in the single band case [3]. In addition, the effective two-band structure entails a pseudospin with a vorticity of +1 at K protected by C 2z T . In the vicinity of the higher-order VHS, K hosts a +1 Dirac cone surrounded by six additional Dirac cones [90]. The side cones are organized in two C 3z -symmetric triplets along K K (ΓK) with vortic-ity +1 (−1). The higher-order VHS occurs precisely as all Dirac cones meet at K, leading to cubic band touching with +1 vorticity. A single Dirac cone at K remains above the critical value D > D c .
In the one-band limit of vanishing ξ, the semiclassical orbits are fully open with an anisotropic elliptic umbilic (D − 4 ) structure [5,6] separated by C 3z -symmetric lines and a cusp at q = 0. Quite on the contrary, γ = 0 predicts closed and isotropic semiclassical orbits. The corresponding iso-energy contours are illustrated in Fig. 3 (a-c). This difference in topology indicates that a Lifshitz transition is expected to occur as function of the ratio 2γ/ξ. This is visible in the divergence in the prefactor F of Eq. (4) shown in Fig.2 (b). Indeed, for ξ = 2γ, the spectrum at the higher-order VHS is E ± (q) = ξq 3 [±1 + cos(3φ)], where φ is the polar angle between q and the x-axis. It predicts zero-energy lines along the axis at φ = ±π/3, π[π], shown in Fig. 3 (b), which are responsible for an even stronger divergence in the density of states. The zero-energy lines are lifted by further expanding the k · p approach in Eq. (2) to the first non-vanishing order. The fourth order being zero, due to the symmetries C 3z and M z C 2x P , the dispersion finally takes the form [90] E + (q) = ξq 3 [1 + cos(3φ)] + Λq 5 (5) in the vicinity of the critical line at φ (the other lines are inferred by symmetry), where Λ is a positive coefficient. Thanks to this term, the divergence in the density of states is cut off at the angle δφ ∼ Λ/ξ(ω/Λ) 1/5 , where δφ = φ − π, and a new power law is obtained stronger than the naive scaling dimension −1/3 of the higher-order VHS. Coming back to the original trilayer model in Eq. (1), the ratio 2γ/ξ along the (red) line of higher-order VHS is represented in Fig. 2 (c). We then find that the critical point ξ = 2γ is reached when the corrugation parameter is r 0.34, corresponding to the twist angle θ 1.543 • and the displacement field D/ε 0 0.086 V/nm. At this point, the DOS critical exponent thus changes to −2/5 while it is −1/3 along the rest of the line, indicating a change of topology of the iso-energy contours, from local open orbits when r < 0.34, to closed ones when r > 0.34.
The change of Fermi surface topology is reflected in the Landau level (LL) spectrum [96] calculated in Fig. 3 (d) for 0.01 T. It is obtained from minimal coupling applied to the effective model in Eq. (2), at vanishing (v, η), where the q 4 terms have been added to regularize the region with open orbits. We find markedly different behavior of both sides of the Lifshitz transition: a regularly spaced structure for closed orbits, corresponding to 2γ/ξ < 1 (and r > 0.34), which evolves into weak oscillations in the regime of open orbits. The LL energy scaling ∝ (nB) 3/2 [3] of the higher-order VHS, with LL index n and magnetic field B, becomes ∝ (nB) 5/3 at the transition due the anomalous −2/5 DOS exponent. This yields in Fig. 3 (d) a relative collapse of LL at the transition, even more pronounced for decreasing B.
Conclusions.-We derived the existence of a higherorder VHS in mirror-symmetric TTG associated with a strong zero-energy peak in the density of states. In contrast with twisted bilayer graphene, the HOVHS arises from the fusion between the standard finite energy VHS and the Dirac cone at K, and it is protected by C 3z and M z C 2x P . As long as these two symmetries are not broken, the HOVHS will still be present with effective parameters which are renormalized by the electron-electron interaction. We emphasize that the peak in the DOS is strongly enhanced in the vicinity of the HOVHS, as shown in Fig. 1(a), and that the position shifts toward charge neutrality. This has two consequences: there will be a larger range of doping for which the Stoner criterion is fulfilled, and this criterion is satisfied for smaller doping to charge neutrality. We argued that current experiments at the magic angle can be brought close to this strong singularity by electric gating, although very strong fields might be needed to compensate screening. This finding has far-reaching consequences, as this strong divergence in the density of states, will result in the emergence of a plethora of many-body phenomena that have barely begun to be studied [75]. Moreover, the type of van Hove singularity that we found differs fundamentally from the standard van Hove singularity, such as the one classified in Refs. [5,6], as it also involves a Dirac point. This new structure gives a topological Lifshitz transition with anomalous exponent −2/5 that in the regime of small magnetic fields results in a Landau level energy collapse. Our work opens a new horizon in the study of the interplay between new correlated or superconducting phases and exotic van Hove singularities that are realized in Moiré superlattice materials.
Acknowledgments.-We have benefited from discussions with K. Kolář and F. von Oppen. This work was supported by the French National Research Agency (project SIMCIRCUIT, ANR-18-CE47-0014-01). The Flatiron Institute is a division of the Simons Foundation.

I. THE CONTINUUM MODEL FOR STTG
We use an effective continuum model for the low-energy electronic bands of the symmetric twisted trilayer graphene. Following Refs. [1][2][3][4] the Hamiltonian reads: wherek = −i∇ r , T (r) = 3 j=1 T j e iqj ·r , q 1 = k θ (0, 1), q 2/3 = k θ (∓ √ 3/2, −1/2), k θ = 4π/3 L s with a the graphene lattice parameter and L s = a/(2 sin θ/2) a/θ the linear size of the Moiré unit cell, U is the gate potential applied on the top and bottom layers and σ i ±θ/2 = j=x,y R ij (±θ/2) σ j with R ij (φ) rotation of angle φ around the z axis. In the Hamiltonian (1) we have also introduced the interlayer tunneling matrix where j = 0, 1, 2, w AB = w and w AA = w r with r < 1 due to lattice relaxation effects [5,6]. In the absence of the displacement field U = 0 the Hamiltonian (1) is invariant under the unitary transformation: which leaves the middle layer invariant while exchanging bottom and top layers. The mirror symmetry M z is characterized by two eigenstates | e 2 = (| 1 + | 3 )/ √ 2 and | e 1 =| 2 with eigenvalue +1, and one eigenstate | o = (| 1 − | 3 )/ √ 2 with −1. In the basis of the eigenstates of M z the Hamiltonian H(r) becomes: which is the staring point in the main text. In this basis we easily find that for U = 0 the model (4) is the sum of the TBG continuum model with interlayer hopping multiplied by √ 2 and an isolated Dirac cone. A finite displacement field, U = 0, breaks the M z symmetry and couples the additional Dirac cone with the bands of TBG. Assuming a small twist angle θ, we can neglect the phase factor in the Pauli matrices σ ±θ/2 → σ and the Hamiltonian (1) in the Moiré reciprocal lattice takes the form: In the previous expression (5) α and β denote the sublattice indices (α, β = A, B), MBZ stands for the Moiré Brillouin zone (grey shaded region in Fig. 1), l = 1, 2, 3 denotes the bottom, middle and top layers, Q ± is the triangular lattice for layer l = 2 and l = 1, 3, respectively, depicted as blue and red dots in Fig. 1 Figure 1: the k-space lattice generated by q 1 , q 2 and q 3 . The honeycomb lattice is composed by two interpenetrating triangular sublattices Q + and Q − corresponding to the middle (red dots) and the bottom/top (blue dots) layers, respectively. The grey shaded region centered at Γ depicts the MBZ, K and K are the high-symmetry points at the corner of the MBZ, while the M point is the middle point between K and K . The vectors b 1 = q 1 − q 2 and b 2 = q 1 − q 3 are the primitive vectors of the triangular lattice. Given Q 0 = Zb 1 + Zb 2 with Z the set of relative integers, the two sublattices Q + and Q − are obtained as Q + = Q 0 + q 2 and Q + = Q 0 − q 3 .
the eigenstates and eigenvalues of H(k) (5). Notice that H (1) has only three dimensionless parameter α = w/ v 0 k θ , r and u = U/ v 0 k θ . Before concluding the section we briefly list the symmetries of the Hamiltonian (1) that we are going to employ in the following. A detailed description of the symmetries of the model (1) is presented in Ref. [4]. At vanishing displacement field the Hamiltonian (1) commutes with the mirror symmetry M z , the threefold rotation C 3z around the z axis, the twofold rotation C 2x around the in-plane x axis and the C 2z T symmetry. In addition there are three useful transformations P , C 2x P and M z C 2x P . The first two, P and C 2x P , act differently on the even sector {| 2 , (| 1 + | 3 )/ √ 2} (upper 2 × 2 block of the matrix in Eq. (4)) and on the odd one {(| 1 − | 3 )/ √ 2} (lower one-dimensional block of the matrix in Eq. (4)). In particular we find that P anticommutes with the projection of the Hamiltonian in the even sector, while C 2x P anticommutes with both the even and odd mirror sectors. Moreover, C 2x P commutes with the displacement field term coupling the even (| 1 + | 3 )/ √ 2 and the odd (| 1 − | 3 )/ √ 2 states. From the previous considerations, we readily realize that the anticommuting symmetry of the Hamiltonian at finite U is obtained by combining M z and C 2x P to form M z C 2x P ({M z C 2x P, H} = 0). The introduction of a displacement field breaks M z and C 2x symmetries, while C 3z , C 2z T and M z C 2x P remain proper symmetries of the Hamiltonian (1).

The evolution of the Dirac cones as a function of the displacement field
In this section we provide a simple picture for the generation of electric field induced higher-order VHS by the fusion of zero-energy Dirac cones. The evolution of the zero-energy Dirac cones in the MBZ as a function of the applied external electric displacement field is shown in Fig. 2. The initial situation at zero displacement D = 0 consists of a total of three zero-energy Dirac cones with the same vorticity +1: two located at K and one at K. The Dirac cone at K is protected by C 2z T (forbids a gap opening) and M z C 2x P (pins it at zero-energy by enforcing a local particle-hole symmetric spectrum). None of these symmetries is broken by the displacement field such there is always a Dirac cone at K regardless of D. The situation is different at K where the two zero-energy Dirac cones are shifted in momentum space as D departs from zero, one along the K K line (red dots) and the other one along the K Γ line (violet dots). Further increasing D, six new zero-energy Dirac cones appear at a critical D ≥ 0.293V/nm in the middle of the Moiré Brillouin zone. There can be decomposed into two sets of threefold C 3z -symmetric cones with opposite vorticities ±1 (+1 as orange dots, −1 as blue dots). The +1 Dirac cones further move towards Γ as D increases, and the −1 Dirac  cones towards K. Eventually, in the region close to the higher-order VHS where the velocity v and the curvature η are small (D ∼ 0.34 − 0.366 V/nm for r = 0.8 and θ = 1.59 • ), the landscape of Dirac cones close to K is illustrated in Fig. 2(b), with the symmetry-protected +1 cone at K surrounded by six cones: three (red dots) of vorticity +1 originating from K and three (blue dots) of opposite vorticity. The higher-order VHS thus corresponds exactly to the merging of all six side Dirac cones with the central one at the K point. Above the critical value, D > D c , only the isolated Dirac cone located K with vorticity +1. Finally, we report in Fig. 3 the evolution of the density of states (DOS) as a function of the energy and of the displacement field for r = 0.8 and θ = 1.59 • . In correspondence of the point where the Fermi velocity vanishes D/ 0 0.366V/nm we see a strong enhancement of the DOS close to charge neutrality.
II. THE k · p HAMILTONIAN AT K The k · p Hamiltonian at K is obtained by performing degenerate perturbation theory on the zero energy doublet {| Φ ω , | Φ ω * } in the small perturbation δH(q) = H(K + q) − H(K) where H(k) is given in Eq. (5), | Φ ω and | Φ ω * are eigenstates of H(K) and C 3z with eigenvalues ω = e 2πi/3 and ω * . In the following we derive the k · p Hamiltonian to third order in the small deviation q = k − K. To this aim we are going to show for each order in q the terms allowed by the symmetries listed in the last part of Section I. Among them particularly important are C 3z = e 2πτ z i/3 , C 2z T = τ x K with K complex conjugation and M z C 2x P = τ x that applied to the k · p Hamiltonian at K give: At first order in perturbation theory we find that the only terms invariant under C 3z are q − τ + and q + τ − where q ± = q x ± iq y , τ ± = (τ x ± iτ y )/2, and τ 0,x,y,z are the Pauli matrices acting on the two-dimensional degenerate subspace {| Φ ω , | Φ ω * }. Moreover, any linear combination of the form r q + τ − + r * q − τ + with r ∈ C is invariant under the symmetry C 2z T . We employ the gauge freedom to choose the phases of the wave function | Φ ω and | Φ ω * so that r = v ∈ R. The latter constraint plays the role of a gauge-fixing condition and gives the first order correction The second order terms that are allowed by the C 3z symmetry are q 2 − τ − , q 2 + τ + , and q − q + τ 0 . Let us start by looking at the first two contributions that give rise to the quadratic correction f * q 2 x − q 2 y ) τ y + 2q x q y τ x . From the particle-hole symmetry M 2z C 2x P (6) it follows that Ref = 0 and f = −iη. Additionally, as a consequence of M 2z C 2x P the coefficient of q − q + τ 0 vanishes. Finally, we find that the second order correction reads: As a consequence of the C 3z symmetry at third order in the wave vector q we have q 3 − τ 0 , q 3 + τ 0 , q 2 − q + τ + and q 2 + q − τ − . The latter two terms give rise to the cubic contribution g q 2 + q − τ − +g * q 2 − q + τ + = Reg q 2 (q · τ )+Img q 2 (q x τ y − q y τ x ), that is invariant under C 2z T . Following the same line of reasoning used for the quadratic term, only q 2 (q · τ ) with q 2 = q 2 x + q 2 y is allowed by the particle-hole symmetry M z C 2x P . We conclude that Img = 0 and g = ξ ∈ R. Moreover, the same M z C 2x P symmetry selects the combination γ(q 3 − + q 3 + ) τ 0 . Therefore, the third order correction reads: By collecting the contributions the final Hamiltonian is, as given in the main text, where the coefficients v, η, γ and ξ are obtained by computing the expectation values: where Σ − Qα,Q β = σ − αβ δ Q,Q and we have introduced the projector in the high-energy states with L = {| Φ ω , | Φ ω * }, | u n =| u nK and n = nK . Before going on we remark that the specific form of the k · p model at K (10) is determined by the C 3z and the M z C 2x P symmetries. The breaking of one of these two symmetries would introduce additional q dependent corrections to H K (q) that spoil the higher-order VHS at θ c and D c . Indeed, in the presence of additional quadratic corrections it is not possible to induce an higher-order VHS with the only two tuning parameters θ and D. The evolution of the coefficients {v, ξ, η, γ} for the critical value of the twist-angle θ c 1.631 • as a function of the displacement field D at r = 0.8 is shown in Fig. 4. We find that when D/ 0 = D c / 0 0.5 V/nm the K point is an higher-order VHS where the first and second order coefficients both vanish v = η = 0 as highlighted by the inset of Fig. 4. Fig. 5 shows the trajectory of points where the higher-order VHS takes place together with its projections on the (r, θ) and (r, D) planes. As D increase, the trajectory extends till the end point at r 0.842 above which it is not possible to have v and η both vanishing. In the chiral limit r = 0 the trajectory originates from the first magic angle θ CH 1.536 • at vanishing displacement field. By reducing the twist-angle we find an additional trajectory of higher-order VHS, see Fig. 6, which starts from the second magic angle θ  By assuming that ξ, γ > 0 near the higher-order VHS the density of states reads: where φ + = arccos(ξ/2γ) if ξ/2γ ≤ 1 else φ + = 0, while φ − = arccos(−ξ/2γ) if ξ/2γ ≤ 1 else φ − = π. We find the DOS is power law divergent with ν = −1/3 and particle-hole symmetric. Interestingly, at the higher-order VHS there are two parameters γ and ξ that determines the properties of the iso-energy lines around K. The qualitative evolution of the energy contours as a function of the only dimensionless parameter 2γ/ξ at the higher-order VHS (v = η = 0) is shown in Fig. 7. We observe that in the regime of 2γ/ξ 1, see Fig. 7a calculated at 2γ/ξ = 0.1, the iso-energy lines are closed orbits almost rotationally symmetric winding around K. By taking the limit γ → 0 in Eq. (13) the DOS reads ρ(ω) = |ω| −1/3 /(6πξ 2/3 ). As we increase the ratio 2γ/ξ the iso-energy lines deform and the U (1) rotational symmetry is lowered to the discrete C 3z symmetry. Interestingly, by further increasing the ratio 2γ/ξ the effective model (10) shows a Lifshitz transition at 2γ/ξ = 1 where the iso-energy lines are open hyperbolic orbits with asymptotes forming an angle π/3 + (j − 1)2π/3, j = 1, 2, 3, with the x axis. Correspondingly, as shown in the main text the prefactor F (2γ/ξ) defined in Eq. (13) diverges and the power-law singularity becomes stronger than −1/3. In the next section we will show that higher-order corrections fix to ν = −2/5 the power-law singularity at the Lifshitz transition. Finally, for 2γ/ξ 1, see Fig. 7c, the energy contours are hyperbolic orbits characterized by asymptotes π/6 + (j − 1)π/3, j = 1, 2, 3, 4, 5, 6, with the x axis. These iso-energy contours have also been predicted in ABC stacking trilayer graphene on boron nitride in Ref. [7] and in biased bilayer graphene in Ref. [8]. By taking the limit ξ/2γ → 0 in Eq. (13) we obtain ρ(ω) = Γ(7/6)|ω| −1/3 /[π 3/2 (2γ) 2/3 Γ(2/3)] where Γ(z) is the Euler Gamma function. The evolution of ξ and 2γ along the trajectory of higher-order VHS originating from the first magic angle at r = 0 is shown in the main text. Starting from r 0.81 we find at first that the value of ξ is much larger that γ, 2γ/ξ 1. Then, by reducing the atomic corrugation r the ratio 2γ/ξ displayed in the main text grows monotonically till the value r 0.34 (θ c 1.543 • , D/ 0 0.086 V/nm) where 2γ/ξ crosses the critical point 1. Here, the Lifshitz transition takes place and the power-law singularity of the density of states changes from ν = −1/3 to −2/5. Finally, for r < 0.34, 2γ/ξ becomes larger than 1 and increases approaching the chiral limit r = 0. Thus, by changing the value of r the constant energy lines at the higher-order VHS evolves from 2γ/ξ 1 to 2γ/ξ 1.

B. The density of states at the Lifshitz transition
At the Lifshitz transition 2γ = ξ and to leading order in the k · p expansion the dispersion relation reads: where q = q 2 x + q 2 y ≥ 0. The vanishing of the dispersion relation E + (q) and E − (q) along the lines φ + j = π/3 + 2π(j − 1)/3 and φ − j = 2π(j − 1)/3 with j = 1, 2, 3 gives rise to the divergence of F(2γ/ξ) in Eq. (13) as displayed in the main text. In order to characterize the singularity we include next to leading order corrections in the wave vector q to the k · p Hamiltonian (10). By employing the symmetries of the model (6) at K we find that the q 4 correction reads: while the q 5 one is: By taking into account H K (q) and H K (q), and expanding to the fifth order the dispersion relation, we find: E ± (q) = ξq 3 (±1 + cos 3φ) ∓ (µ + ν) q 4 sin 3φ + 2ζ cos 3φ ± (µ − ν) 2 + 4wξ + (4ξΥ + (µ − ν) 2 ) cos 6φ 4ξ q 5 . (17) From the latter expression, we realize that the fourth order term, proportional to ∝ q 4 sin 3φ, vanishes along the directions φ + j and φ − j where the singularity in the angular integral takes place and the first non-vanishing correction is q 5 . The bottom panel of Fig. 8 shows the lowest energy electronic bands ±k along the line k = K + κ with κ k x obtained by diagonalizing the model (1) at the higher-order VHS (2γ = ξ). In agreement with the results of the k · p expansion the top panel of Fig. 8 shows that for κ = −q < 0 the dispersion +K−qx goes as ∼ q 5 while for positive wave vectors deviations κ = q > 0 we find +K+qx ∼ q 3 .
To determine the power-law singularity of the density of states we expand the dispersion relation E ± (q) for angles φ close to the singular lines φ ± j keeping terms up to the fifth order: where as shown in Fig. 8 the parameter Λ is positive Λ > 0. By considering ω > 0 the wave vector q 0 that satisfies the delta function constraint in Eq. (13) is now solution of ω = ξ q 3 δφ 2 /2 + (µ + ν) δφ q 4 + Λ q 5 .
where we have projected the terms (30) and (32) in the basis of eigenstates of the harmonic oscillator | n with n = 0, 1, · · · . The fourth order contribution is reduced by a factor 1/g with respect to the third order one. For a magnetic field that ranges from B = 0.005 − 1.0 T and in the magic-angle region k θ 0.47 nm −1 (L s 8.8 nm) the value of g goes from ∼ 170 to ∼ 12. In order to find the spectrum of the Hamiltonian in Eq. (34) we truncate the basis of harmonic oscillator states to a given integer N s , i.e. n = 0, 1, · · · , N s − 1. The number of states N s needed to sample the turning points produced by the q 4 correction goes as N s ∼ (k θ l B ) 2 = g 2 . Therefore, we conclude that N s increases like 1/B as we decreases the value of B. For a magnetic field of B = 0.5 T the value of g in the magic angle region is g 16.6 and we expect to achieve convergence for N s > 275. This is confirmed by Fig. 10 where we show the Landau level energy as a function of N s for 2γ/ξ = 1.73. Finally, Fig. 11 shows the Landau level spectrum across the Lifshitz transition 2γ/ξ = 1 for a magnetic field B = 0.5 T. states. The calculation has been performed by keeping N s = 393 harmonic oscillator states for | Φ ω and | Φ ω * .