Magic of high-order van Hove singularity

The van Hove singularity in density of states generally exists in periodic systems due to the presence of saddle points of energy dispersion in momentum space. We introduce a new type of van Hove singularity in two dimensions, resulting from high-order saddle points and exhibiting power-law divergent density of states. We show that high-order van Hove singularity can be generally achieved by tuning the band structure with a single parameter in moiré superlattices, such as twisted bilayer graphene by tuning twist angle or applying pressure, and trilayer graphene by applying vertical electric field. Correlation effects from high-order van Hove singularity near Fermi level are also discussed.

T he recent discovery of unconventional insulating states and superconductivity in twisted bilayer graphene (TBG) 1,2 has attracted enormous interest. At ambient pressure, these intriguing phenomena of correlated electrons only occur near a specific twist angle θ % 1:1 , widely referred to as the magic angle. The existence of such a magic angle was first predicted from band structure calculations. Based on a continuum model 3,4 , Bistritzer and MacDonald showed 5 that the lowest moiré bands in TBG become exceptionally flat at this magic angle, and are hence expected to exhibit strong electron correlation. This pioneering work inspired a large body of experimental works in recent years, which culminated in the discovery reported in refs. 1,2 . The origin of magic-angle phenomena, i.e., the emergence of superconductivity and correlated insulator, is now under intensive study [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] .
While the reduction of bandwidth is undoubtedly important, it is evident that the phenomenology of magic-angle TBG cannot be ascribed to a completely flat band. Quantum oscillation at low magnetic fields reveals doping-dependent Fermi surfaces within the narrow moiré band. Detailed scanning tunneling spectroscopy (STS) studies [26][27][28] and compressibility measurements 29 show that the moiré bandwidth at a magic angle is a few tens of meV, larger than what previous calculations found 5,[30][31][32][33][34][35] . Moreover, there is no direct evidence from existing STS measurements at various twist angles that the moiré band is exceptionally flat right at the magic angle (see also ref. 28 ). These experimental results motivated us to consider additional feature of electronic structure which may create favorable condition for correlated electron phenomena in magic-angle TBG.
One such feature is the van Hove singularity (VHS) in moiré bands. The importance of VHS has already been recognized in our theory of correlated TBG 36,37 and related studies [38][39][40][41][42][43][44] using the weak coupling approach. Generally speaking, VHS with divergent density of states (DOS) in two-dimensional systems are associated with saddle points of energy dispersion in k space. When a VHS is close to Fermi energy, the increased DOS amplifies electron correlation, resulting in various ordering instabilities, such as density wave and superconductivity at low temperature 36 . Indeed, the two recent STS measurements on gatetunable TBG around magic angle 26,27 find that when the VHS shifts to the Fermi level under gating (the corresponding density is within 10% of half filling), the VHS peak in tunneling density of states splits into two new peaks (see also ref. 45 ). These findings clearly demonstrate the prominent role of VHS in magic-angle TBG. However, it is also known from general consideration and previous STS measurements [46][47][48] that VHS are present at all twist angles. It is therefore unclear whether VHS at magic angle is anything special.
In this work, we propose a new perspective that relates magic angle to a new type of VHS in the single-particle energy spectrum of TBG. We show that as the twist angle decreases below a critical value θ c , the van Hove saddle point-which marks the change of topology in Fermi surface (Lifshitz transition)-undergoes a topological transition whereby a single saddle point splits into two new ones. Right at θ c , the saddle point changes from secondorder to higher-order; the DOS at VHS is significantly enhanced from logarithmic to power-law divergence, which promises stronger electron correlation. We propose that proximity to such "high-order van Hove singularity", which requires tuning to the critical twist angle or pressure, is an important factor responsible for correlated electron phenomena in TBG near half filling.
We demonstrate by topological argument that high-order VHS can be generally achieved by tuning the band structure with just a single control parameter. For TBG, besides twist angle, pressure can also induce superconductivity at a twist angle larger than 1:1 49 . By tracking the evolution of the moiré band with twist angle and pressure, we locate the high-order VHS in experimentally relevant parameter regime and predict its key feature, a distinctive asymmetric peak in DOS. This feature compares well with the experimentally observed VHS peak in magic-angle TBG, as shown in Fig. 1. We also discuss the splitting of VHS peak by nematic or density wave order.

