Magic of high order van Hove singularity

TL;DR

This paper introduces high-order van Hove singularities in two-dimensional materials, where saddle points lead to power-law divergences in the density of states, potentially enhancing electron correlation effects.

Contribution

It presents the concept of high-order van Hove singularities and demonstrates how they can be realized in moiré superlattices through tuning parameters like twist angle or pressure.

Findings

High-order van Hove singularities cause power-law divergence in density of states.

Such singularities can be achieved in moiré materials by tuning a single parameter.

Enhanced electron correlation effects are expected near these singularities.

Abstract

We introduce a new type of van Hove singularity in two dimensions, where a saddle point in momentum space is changed from second-order to high-order. Correspondingly, the density of states near such ``high-order van Hove singularity'' is significantly enhanced from logarithmic to power-law divergence, which promises stronger electron correlation effects. High-order van Hove singularity can be generally achieved by tuning the band structure with a single parameter in moir\'e superlattices, such as twisted bilayer graphene by tuning twist angle or applying pressure, and trilayer graphene by applying vertical electric field.

Figures

a of maintext), while at $g=2$ there are six ordinary VHS on two sides of $\Gamma M$ (Fig. 2c of maintext). According to our analysis above, this implies a topological transition of VHS occurs at $g=g_{c}\in(1,2)$. To determine the value $g_{c}$, we calculate the Taylor coefficients $\alpha,\beta,\gamma,\kappa$ by numerical derivatives of energy dispersion with respect to momentum components $k_{x},k_{y}$ at VHS point. As shown in Fig. 3b of maintext, $\alpha,\beta,\gamma,|\kappa|$ all decrease as $g$ increases, and across $g_{c}\approx 1.995$, $\beta$ changes sign (Fig. 3b inset) while $\alph

c of maintext is plotted, and the order parameter is $\Delta=0.01E_{0}$ with broadening $\eta=5\times 10^{-3}E_{0}$.

a of maintext), while at $g=2$ there are six ordinary VHS on two sides of $\Gamma M$ (Fig. 2c of maintext). According to our analysis above, this implies a topological transition of VHS occurs at $g=g_{c}\in(1,2)$. To determine the value $g_{c}$, we calculate the Taylor coefficients $\alpha,\beta,\gamma,\kappa$ by numerical derivatives of energy dispersion with respect to momentum components $k_{x},k_{y}$ at VHS point. As shown in Fig. 3b of maintext, $\alpha,\beta,\gamma,|\kappa|$ all decrease as $g$ increases, and across $g_{c}\approx 1.995$, $\beta$ changes sign (Fig. 3b inset) while $\alph

a of maintext), while at $g=2$ there are six ordinary VHS on two sides of $\Gamma M$ (Fig. 2c of maintext). According to our analysis above, this implies a topological transition of VHS occurs at $g=g_{c}\in(1,2)$. To determine the value $g_{c}$, we calculate the Taylor coefficients $\alpha,\beta,\gamma,\kappa$ by numerical derivatives of energy dispersion with respect to momentum components $k_{x},k_{y}$ at VHS point. As shown in Fig. 3b of maintext, $\alpha,\beta,\gamma,|\kappa|$ all decrease as $g$ increases, and across $g_{c}\approx 1.995$, $\beta$ changes sign (Fig. 3b inset) while $\alph

Tables

Table 1. Table 1: Fitting parameters of tunneling conductance data at twist angle θ 𝜃 \theta . In the last column ± plus-or-minus \pm denote positive/negative energy band respectively. For θ = 2.02 ∘ 𝜃 superscript 2.02 \theta=2.02^{\circ} , G m subscript 𝐺 m G_{\rm m} denotes the maximal value of tunneling conductance data.

$\theta$	$G_0$ (nS)	$G_c$ (nS)	$E_v$ (meV)	$\eta$ (meV)	$w$ (meV)	band
2.02∘	6.112$G_m$	3.108$G_m$	-94.15	23.29	49.69	$-$
1.10∘	98.22	67.57	-36.05	2.786	3.360	$-$
	76.34	57.6	16.72	3.556	-0.1442	$+$
0.79∘	30.89	16.22	-3.299	1.822	-4.271	$-$
	32.82	16.60	0.7944	2.074	-5.8787	$+$

Equations

$${\nabla }_{\mathbf{k}}E({\mathbf{k}}_s)=0\text{ and }\det D<0,$$

$${\nabla }_{\mathbf{k}}E({\mathbf{k}}_s)=0\text{ and }\det D=0,$$

$$E-E_{\mathrm{v}}=-\alpha p_x^2+\beta p_y^2+\gamma p_x{p}_y^2+\kappa p_y^4+\dots ,$$

$${\Lambda }_{\pm }=(-\beta \gamma /{\tilde{\gamma }}^2,\pm \sqrt{-2\alpha \beta }/\tilde{\gamma })$$