Results
Ordinary and high-order VHS. In two-dimensional (2D) electron systems with energy dispersion EðkÞ, an ordinary VHS with logarithmically diverging DOS occurs at a saddle point k s , determined by  where D ij 1 2 ∂ i ∂ j E is the 2 2 Hessian matrix of E at k s . Since D is symmetric by definition, we can rotate the axes to diagonalize D and Taylor expand the energy dispersion near k s as E À E v ¼ Àαp 2 x þ βp 2 y , where E v is the VHS energy, the momentum p ¼ k À k s is measured from the saddle point, and the coefficients Àα, β are the two eigenvalues of D with Àαβ ¼ det D < 0. This dispersion describes two pieces of Fermi contours that intersect at the saddle point k s . It is known from a topological consideration that due to the periodicity of EðkÞ on the Brillouin zone (a torus), van Hove saddle points appear quite generally in energy bands of two-dimensional materials 50 .
We now introduce the concept of high-order VHS associated with a high-order saddle point defined by the following condition: which imply αβ ¼ 0. There exist two types of high-order VHS.
The first type corresponds to α ¼ β ¼ 0, i.e., D ¼ 0. The Taylor expansion of EðkÞ around such saddle point then starts from at least the third order, and describes an intersection of three or more Fermi surfaces at a common k point 51 . Such high-order VHS (type-I) is also referred to as multicritical VHS in the recent literature [52][53][54] .
Here we focus instead on a new type of high-order VHS relevant for TBG, where the Hessian matrix D has a single zero eigenvalue (β ¼ 0), while the other eigenvalue is nonzero (α ≠ 0). Then in the Taylor expansion of EðkÞ, besides the single secondorder term α we must also include higher order terms in order to capture the Fermi contours nearby. Importantly, unlike the previous case, type-II high-order VHS still describes the touching of two Fermi surfaces, but they touch tangentially (generally speaking) rather than intersect at a finite angle as in the case of ordinary VHS.

VHS in TBG.
We now turn to VHS in TBG, whose low-energy moiré bands arise from inter-layer coupling of Dirac states on the two graphene layers. Since single-particle scattering between K and K 0 points requires large momentum transfer, it is suppressed at small twist angles due to the long wavelength of moiré potential. In the absence of valley hybridization, the moiré bands from K and K 0 valleys are decoupled [3][4][5]55 , with energy dispersions denoted by E þ ðkÞ and E À ðkÞ, respectively. Time-reversal symmetry implies E À ðkÞ ¼ E þ ðÀkÞ, so it suffices to consider E þ ðkÞ only in the following.
The number and location of VHS points in TBG depend on the twist angle. For example, the band structure calculation for θ ¼ 2 56 reveals the existence of three symmetry-related van Hove saddle points on ΓM lines in the mini-Brillouin zone (MBZ), shown in Fig. 2a. Across the VHS energy, the Fermi contour changes from two disjoint Dirac pockets around MBZ corners to a single pocket enclosing the MBZ center. This leads to the conversion between electron and hole charge carriers, as evidenced by the sign change of Hall coefficient 56 . On the other hand, at θ ¼ 1:05 , on each ΓM line there is a local energy maximum-instead of saddle point-in the moiré valence band. As doping increases, an additional hole pocket (not present at θ ¼ 2 ) emerges out of each energy maximum and eventually intersects two Dirac pockets at two new saddle points on opposite sides of ΓM 32 . In such band structure there are a total of 6 VHS points, shown in Fig. 2c.
By continuity, we deduce from this change in the number of VHS points that a topological transition of saddle points must occur, at a hitherto unknown critical twist angle, in such a way that a saddle point on ΓM (denoted by Λ 0 ) splits into a pair of new ones off ΓM (denoted by Λ ± ). Importantly, the behavior of these VHS points in the vicinity of this transition is solely governed by the local energy dispersion near Λ 0 .
We now expand EðkÞ near Λ 0 to higher orders: where p x (p y ) is parallel (perpendicular) to ΓM line. In the expansion (3), first-order terms vanish because Λ 0 has by definition zero Fermi velocity. Moreover, only even power terms of p y are allowed because of the two-fold rotation symmetry of TBG, which acts within each valley and maps ðp x ; p y Þ to ðp x ; Àp y Þ. The third-order γ term and fourth-order κ term are essential to describe the splitting of VHS across the transition, which we shall show below along with the relation to scaling properties. The behavior of VHS and Fermi contour of the energy dispersion (3) depends crucially on the sign of β. For β > 0, p ¼ 0 (i.e., Λ 0 ) is an ordinary van Hove saddle point with logarithmic divergent DOS. The Fermi contour at the VHS energy consists of two curves that approach straight lines p y =p Þ as shown in Fig. 2d. This behavior corresponds to the VHS in the calculated band structure of TBG at θ ¼ 2 . When β ! 0 þ , φ ! 0 so that the two Fermi contours at the VHS energy tend to touch tangentially. On the other hand, for β < 0, p ¼ 0 becomes a local energy maximum, while two new saddle points appear at momenta . Throughout this manuscript we consider the regime γ 2 þ 4ακ > 0 so thatγ is always real and positive, and two split saddle points Λ ± are well-defined. Λ þ and Λ À are a pair of ordinary VHS points related by two-fold rotation, whose Fermi contours are two parabolas 2αðp x Ç β=γÞ ¼ ðγ ±γÞp 2 y , see Fig. 2f. This behavior corresponds to the VHS in the calculated band structure at θ ¼ 1:05 . In the limit β ! 0 À , Λ ± approaches Λ 0 .
Our unified description of two regimes of ordinary VHS (Λ 0 and Λ ± ) using a single local energy dispersion (3) implies that the transition between them corresponds to the sign change of β. In this process, energy contours at VHS changes from intersecting at one point Λ 0 (Fig. 2d and e) to two points Λ ± (Fig. 2f). Right at β ¼ 0, p ¼ 0 becomes a high-order VHS point. The Fermi contour in its vicinity consists of two parabolas 2αp x ¼ ðγ ±γÞp 2 y , touching each other tangentially at p ¼ 0. In the special case with κ ¼ 0, one of the parabolas become a straight line, as shown in Fig. 2e.
The DOS can be used as an indicator of this topological transition of VHS. For the energy dispersion (3), we find the DOS analytically where z ± ¼ ðβ=αÞ ± ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðβ=αÞ 2 þ ε q and ε ¼γ 2 ðE À E v Þ=α 3 is energy deviation from VHS in dimensionless unit. The integration over p is extended to infinity. Here KðzÞ is the complete elliptic integral of the first kind, sgnðrÞ is the sign function defined as sgnðrÞ ¼ À1 for r < 0 and 1 for r ! 0, and ΘðrÞ 1 2 ½1 À sgnðÀrÞ is the step function so that Θð0Þ ¼ 0. The DOS with different β in Eq. (5) are shown in Fig. 3a. In the following we will discuss the asymptotic behavior of the DOS near ordinary and high-order VHS. When α > 0 and β ≠ 0, the DOS diverges logarithmically at ordinary VHS: where Λ ¼ 64αβ 2 =γ 2 is the high-energy cutoff. In the limit β ! 0, the prefactor in Eq. (6) diverges as 1= ffiffiffiffiffi ffi jβj p , indicating a strong increase of DOS as the VHS approaches a high-order one. Right at the transition point β ¼ 0, we find power-law divergent DOS near high-order VHS where C ¼ ð2πÞ À 5 2 Γð 1 4 Þ 2 ¼ 0:133. This power-law divergence with exponent À1=4 is stronger than the logarithmic divergence at ordinary VHS. Another key difference is that the diverging part of DOS around this high-order VHS is inherently asymmetric with respect to E v , as reflected by the ffiffi ffi 2 p factor in Eq. (7). This is in contrast with ordinary VHS where the DOS above and below are asymptotically symmetric. Using Eq. (5), we can fit the experimental data of tunneling conductance in ref. 26 with finite broadening, as shown in Fig. 1a and Supplementary Material Sec. I.
At the high-order VHS (β ¼ 0), the dispersion has the following scaling invariance with arbitrary λ > 0. Here we set E v ¼ 0 for simplicity. Including higher order terms or the third and fourth order terms other than γ; κ in Taylor expansion (3) will break scaling invariance (8).
Notice that the scaling dimensions along p x and p y directions are different ½p x ¼ 1 2 ; ½p y ¼ 1 4 , with scaling dimension of energy ½E ¼ 1. From the scaling invariance (8) we immediately obtain the scaling dimension of DOS as ½ρ ¼ ½p x þ ½p y À ½E ¼ À 1 4 , the same as Eq. (7). To directly reveal the power-law behavior of experimental data, in Fig. 1b we plot energy deviation jE À E v j and experimental DOS including background contribution both in logarithmic scale. In the log-log plot, two sides of the highorder VHS peak form two parallel lines with the same slope  (3) in d-f. Thick lines are energy contours at VHS energy. In band structure calculation, g 0 ¼ 1:2 g and the parameters are a g ¼ 1, b g ¼ 1:9, and c g ¼ 2.
Besides α and β, the Fermi contours in momentum space depend on both γ and κ, while the DOS in energy domain depends only onγ. This is because the same scaling dimension of p x and p 2 y allows the nonlinear transform under which the dispersion (3) becomes E À E v ¼ Àαp 2 x þ βp 2 y þγ 2p4 y =ð4αÞ. The nonlinear transform (9) changes the shape of Fermi contours while it preserves area element and hence leaves DOS in Eq. (5) invariant.
We have thus established by continuity argument the existence of high-order VHS in the moiré band of TBG at a certain critical twist angle θ c . The value of θ c depends on the model we employ. As an illustrative example, we consider the continuum model [3][4][5]32 of TBG, where approximate particlehole symmetry holds, and the Fermi surface topology is solely determined by two dimensionless parameters g ¼ λu=v; g 0 ¼ λu 0 =v. Here u and u 0 denote interlayer hoppings in AA and AB regions respectively, v is the monolayer Dirac velocity and λ is the moiré wavelength. We further fix the ratio g 0 =g ¼ u 0 =u ¼ 1:2 as calculated for corrugated structure at θ $ 1 32 . At given g 2 ½1; 2, we calculate the Taylor coefficients α; β; γ; κ by numerical derivatives of energy with respect to momentum at VHS on the hole side, and the result is shown in Fig. 3b. It can be found that γ 2 þ 4ακ is always positive when g 2 ½1; 2, and across g c % 1:995, β changes sign (Fig. 3b inset) while α; γ; κ do not. With appropriate and reasonable values of u and v, from g c we can obtain the critical twist angle θ c % 1 , consistent with experiments. Details of continuum model can be found in the Supplementary Material Sec. II. We may include more ingredients such as lattice relaxation 31,57,58 and strain 59 to better describe TBG, thus the Fermi surface topology will depend on more details of the model. Nevertheless, the existence of high-order VHS and critical twist angle θ c is robust 59 .
Many-body phenomena near high-order VHS. When chemical potential is near high-order VHS, many-body effects can be drastic as the strongly diverging DOS leads to divergences in noninteracting susceptibilities of various channels. This signals a strong tendency to various broken symmetry states around van Hove filling, when electron-electron interaction is taken into account. Indeed, the recent STS measurements found that when the Fermi energy approaches the van Hove energy under doping, the VHS peak in DOS splits into two new ones 26,27 .
A detailed analysis of interacting electrons near high-order VHS is beyond the scope of this work. Instead we discuss two possible scenarios for the observed splitting of the VHS peak near the Fermi energy.
First, strains in experimental samples can split DOS peak by breaking rotation symmetries which relate VHS along different directions. Though the energy splitting due to strain may be small at single-particle level, Coulomb interaction U can give rise to the Stoner-type enhancement factor ð1 À Uχ s Þ À1 within the random phase approximation (RPA). Since the noninteracting susceptibility χ s reflects the divergent DOS, strain effect on high-order VHS can be greatly enhanced by interaction.
Second, an intervalley density wave order can split a DOS peak by spontaneously breaking translation symmetry. To show such a density wave instability, we calculate intervalley susceptibility χðqÞ with finite momentum q. There are in total six high-order VHS from two valleys, located at three ΓM lines (Fig. 4a). When we consider two high-order VHS along the same ΓM line but from opposite valleys, χðqÞ is sharply enhanced near wave vector q ¼ Q due to the power-law divergent DOS from both valleys, where Q ¼ 2k s is the momentum separation between the two VHS (Fig. 4a). As a result, when approaching low temperatures, the intervalley susceptibility χðQÞ of high-order VHS diverges more rapidly than logarithm (Fig. 4b), indicating a stronger tendency toward density wave instability than one with ordinary VHS. The DOS peak at VHS energy E v is split into two in the presence of intervalley density wave (Fig. 4c).