$$\rho (E)=\frac{-1}{\pi }\mathrm{Im}\int \frac{d^2\mathbf{p}}{(2\pi {)}^2}\frac{1}{E+i{0}^{+}-E(\mathbf{p})}=\frac{\text{sgn}(\alpha \beta )}{\sqrt{2}\alpha {\pi }^2}\mathrm{Re}\left[\frac{1}{\sqrt{z_{-}}}\mathrm{K}\left(1-\frac{z_{+}}{z_{-}}\right)-\frac{2i}{\sqrt{z_{+}}}\mathrm{K}\left(\frac{z_{-}}{z_{+}}\right)\Theta (-\alpha \beta )\right],$$

$$\rho (E)=\frac{1}{4{\pi }^2\sqrt{\alpha |\beta |}}\times \{\begin{array}{ll}\log \frac{\Lambda }{|E-E_{\mathrm{v}}|}& \beta >0\\ \sqrt{2}\log \frac{\Lambda }{|E-E_{\mathrm{v}}-\delta |}& \beta <0\end{array}$$

$$\rho (E)$$

$$E({\lambda }^{\frac{1}{2}}{p}_x,{\lambda }^{\frac{1}{4}}{p}_y)=\lambda E(p_x,p_y)$$

$${\tilde{p}}_x=p_x-\frac{\gamma }{2\alpha }{p}_y^2,{\tilde{p}}_y=p_y,$$

$$E-E_{\mathrm{v}}=-\alpha p_x^2+\beta p_y^2+\gamma p_x{p}_y^2+\kappa p_y^4,$$

$$\rho (E)=\frac{1}{\sqrt{2}\alpha {\pi }^2}Q(\epsilon ,r),\epsilon =\frac{{\tilde{\gamma }}^2}{{\alpha }^3}(E+i\eta -E_{\mathrm{v}}),r=\frac{\beta }{\alpha },$$

$$f(\epsilon ,r)=\frac{1}{\sqrt{z_{-}}}\mathrm{K}\left(1-\frac{z_{+}}{z_{-}}\right),g(\epsilon ,r)=\frac{2}{\sqrt{z_{+}}}\mathrm{K}\left(\frac{z_{-}}{z_{+}}\right)$$

$$\mathrm{K}(z)={\int }_0^{\pi /2}\frac{d\theta }{\sqrt{1-z{\sin }^2\theta }}.$$

$$G(E)=G_{0}Q(\epsilon ,r)-G_c,$$

$$w(2.02^{\circ })\gg w(1.10^{\circ })\approx 0>w(0.79^{\circ }),$$

$$G_o(E)=A\log |E-E_{\mathrm{v}}|-G_c,G_h(E)=A\mathrm{Re}(E_{\mathrm{v}}-E{)}^{-1/4}-G_c,$$

$$H(\mathbf{k},\mathbf{r})=\left(\begin{array}{cc}v{R}_{\theta /2}{\mathbf{p}}_1\cdot \mathbf{\sigma }& U(\mathbf{r})\\ U^{†}(\mathbf{r})& v{R}_{-\theta /2}{\mathbf{p}}_2\cdot \mathbf{\sigma }\end{array}\right),U(\mathbf{r})=\sum _{j=0}^2\left(\begin{array}{cc}u& u^{\prime }{\omega }^{-j}\\ u^{\prime }{\omega }^j& u\end{array}\right)e^{i{\mathbf{G}}_j\cdot \mathbf{r}},$$

$$H({\lambda }^{-1}\mathbf{k},\lambda \mathbf{r})=\frac{v}{\lambda }\left(\begin{array}{cc}{\mathbf{p}}_1\cdot \mathbf{\sigma }& f(\mathbf{r})\\ f^{†}(\mathbf{r})& {\mathbf{p}}_2\cdot \mathbf{\sigma }\end{array}\right),f(\mathbf{r})=\sum _{j=0}^2\left(\begin{array}{cc}g& g^{\prime }{\omega }^{-j}\\ g^{\prime }{\omega }^j& g\end{array}\right)e^{i{\mathbf{n}}_j\cdot \mathbf{r}},{\mathbf{n}}_0=0,{\mathbf{n}}_{1,2}=(\mp \frac{1}{2},\frac{\sqrt{3}}{2}).$$

$$\mathcal{H}(\mathbf{p})=E_{\mathrm{v}}-\alpha p_x^2+\kappa p_y^4+\gamma p_x{p}_y^2{\tau }_z+\Delta {\tau }_x,$$

$$H(\mathbf{p})=-\mu +a{p}^2+(\gamma +b{p}^2){\sigma }_x+c\left(\begin{array}{cc}0& p_{-}^3\\ p_{+}^3& 0\end{array}\right)$$

$$\mu =0.0027,a=3.2708,\gamma =-0.0083,b=3.2165,c=\mathrm{91.1443.}$$

$$V_{\xi }(\mathbf{r})=V_{\xi }\frac{1+\xi {\sigma }_z}{2}\sum _{j=0}^2\cos ({\mathbf{G}}_j\cdot \mathbf{r}+{\varphi }_{\xi })+V{\sigma }_z$$

$$\delta =-0.017,V_{-}=12.09\mathrm{meV},{\varphi }_{-}=-0.2591\pi .$$

$$E(\mathbf{p}-\mathbf{K})=E_{\mathrm{v}}+\alpha p^2+\gamma (p_x^3-3p_x{p}_y^2)$$

$$H_{tb}=\sum _{⟨ij⟩}{t}_1{c}_i^{†}{e}^{i\varphi {\sigma }_z}{c}_j+\sum _{⟨⟨ij⟩⟩}{t}_2{c}_i^{†}{c}_j+h.c.$$

$$\delta =-0.017,V_{+}=-14.88\mathrm{meV},{\varphi }_{+}=0.2788\pi ,V=70\mathrm{meV},$$

$$t_1=1.20\mathrm{meV},\varphi =-0.48\pi ,t_2=0.45\mathrm{meV}.$$

Full text

The recent discovery of unconventional insulating states and superconductivity in twisted bilayer graphene (TBG) exp1 ; exp2 has attracted enormous interest. At ambient pressure, these intriguing phenomena of correlated electrons only occur near a specific twist angle $\theta\approx 1.1^{\circ}$, 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 Neto ; Neto2 , Bistritzer and MacDonald showed MacDonald 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 Ref. exp1 ; exp2 . The origin of magic-angle phenomena, i.e., the emergence of superconductivity and correlated insulator, is now under intensive study Xu ; Volovik ; Yuan ; Po1 ; Kivelson ; Phillips ; Lee ; Oskar ; Thomson ; Biao ; Venderbos ; Fidrysiak ; Ying ; Guinea ; Ochi ; XCWu ; FCWu ; FCWu1 ; Ray ; Zhang .

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 Abhay ; Stevan ; Yazdani ; LeRoy and compressibility measurements Ashoori show that the moiré bandwidth at a magic angle is a few tens of meV, larger than what previous calculations found MacDonald ; NamKoshino ; KoshinoFu ; Shiang ; BS1 ; BS2 ; BS3 . 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. LeRoy ). 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.

In this work, we propose a new perspective that relates magic angle to a new type of van Hove singularity (VHS) in the single-particle energy spectrum of TBG. The importance of VHS has already been recognized in our theory of correlated TBG Isobe ; Kozii and related studies Nandkishore ; VH1 ; VH2 ; VH3 ; VH4 ; VH5 ; VH6 using the weak coupling approach. Generally speaking, VHS in two-dimensional systems are associated with saddle points of energy dispersion in $\bm{k}$ space. When a VHS is close to Fermi energy, the increased density of states (DOS) amplifies electron correlation, resulting in various ordering instabilities such as density wave and superconductivity at low temperature Isobe . Indeed, the two recent STS measurements on gate-tunable TBG around magic angle Abhay ; Stevan 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. Lin ). These findings clearly demonstrate the prominent role of VHS in magic-angle TBG.

On the other hand, it is also known from general consideration and previous STS measurements stm1 ; stm2 ; stm3 that VHS are present at all twist angles. It is therefore unclear whether VHS at magic angle is anything special.

We propose the following scenario. As the twist angle decreases below a critical value $\theta_{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 $\theta_{c}$, the saddle point changes from second-order 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^{\circ}$ Dean . 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.

*Ordinary and high-order VHS. * In two-dimensional (2D) electron systems with energy dispersion $E(\bm{k})$, an ordinary VHS with logarithmically diverging DOS occurs at a saddle point ${\bm{k}}_{s}$, determined by

[TABLE]

where $D_{ij}\equiv\frac{1}{2}\partial_{i}\partial_{j}E$ is the $2\times 2$ Hessian matrix of $E$ at ${\bm{k}}_{s}$. Since $D$ is symmetric by definition, we can rotate the axes to diagonalize $D$ and Taylor expand the energy dispersion near ${\bm{k}}_{s}$ as $E-E_{\rm v}=-\alpha p_{x}^{2}+\beta p_{y}^{2}$, where $E_{\rm v}$ is the VHS energy, the momentum $\bm{p}={\bm{k}}-{\bm{k}}_{s}$ is measured from the saddle point, and the coefficients $-\alpha$, $\beta$ are the two eigenvalues of $D$ with $-\alpha\beta=\det D<0$. This dispersion describes two pieces of Fermi contours that intersect at the saddle point ${\bm{k}}_{s}$. It is known from a topological consideration that due to the periodicity of $E(\bm{k})$ on the Brillouin zone (a torus), van Hove saddle points appear quite generally in energy bands of two-dimensional materials vanHove .

We now introduce the notion of high-order VHS by the following conditions at a saddle point

[TABLE]

which imply $\alpha\beta=0$. There exist two types of high-order VHS. The first type corresponds to $\alpha=\beta=0$, i.e., $D=\bm{0}$. The Taylor expansion of $E(\bm{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 $\bm{k}$ point FuSachdev . Such high-order VHS (type-I) is also referred to as multicritical VHS in the recent literature multicritical1 ; multicritical2 ; ruthenate .

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 ($\beta=0$), while the other eigenvalue is nonzero ($\alpha\neq 0$). Then in the Taylor expansion of $E(\bm{k})$, besides the single second-order term $\alpha$ 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 $\bm{K}$ and $\bm{K}^{\prime}$ 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 $\bm{K}$ and $\bm{K}^{\prime}$ valleys are decoupled Neto ; Neto2 ; MacDonald ; Pablo , with energy dispersions denoted by $E_{+}(\bm{k})$ and $E_{-}(\bm{k})$ respectively. Time-reversal symmetry implies $E_{-}(\bm{k})=E_{+}(-\bm{k})$, so it suffices to consider $E_{+}(\bm{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 $\theta=2^{\circ}$ smet reveals the existence of three symmetry-related van Hove saddle points on $\Gamma 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 smet . On the other hand, at $\theta=1.05^{\circ}$, on each $\Gamma 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 $\theta=2^{\circ}$) emerges out of each energy maximum and eventually intersects two Dirac pockets at two new saddle points on opposite sides of $\Gamma M$ KoshinoFu . 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 $\Gamma M$ (denoted by $\Lambda_{0}$) splits into a pair of new ones off $\Gamma M$ (denoted by $\Lambda_{\pm}$). Importantly, the behavior of these VHS points in the vicinity of this transition is solely governed by the local energy dispersion near $\Lambda_{0}$.

We now expand $E(\bm{k})$ near $\Lambda_{0}$ to higher orders:

[TABLE]

where $p_{x}$ ($p_{y}$) is parallel (perpendicular) to $\Gamma M$ line. In the expansion (3), first-order terms vanish because $\Lambda_{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 $\gamma$ term and fourth-order $\kappa$ 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 $\alpha\beta$. Without loss of generality we assume $\alpha>0$ in the following. For $\beta>0$, ${\bm{p}}=\bm{0}$ (i.e., $\Lambda_{0}$) is an ordinary van Hove saddle point with logarithmic divergent DOS. The Fermi contour at the VHS energy consists of two curves that approaches straight lines $p_{y}/p_{x}=\pm\sqrt{\alpha/\beta}$ as $\bm{p}\rightarrow\bm{0}$, intersecting at a finite angle $\varphi=2\arctan(\sqrt{\beta/\alpha})$ as shown in Fig. 2d. This behavior corresponds to the VHS in the calculated band structure of TBG at $\theta=2^{\circ}$. When $\beta\to 0^{+}$, $\varphi\to 0$ so that the two Fermi contours at the VHS energy tend to touch tangentially. On the other hand, for $\beta<0$, $\bm{p}=\bm{0}$ becomes a local energy maximum, while two new saddle points appear at momenta

[TABLE]

whose energy is shifted from $E_{\rm v}$ by $\delta=-\alpha\beta^{2}/\tilde{\gamma}^{2}<0$, where $\tilde{\gamma}\equiv\sqrt{\gamma^{2}+4\alpha\kappa}$. Throughout this manuscript we consider the regime $\gamma^{2}+4\alpha\kappa>0$ so that $\tilde{\gamma}$ is always real and positive, and two split saddle points $\Lambda_{\pm}$ are well-defined. $\Lambda_{+}$ and $\Lambda_{-}$ are a pair of ordinary VHS points related by two-fold rotation, whose Fermi contours are two parabolas $2\alpha(p_{x}\mp\beta/\tilde{\gamma})=(\gamma\pm\tilde{\gamma})p_{y}^{2}$, see Fig. 2f. This behavior corresponds to the VHS in the calculated band structure at $\theta=1.05^{\circ}$. In the limit $\beta\rightarrow 0^{-}$, $\Lambda_{\pm}$ approaches $\Lambda_{0}$.

Our unified description of two regimes of ordinary VHS ($\Lambda_{0}$ and $\Lambda_{\pm}$) using a single local energy dispersion (3) implies that the transition between them corresponds to the sign change of $\beta$. In this process, energy contours at VHS changes from intersecting at one point $\Lambda_{0}$ (Fig. 2d and e) to two points $\Lambda_{\pm}$ (Fig. 2f). Right at $\beta=0$, $\bm{p}=\bm{0}$ becomes a high-order VHS point. The Fermi contour in its vicinity consists of two parabolas $2\alpha p_{x}=(\gamma\pm\tilde{\gamma})p_{y}^{2}$, touching each other tangentially at $\bm{p}=\bm{0}$. In the special case with $\kappa=0$, one of the parabolas becomes 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

[TABLE]

where $z_{\pm}=(\beta/\alpha)\pm\sqrt{(\beta/\alpha)^{2}+\varepsilon}$ and $\varepsilon=\tilde{\gamma}^{2}(E-E_{\rm v})/\alpha^{3}$ is energy deviation from VHS in dimensionless unit. The integration over $\bm{p}$ is extended to infinity. Here ${\rm K}(z)$ is the complete elliptic integral of the first kind, ${\rm sgn}(r)$ is the sign function defined as ${\rm sgn}(r)=-1$ for $r<0$ and $1$ for $r\geq 0$, and $\Theta(r)\equiv\frac{1}{2}[1-{\rm sgn}(-r)]$ is the step function so that $\Theta(0)=0$.

The DOS with different $\beta$ 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 $\beta\neq 0$, the DOS diverges logarithmically at ordinary VHS:

[TABLE]

where $\Lambda=64\alpha\beta^{2}/\tilde{\gamma}^{2}$ is the high-energy cutoff. In the limit $\beta\rightarrow 0$, the prefactor in Eq. (6) diverges as $1/\sqrt{|\beta|}$, indicating a strong increase of DOS as the VHS approaches a high-order one. Right at the transition point $\beta=0$, we find power-law divergent DOS near high-order VHS

[TABLE]

where $C=(2\pi)^{-\frac{5}{2}}\Gamma(\frac{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_{\rm v}$, as reflected by the $\sqrt{2}$ 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. Abhay with finite broadening, as shown in Fig. 1a and Supplemental Material.

At the high-order VHS ($\beta=0$), the dispersion has the following scaling invariance

[TABLE]

with arbitrary $\lambda>0$. Here we set $E_{\rm v}=0$ for simplicity. Including higher order terms or the third and fourth order terms other than $\gamma,\kappa$ 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}]=\frac{1}{2},[p_{y}]=\frac{1}{4}$, with scaling dimension of energy $[E]=1$. From the scaling invariance (8) we immediately obtain the scaling dimension of DOS as $[\rho]=[p_{x}]+[p_{y}]-[E]=-\frac{1}{4}$, the same as Eq. (7). To directly reveal the power-law behavior of experimental data, in Fig. 1b we plot energy deviation $|E-E_{\rm v}|$ and experimental DOS including background contribution both in logarithmic scale. In the log-log plot, two sides of the high-order VHS peak form two parallel lines with the same slope $-1/4$, and the particle-hole asymmetry ratio is 1.32 close to $\sqrt{2}$, which are both consistent with Eq. (7).

Besides $\alpha$ and $\beta$, the Fermi contours in momentum space depend on both $\gamma$ and $\kappa$, while the DOS in energy domain depends only on $\tilde{\gamma}$. This is because the same scaling dimension of $p_{x}$ and $p_{y}^{2}$ allows the nonlinear transform

[TABLE]

under which the dispersion (3) becomes $E-E_{\rm v}=-\alpha\tilde{p}_{x}^{2}+\beta\tilde{p}_{y}^{2}+\tilde{\gamma}^{2}\tilde{p}_{y}^{4}/(4\alpha)$. The nonlinear transform (9) changes the shape of Fermi contours while 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 $\theta_{c}$. The value of $\theta_{c}$ depends on the model we employ. As an illustrative example, we consider the continuum model Neto ; Neto2 ; MacDonald ; KoshinoFu of TBG, where approximate particle-hole symmetry holds, and the Fermi surface topology is solely determined by two dimensionless parameters $g=\lambda u/v,g^{\prime}=\lambda u^{\prime}/v$. Here $u$ and $u^{\prime}$ denote interlayer hoppings, $v$ is the monolayer Dirac velocity and $\lambda$ is the moiré wavelength. We further fix the ratio $g^{\prime}/g=u^{\prime}/u=1.2$ as calculated for corrugated structure at $\theta\sim 1^{\circ}$ KoshinoFu . At given $g\in[1,2]$, we calculate the Taylor coefficients $\alpha,\beta,\gamma,\kappa$ 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 $\gamma^{2}+4\alpha\kappa$ is always positive when $g\in[1,2]$, and across $g_{c}\approx 1.995$, $\beta$ changes sign (Fig. 3b inset) while $\alpha,\gamma,\kappa$ do not. With appropriate and reasonable values of $u$ and $v$, from $g_{c}$ we can obtain the critical twist angle $\theta_{c}\approx 1^{\circ}$, consistent with experiments. Details of continuum model can be found in the Supplemental Material. We may include more ingredients such as lattice relaxation NamKoshino ; Stephen and strain BiZhen 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 $\theta_{c}$ is robust BiZhen .

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 Abhay ; Stevan .

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\chi_{s})^{-1}$ within the random phase approximation (RPA). Since the noninteracting susceptibility $\chi_{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 $\chi(\bm{q})$ with finite momentum $\bm{q}$. There are in total six high-order VHS from two valleys, located at three $\Gamma M$ lines (Fig. 4a). When we consider two high-order VHS along the same $\Gamma M$ line but from opposite valleys, $\chi(\bm{q})$ is sharply enhanced near wave vector $\bm{q}=\bm{Q}$ due to the power-law divergent DOS from both valleys, where $\bm{Q}=2\bm{k}_{s}$ is the momentum separation between the two VHS (Fig. 4a). As a result, when approaching low temperatures, the intervalley susceptibility $\chi(\bm{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_{\rm v}$ is split into two in the presence of intervalley density wave (Fig. 4c).

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 MacDonald 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 Vishwanath . 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 many-body phenomena. Such examples include TBG discussed in the maintext and trilayer graphene on boron nitride discussed in Supplemental Material. The time is coming for creating strong electron correlation by designing single-particle dispersion around VHS.

Acknowledgment. We thank Alex Kerelsky and Abhay Narayan Pasupathy for sharing experimental data with us. This work is supported by DOE Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award DE-SC0010526. LF is partly supported by the David and Lucile Packard Foundation.

Supplemental Material for “Magic of high-order van Hove singularity”

.1 Theoretical fitting of tunneling conductance in twisted bilayer graphene

.1.1 Model-dependent fitting

The dispersion near a van Hove singularity (VHS) is

[TABLE]

and the corresponding density of states (DOS) is

[TABLE]

where $Q\equiv{\rm sgn}(r)\left[{\rm Re}f+\Theta(-r){\rm Im}g\right]$, $\tilde{\gamma}=\sqrt{\gamma^{2}+4\alpha\kappa}$ and

[TABLE]

with $z_{\pm}=r\pm\sqrt{r^{2}+\varepsilon}$. The elliptic integral we use is defined as

[TABLE]

In Eq. (S2) we assume the dispersion (S1) extends to the whole momentum space, while for realistic systems Eq. (S1) holds only for a finite range $|E-E_{\rm v}|<\Omega$. To account for the high-energy cutoff $\Omega$, we introduce a negative background to the DOS expression (S2), and hence tunneling conductance becomes

[TABLE]

with two additional parameters $G_{0}$ and $G_{c}$. Here $G_{0}$ is due to the tunneling matrix element between the sample and the tip in STS experiments, and $G_{c}$ is due to the high-energy cutoff $\Omega$. Notice that Eq. (S5) only applies to energy range $|E-E_{\rm v}|<\Omega$.

In total we have 6 parameters $G_{0},G_{c},E_{\rm v},\eta,r,E_{0}\equiv\alpha^{3}/\tilde{\gamma}^{2}$ in the fitting of tunneling conductance. However, due to the scaling property $Q(\lambda^{2}\varepsilon,\lambda r)=\lambda^{-\frac{1}{2}}Q(\varepsilon,r)$, there are only 5 independent parameters $G_{0},G_{c},E_{\rm v},\eta$ and $w\equiv E_{0}r^{2}{\rm sgn}(r)=\alpha\beta^{2}{\rm sgn}(\beta)/\tilde{\gamma}^{2}$. Here $|w|$ is the low-energy cutoff of the high-order VHS, namely in energy range $|w|<|E-E_{\rm v}|<\Omega$ the VHS can be treated approximately high-order. When $w=0$ the VHS is exactly high-order.

The optimal fitting parameters in least squares fitting are given below in Table. 1 for different samples, and fitting results are shown in Fig. S1. We can treat $w$ as the single indicator of the topological transition of VHS similar to $\beta$ in Eq. (S1). From Table. 1, it can be found

[TABLE]

indicating a topological transition of VHS.

.1.2 Model-independent fitting

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)

[TABLE]

each with three parameters $A,E_{\rm v}$ and $G_{c}$. Among them, $E_{\rm v}$ denotes the energy of conductance peak, and $G_{c}$ is due to high-energy cutoff.

In Fig. S2, we plot conductance data $G$ near the peaks at two different angles $\theta=2.02^{\circ}$ and $1.10^{\circ}$ respectively, where energy is plot in logarithmic scale. When $\theta=2.02^{\circ}$, the conductance is plot in original scale. Since the VHS is ordinary, conductance peak is symmetric and logarithmic, two sides of the peak will collapse into the same line. When $\theta=1.10^{\circ}$, the conductance is 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.

.2 Continuum Model

To determine the value of $\theta_{c}$, we calculate the moiré band structure and track its evolution with twist angle. It is important to note that our calculation of TBG band structure is only to illustrate the Lifshitz transition of VHS and existence of high-order VHS in TBG, thus many ingredients such as strain and in-plane lattice relaxation are neglected.

Our calculation uses a generalized continuum model KoshinoFu

[TABLE]

with $\omega=e^{2\pi i/3}$. Here $\bm{\sigma}$ acts on the two components of massless Dirac fermion at $\bm{K}$ point of each layer. Due to the relative orientation of the two layers, the Dirac spinors on the two layers are rotated by $\pm\theta/2$ respectively under the orthogonal rotation matrix $R_{\pm\theta/2}$. $\bm{p}_{l}=\bm{k}-\bm{K}_{l}$ denotes momentum deviation from the $\bm{K}$ point in layer $l=1,2$. The $2\times 2$ matrix $U(\bm{r})$ denotes interlayer tunneling and is periodic in real space $U(\bm{r}+\lambda\bm{e}_{\pm})=U(\bm{r})$, where $\lambda\bm{e}_{\pm}$ are two primitive vectors of TBG moiré superlattice, and $\lambda=a/(2\sin\frac{\theta}{2})$ is the moiré wavelength at twist angle $\theta$ with graphene lattice constant $a$. $U(\bm{r})$ contains two independent parameters $u$ and $u^{\prime}$ associated with interlayer tunneling within the same sublattice and between the two sublattices respectively. Due to the out-of-plane lattice relaxation, AA regions have larger interlayer distance than AB regions Abhay , resulting in $u^{\prime}<u$ in improved continuum model of TBG band structure KoshinoFu .

At small $\theta$, Fermi surface topology and hence VHS in the continuum model (S8) are solely determined by two dimensionless parameters $g\equiv\lambda u/v$ and $g^{\prime}\equiv\lambda u^{\prime}/v$ ,where $v$ is the Dirac velocity and $\lambda$ is the moiré wavelength. This can be seen by following scaling relation

[TABLE]

Dimensionless parameters $g,g^{\prime}$ increase with the moiré wavelength as $\theta$ decreases; they also increase with interlayer tunneling under pressure. For simplicity, in the following we fix the ratio $g^{\prime}/g=u^{\prime}/u=1.2$ as calculated for relaxed structure at $\theta\sim 1^{\circ}$ KoshinoFu and study VHS evolution with the single parameter $g$. Since electron and hole sides are approximately symmetric, we focus on VHS in the hole side.

At $g=1$ there are three ordinary VHS on $\Gamma M$ (Fig. 2a of maintext), while at $g=2$ there are six ordinary VHS on two sides of $\Gamma M$ (Fig. 2c of maintext). According to our analysis above, this implies a topological transition of VHS occurs at $g=g_{c}\in(1,2)$. To determine the value $g_{c}$, we calculate the Taylor coefficients $\alpha,\beta,\gamma,\kappa$ by numerical derivatives of energy dispersion with respect to momentum components $k_{x},k_{y}$ at VHS point. As shown in Fig. 3b of maintext, $\alpha,\beta,\gamma,|\kappa|$ all decrease as $g$ increases, and across $g_{c}\approx 1.995$, $\beta$ changes sign (Fig. 3b inset) while $\alpha,\gamma$ stay positive and $\kappa$ stays negative. It can be verified that through the whole process $\gamma^{2}+4\alpha\kappa$ is positive so that the asymptotic behavior is always saddle point.

At $g<g_{c}$, we find that over an extended parameter range $g\in[1.6,1.8]$, the van Hove filling is close to two electrons per unit cell (where correlated electron phenomena are observed), in agreement with the observed sign change of Hall coefficient under doping and pressure Dean . Importantly, in this range $0<\beta\ll\alpha$ and VHS can still be treated approximately as high-order. At small $g-g_{c}>0$, the bandwidth decreases rapidly; the band structure calculated from continuum model becomes rather complex above $g_{c}$ Kasra . It is possible that additional effects not captured in Eq. (S8), such as in-plane lattice distortion, become important in this regime NamKoshino ; Stephen ; BiZhen . Nonetheless, the existence of topological transition of VHS as a function of $\theta$ is a robust feature.

To translate the parameter $g$ into the actual twist angle requires the knowledge of Dirac velocity $v$ and interlayer hoppings $u,u^{\prime}$, which are subject to some uncertainty. It is known that Coulomb interaction increases the Dirac velocity substantially at low energy Geim ; Elias ; Schliemann . It is recently proposed that this increased velocity may account for the discrepancy between the moiré bandwidth found in STS measurements and in previous calculations Abhay . When we use an increased Dirac velocity $va^{-1}=2.41$eV, 13% larger than DFT value and adopt the interlayer hopping parameters $u=79.7$meV, $u^{\prime}=97.5$meV from Ref. KoshinoFu , we find $\theta_{c}=0.95^{\circ}$, and the experimentally established range of magic angle $\theta_{\rm exp}\in[1^{\circ},1.2^{\circ}]$ corresponds to the range $g_{\rm exp}\in[1.6,1.9]$ where $0<\beta\ll\alpha$ and VHS is approximately high-order. We can also apply pressure to change interlayer couplings $u,u^{\prime}$ and hence $g,g^{\prime}$. In this way, high-order VHS can be achieved at larger twist angles.

.2.1 Total density of states from continuum model

We can also use continuum model instaed of effective model (S1) to compute the total density of states. In Fig. S3 we plot the Fermi contours and corresponding total density of states near van Hove singularity in continuum model with different coupling constant $g$.

.2.2 Local density of states from continuum model

In the maintext and also previous sections, we calculate the total density of states which is the averaged DOS over the whole sample. In this section we show that at magic angle, even local density of states (LDOS) also exhibits behaviors of high-order VHS, both in AA and AB regions.

In Fig. S4, by continuum model we calculate LDOS in AA and AB regions of TBG at magic angle as functions of energy, where both LDOS and energy are plot in log scale and $E_{\rm v}$ denotes the energy of conductance peak. Red and blue colors denote peaks at positive and negative energy sides respectively. We find that LDOS in both AA and AB regions show asymmetric power-law divergent peaks, while LDOS in AA region is much larger than LDOS in AB region, which is consistent with experimental data in arXiv:1812.08776 by A. Kerelsky et. al.

.3 Mean-field Hamiltonian of interacting high-order VHS

With finite intervalley density wave order parameter $\Delta$, the low-energy Hamiltonian of two high-order VHS on the same $\Gamma M$ line is

[TABLE]

where Pauli matrices $\bm{\tau}$ act on the valley indices. From this Hamiltonian Fig. 4c of maintext is plotted, and the order parameter is $\Delta=0.01E_{0}$ with broadening $\eta=5\times 10^{-3}E_{0}$.

.4 Trilayer graphene on boron nitride

The low-energy physics of ABC-stacked trilayer graphene (ABC-TLG) is described by states at $A_{1}$ and $B_{3}$ sites near $\pm\bm{K}$ points of Brillouin zone shown in Fig. S5c. The point group at $\bm{K}$ point is $D_{3}$ with two generators $C_{3z}$ and $C_{2x}$. Since $A_{1}$ and $B_{3}$ are at rotation center $O$ and equivalent points, Bloch waves with momentum $\bm{K}$ at $A_{1}$ and $B_{3}$ will carry zero angular momentum under $C_{3z}$ and furnish 1D representations $A_{1}$ and $A_{2}$ of $D_{3}$. As a result, the band spectrum at $\bm{K}$ point is gapped. In the basis of $A_{1}$ and $B_{3}$, the 2 by 2 effective $k\cdot p$ Hamiltonian near $+\bm{K}$ point reads

[TABLE]

up to the third order in $\bm{p}$, where $p_{\pm}=p_{x}+ip_{y}$ and $p=|\bm{p}|=|p_{\pm}|$. Here the parameters are (unit: eV)

[TABLE]

Denote the lattice constants of ABC-TLG and hexagonal boron nitride (h-BN) as $a_{G}$ and $a_{BN}$ respectively. For simplicity, we assume there exist two coprime positive integers $m,n$ such that $ma_{G}=na_{BN}=L_{\rm M}$ is the superlattice constant. Denote $\delta=(a_{G}-a_{BN})/a_{BN}$ as the lattice constant mismatch, then moiré superlattice constant is $L_{\rm M}=a_{G}/|\delta|$.

The substrate h-BN will create a periodic potential on the moiré scale

[TABLE]

where $\bm{G}_{j}=4\pi/(\sqrt{3}L_{\rm M})(-\sin\frac{2\pi j}{3},\cos\frac{2\pi j}{3})$, and $\xi=\pm$ denotes different alignments of TLG and h-BN.

.4.1 High-order van Hove singularity

In this section we consider alignment scheme $\xi=-1$ where only B sublattice is affected by periodic potential of h-BN. The moiré parameters are

[TABLE]

We focus on the valence bands in Fig. S6. Near $K^{\prime}$ point, we can expand the dispersion of the upper valence band as

[TABLE]

where the coefficients $\alpha,\gamma$ and $\kappa$ are functions of electric field $V$. Near transition field $V_{c}=0.06$eV, the corresponding Fermi contours are shown in Fig. S7. We find that when $V=V_{c}$, $K^{\prime}$ point is a type-I high-order VHS where second order coefficients all vanish and three Fermi surfaces merge together.

.4.2 Effective tight-binding model

We use the following tight-binding model to fit the two trivial bands in the continuum model spectrum

[TABLE]

For moiré parameters in alignment scheme $\xi=+1$ where only A sublattice is affected by periodic potential of h-BN

[TABLE]

the corresponding fitting parameters of the lowest conduction bands are

[TABLE]

Since $\phi$ is close to $\pi/2$, the nearest neighbor term mostly contributes to the warping effect. Fittings both along specific path and in the entire MBZ are shown in Figs. S8 and S9 respectively.