Discussion
To summarize, we propose that the proximity to high-order VHS underlies correlated electron phenomena in TBG near magic angle. We reveal a distinctive feature of high-order VHS -power-law divergent DOS with an asymmetric peak, which compares well with the recent STS data on magic-angle TBG. We also discuss nematic and density wave instabilities due to electron-electron interaction near van Hove filling, and illustrate the splitting of VHS in these broken symmetry states.
It is worth examining the similarity and difference between the two scenarios favoring correlated electron phenomena in TBG: flat band 5 and (high-order) VHS. Both scenarios rely on large DOS enhancing electron correlation. Obviously, the largest possible DOS is realized in the limit of a completely flat band 60 . However, a completely flat band is difficult to achieve under realistic conditions. On the other hand, the VHS scenario is more likely to be relevant when the bandwidth is not so small in comparison to electron interaction. In this scenario, the DOS near Fermi energy matters more than those far away, and it is further increased by proximity to high-order VHS. Importantly, to achieve high-order VHS generally requires tuning the band structure with only a single parameter, whereas to achieve a completely or extremely flat band usually requires more fine tuning.
Broadly speaking, moiré superlattices in 2D materials offer an unprecedented platform for studying correlated electron phenomena near VHS. Electrostatic gating enables doping these systems to van Hove filling without introducing disorder. The great tunability of moiré band structure by twist angle, electric field, strain and other means grants access to high-order VHS, where the strongly divergent DOS promises a plethora of manybody phenomena. Such examples include TBG discussed in the maintext and trilayer graphene on boron nitride discussed in Supplementary Material Sec. III. The time is coming for creating strong electron correlation by designing single-particle dispersion around VHS.

Methods
Fitting of tunneling conductance peaks. In Eq. (5) we assume the dispersion (3) extends to the whole momentum space, which holds only for a finite momentum range in realistic systems. Assuming the dispersion over the whole momentum space as E p , whose Taylor expansion is E p in the finite momentum range jpj < Λ, we find the total DOS of such dispersion can be written as where ρðEÞ is Eq. (5) if the Taylor expansion is Eq. (5), and ρ c is the background contribution When DOS of the tip is featureless, the tunneling conductance from the tip to the sample will be proportional to the total DOS ρ tot of the sample. There are two distinct types of divergent DOS and hence conductance peaks in this work and also in experimental data. The first type is due to ordinary VHS, which is symmetric and logarithmic, and another type is due to high-order VHS, which is intrinsically asymmetric and power-law. These two types of conductance peaks follow different functional forms as follows (o denotes ordinary VHS and h denotes high-order VHS) each with three parameters A; E v and G c . Among them, E v denotes the energy of conductance peak, and G c is the background contribution. In Fig. 1b, we plot conductance data G near the peaks at twist angle 1:10 , where energy and conductance are both plot in logarithmic scale. Since the VHS is high-order, conductance peak is asymmetric and power-law, two sides of the peak will follow two parallel lines with the same slope −1/4 but different vertical intercepts.
We make the substitution E ! E þ iη in the expressions above of tunneling conductance to describe broadening effect. The two sides of the DOS peak of ordinary VHS are broadened in the same way, while the two sides of the DOS peak of high-order VHS are broadened differently due to the intrinsic particle-hole asymmetry. As a result, after broadening, E v coincides with the DOS peak energy E m when VHS is ordinary, while for high-order VHS, E v can deviate from E m as shown in Fig. 1a.
Mean-field Hamiltonian of interacting high-order VHS. With finite intervalley density wave order parameter Δ, the low-energy Hamiltonian of two high-order VHS on the same ΓM line can be written as HðpÞ ¼ E v À αp 2 x þ κp 4 y þ γp x p 2 y τ z þ Δτ x ; where Pauli matrices τ act on the valley indices. From this Hamiltonian Fig. 4c is plotted, with order parameter Δ ¼ 0:01E 0 and broadening η ¼ 5 10 À3 E 0 .

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.