Introduction

Twisted bilayer graphene (TBG) is a two-dimensional correlated electronic system, which exhibits superconductivity1,2,3 and correlated phases4,5,6,7,8,9,10,11,12,13. The focus of our work is the analysis of a cascade of phase transitions near-integer fillings ∣n∣ = 1, 2, 3, 4, detected in STM and electronic compressibility measurements5,14,15 (panels (a)–(e) in Fig. 1). Compressibility measurements show sharp seesaw features of dμ/dn near integer ∣n∣, and STM data show that around each of these n a peak in the density of states splits, and one of its components appears on the other side of the Fermi level. For the interpretation, the authors of ref. 14 adopted a strong-coupling approach and associated the observed STM peaks with narrow sub-bands. They argued that at each transition one sub-band crosses the Fermi level, moves away from it, and becomes incoherent. The authors of ref. 15 interpreted compressibility data within a moderate coupling scenario of a 4-fold spin/isospin degenerate band and argued that the cascade can be understood as a series of interaction-driven transitions. They conjectured that at, e.g., electronic doping one of the bands gets completely filled at each transition, while the occupation of the remaining ones gets depleted; mirror symmetric behavior holds for hole doping.

Fig. 1: The cascade of electronic transitions in twisted bilayer graphene.
figure 1

a–d the proposed splitting of the initially 4-fold (spin and valley) degenerate vH peak upon raising the electron filling n (panels (a), (b)) and hole filling (panels (c), (d)), and the corresponding STM data for twist angle of 1.06∘, reproduced with permission from the authors of ref. 14. The van Hove peaks in the conduction (valence) bands are labeled by I(II), and subscripts a,b,c,d label the 4 peak components. e Experimental data for inverse compressibility with seesaw features, interpreted as a cascade of phase transitions. Reproduced with permission from the authors of ref. 15. The data are for twist angle θ = 1.13∘, f The schematic phase diagram of TBG upon electron or hole doping, based on comparison between the theory and the STM data. Within our model, we obtained SU(3) × U(1) (U(1) × SU(3)) symmetry for both electron and hole doping at 1 < ∣n∣ < 2 (3 < ∣n∣ < 4), corresponding to 3-1 (1-3) splitting, but it is reduced to SU(2) × U(1) × U(1) (U(1) × U(1) × SU(2)) if there is an intermediate phase with 2-1-1 (1-1-2) splitting, as STM data for hole doping likely indicate.

Full size image

In this communication we propose the scenario in which the cascade of transitions is caused by the development of particle-hole orders, like in ref. 15, but we specifically identify the STM peaks with van Hove (vH) singularities. We argue that the components of the initially 4-fold degenerate vH peak move through the Fermi level one by one, but remain close to it. The split peaks recombine into a single 4-fold peak at ∣n∣ ≲ 4, when electronic order vanishes. Our scenario is illustrated in panels (a), (c), and (f) in Fig. 1.

A cascade of transitions has been observed near van Hove doping in less correlated non-twisted Bernal bilayer (BBG) and rhombohedral trilayer graphene (RTG)16,17,18,19. We show that our vH scenario equally explains the sequence of transitions in these materials. We believe that the similarity between the ordered states and electronic reconstruction in BBG/RTG and in TBG supports a moderate coupling vH-based approach. We emphasize, however, that we use this approach specifically to describe the cascade of phase transitions with doping. A strong-coupling approach is needed for explaining the insulating behavior of TBG near-integer fillings.

We further emphasize that (i) vH peaks have been observed in TBG at different twisting angles6,20, (ii) are present in the electronic dispersion, obtained in first-principle calculations, and in the one renormalized by the interaction, even if the bottom of the dispersion moves away from Dirac points21, and (iii) the cascade of transitions, observed in magic-angle twisted trilayer graphene, has been argued to be triggered by vH peaks, at least at high displacement fields22. The vH scenario has been also discussed in context of chiral density wave and superconductivity in TBG (see e.g., refs. 23,24,25,26,27).

The summary of our results for TBG is presented in Fig. 1 along with the experimental data from refs. 14,15. We label the vH peaks in conduction (valence) bands in by I (II) and label peak components by a,b,c,d. Our interpretation of the STM data from ref. 14 for electron doping (panel (b)) is the following: as the system moves away from charge neutrality, the 4-fold degenerate peak Ia+b+c+d approaches μ from above, and at n ≤ 1, splits in a 3-1 fashion: three components Ia+b+c stay above μ, and one component, Id, jumps to below the Fermi level, but remains close to it. At n ≤ 2, the three components again come close to μ, and the vH peak Ia+b+c splits in 2-2 fashion, such that Ia+b moves back, while Ic jumps to below the Fermi level and merges with Id into Ic+d. At n ≲ 3, Ia+b splits and Ib jumps to below the Fermi level and merges with Ic+d into Ib+c+d. Finally, at n ≲ 4, the last component Ia jumps across the Fermi level and merges with three other components into 4-fold degenerate Ia+b+c+d. For hole doping the overall evolution is the same, but the data seem to show a more gradual behavior: the components of peak II in panel (d) cross the Fermi level one-by-one, indicating the presence of an intermediate state between 3-1 and 1-3 ones.

We consider these data as evidence that once the 4-fold degenerate vH peak gets close to μ at ∣n∣ ~ 1 (peak I for n > 0 and peak II for n < 0), the system develops a vH-induced particle-hole order. The order exists between ∣n∣ ≲ 1 and ∣n∣ ≲ 4 and reconstructs the fermionic spectra, pushing vH peaks in some bands above μ and in other band(s) below μ. The structure of the order changes near ∣n∣ = 2 and ∣n∣ = 3, via first-order transitions, and this changes the splitting of the vH peak and simultaneously gives rise to sharp changes in the compressibility (see Fig. 1e).

Here, we present the theoretical description of this scenario within the model of interacting itinerant electrons whose band structure has a vH singularity near μ. We use as an input our earlier result28,29 that the increased density of states near the vH singularity enables a spontaneous symmetry breaking in spin and valley spaces (for particle-hole orders with zero transferred momentum electronic DOS has to exceed a threshold in order to satisfy the Stoner-type criterion). For a model with intra-site (Hubbard) and assisted hopping interactions within a given hexagon30, we found 15 particle-hole order parameters, for which the couplings are attractive, near-equal, and larger than for other order parameters. They describe intra-valley order at zero momentum (Q = 0) and inter-valley density waves (Q ≠ 0). These 15 order parameters are described by 4 × 4 matrices in spin and valley spaces, specified by spin σ and valley isospin τ and form the adjoint representation of the SU(4) group. For shortness, we call an order in (σ, τ) space a spin/isospin order. An SU(4)-symmetric order parameter manifold and the interplay between spin and isospin orders have been recently discussed in refs. 30,31,32,33,34,35,36,37—these works provide additional motivation for us. Another input for our analysis are band structure calculations38,39, which reported the pinning of the vH singularity to the chemical potential over the range of n.

We derive and analyze the Landau free energy for an SU(4)-symmetric fermionic model. We argue that there are three sets of ordered states, which split the 4-fold degenerate vH peak in three different ways. The first vH-induced instability splits the vH peak in a 3-1 fashion, with peaks for 3 degenerate bands shifting towards charge neutrality, and the remaining peak moving to below μ. As the magnitude of a spin/isospin order increases, the system undergoes a transition into a different ordered state, which splits the vH peak into a 2-2 fashion. The transition is first-order in our model, but in reality may occur via an intermediate phase with 2-1-1 splitting. At ∣n∣ ~ 3 the magnitude of the order starts decreasing as some of vH peak components move further away from μ, and the system behavior goes in reverse—first the system undergoes a transition into an ordered state, which gives rise to 1-3 vH peak splitting, again either via a direct first-order transition, or via an intermediate phase with 1-1-2 splitting, and then, at even larger ∣n∣ ≤ 4, the order vanishes, and all 4 vH peak components merge into a single vH peak below μ.

We use the same approach for BBG and RTG. The bands structures and Fermi surfaces of BBG and RTG are very similar, and we model both systems by an effective patch model of fermions, located in the vicinity of K and \({{{\bf{K}}}}^{\prime}\) points in the BZ. We find that the 15 leading instabilities are analogous to TBG: towards valley polarization and intra-valley spin order (both with Q = 0) and towards inter-valley charge and spin density wave orders with \({{{\bf{Q}}}}={{{\bf{K}}}}-{{{\bf{K}}}}^{\prime}\). We find that for a Hubbard interaction, these orders are described by the same SU(4)-symmetric Landau free energy functional, Eqs. (1) and (2), as in TBG. Like in TBG, the first vH-induced transition is into a state with valley polarization and ferromagnetism in a single valley. This order gives rise to 3–1 splitting, which in the case of BBG/RTG gives rise to one larger and three smaller Fermi pockets. This splitting is analogous to the one observed in the IF1 state in the notations of ref. 16. The subsequent transition upon doping is into a state with either pure valley charge order or ferromagnetic order in both valleys. This state gives rise to 2–2 splitting, which in BBG/RTG gives rise to two larger and two smaller Fermi pockets. This is analogous to PIP2 state16. A potential intermediate state with 2–1–1 splitting is analogous to PIP1 state in ref. 16.

Results

Cascade of transitions in TBG

Band structure calculations show that there are eight bands within the flat-band regime of TBG, accounting for two spin projections, two valley degrees of freedom from the original graphene layers, and two sublattices of the moiré superlattice. Four bands are with upward and four with downward dispersion, merging at Dirac points K and \({{{\bf{K}}}}^{\prime}\). Upon electron (hole) doping the chemical potential moves up (down), simultaneously changing the filling of four bands. Each band displays a vH singularity. The vH singularities for four conduction (four valence) bands are at the same energy. It was argued that strong-coupling renormalizations may shift the minimum of electron band to the Γ point21,38, but vH singularities remain even for the renormalized dispersion38,40.

We study the cascade of phase transitions by analyzing the Landau free energy for the ordered phases of fermions with vH singularity near μ. The order parameters are expectation values of fermionic bilinears, and the free energy can be obtained by departing from a microscopic model of vH fermions with 4-fermion Hubbard and assisted hopping interactions30 and integrating out fermions after performing a Hubbard-Stratonovich transformation. Alternatively, one can write down the Landau free energy solely based on symmetries and fix parameters phenomenologically through comparison with experiments. In an earlier study29 we found that out of a large number of possible fermionic bilinears (143 in the 6-patch vH model and even larger number in 12-patch model) there are 15, for which the couplings are attractive and the largest by magnitude. The set of 15 is composed of two subsets of 7 and 8 bilinears with a single coupling within each subset. 7 bilinears with coupling λ7 are intra-valley with transferred momentum Q = 0, and 8 with coupling λ8 are inter-valley with a finite Q (Kekule-type states considered in ref. 41).

The couplings λ7 and λ8 are not identical, but are close to each other28. In our analysis we treat λ7 and λ8 as equal, in which case the 15 bilinears form an adjoint representation of SU(4). We checked that the cascade of transitions and the sequence of vH peak splitting is the same in the model with only 7 bilinears (the case λ7 > λ8). In the model with 8 bilinears there is a single ordered phase and no cascade.

For the SU(4) case, the Landau free energy up to fourth order is29 (note, that the expression here uses a slightly different definition of prefactors in the free energy than the one in ref. 29)

$${{{\mathcal{F}}}}=-\frac{\alpha }{2}{{{\rm{Tr}}}}\left[{\hat{\Phi }}^{2}\right]+\frac{\gamma }{3}{{{\rm{Tr}}}}\left[{\hat{\Phi }}^{3}\right]+\frac{\beta }{4}{{{\rm{Tr}}}}\left[{\hat{\Phi }}^{4}\right]+\frac{\beta ^{\prime} }{4}{{{\rm{Tr}}}}{\left[{\hat{\Phi }}^{2}\right]}^{2},$$
(1)

where \(\hat{\Phi }=\mathop{\sum }\nolimits_{j = 1}^{15}{\phi }_{j}{T}^{j}\), Tj are generators of SU(4), and ϕj ~ f†Tjf are fermionic bilinears, which we treat as Hubbard-Stratonovich fields (f and f(†) are operators of electrons near vH points). The term \(\beta ^{\prime}\) does not appear within Hubbard-Stratonovich but is allowed by symmetry, and we keep it for generality.

By construction, \(\hat{\Phi }\) can be represented by a traceless matrix42. In the diagonal basis

$$\hat{\Phi }={{{\rm{diag}}}}({\lambda }_{1},{\lambda }_{2},{\lambda }_{3},-({\lambda }_{1}+{\lambda }_{2}+{\lambda }_{3})),$$
(2)

and the free energy is

$$\begin{array}{ll}{{{\mathcal{F}}}}\,=-\frac{\alpha }{2}\left(\mathop{\sum }\limits_{j}^{3}{\lambda }_{j}^{2}+{\left(\mathop{\sum }\limits_{j}^{3}{\lambda }_{j}\right)}^{2}\right)+\frac{\gamma }{3}\left(\mathop{\sum }\limits_{j}^{3}{\lambda }_{j}^{3}-{\left(\mathop{\sum }\limits_{j}^{3}{\lambda }_{j}\right)}^{3}\right)\\ \qquad\quad+\,\,\frac{\beta }{4}\left(\mathop{\sum }\limits_{j}^{3}{\lambda }_{j}^{4}+{\left(\mathop{\sum }\limits_{j}^{3}{\lambda }_{j}\right)}^{4}\right)+\frac{\beta ^{\prime} }{4}{\left(\mathop{\sum }\limits_{j}^{3}{\lambda }_{j}^{2}+{\left(\mathop{\sum }\limits_{j}^{3}{\lambda }_{j}\right)}^{2}\right)}^{2}.\end{array}$$
(3)

At γ = 0, the order develops continuously when α changes sign and becomes positive. At a finite γ, the transition is necessarily first order and occurs already when α is negative. Below we restrict to α > 0, when the order is already finite and also set β > 0, consistent with the Hubbard-Stratonovich analysis and the calculation of α, β for the tight-binding model near a vH singularity29. We discuss the behavior of γ below and in Supplementary Discussion V. Minimizing \({{{\mathcal{F}}}}\) with respect to λj (j = 1, 2, 3), we find three solutions (up to permutations of λj): (i) λ1 = λ2 = λ3; (ii) λ1 = λ2 = − λ3; (iii) λ1 = λ2 ≠ λ3 (see Supplementary Discussion I for details). For the first solution, \(\hat{\Phi }={{{\rm{diag}}}}(\lambda ,\lambda ,\lambda ,-3\lambda )\), and the broken symmetry is described by the coset SU(4)/[SU(3) × U(1)], where SU(3) corresponds to the transformation within the subset of the first three components of \(\hat{\Phi }\), and U(1) to a rotation of the last component relative to the other three. The ordered states in terms of expectation values of fermionic bilinears 〈ϕi〉 are mixtures of spin/isospin order with particular ratios of spin and isospin components29. For example, a pure intra-valley order is a part of this set, but a pure inter-valley order is not. The order parameter manifold has 15 – 8 – 1 = 6 Goldstone modes. The feedback of this order on fermions is 3-1 or 1-3 splitting of vH peaks, depending on the sign of γ. For the second solution, \(\hat{\Phi }={{{\rm{diag}}}}(\lambda ,\lambda ,-\lambda ,-\lambda )\), and the broken symmetry is SU(4)/[SU(2) × SU(2) × U(1)], where the two SU(2)’s correspond to rotations within the subsets of the first two and the last two components of \(\hat{\Phi }\), and U(1) corresponds to a rotation of one subset relative to the other. The ordered states in terms of 〈ϕi〉 include pure spin and isospin orders, e.g., intra-valley ferromagnetism and valley polarization, and various inter-valley density waves29. This manifold has 15 – 6 – 1 = 8 Goldstone modes. The feedback from such order on fermions leads to 2-2 splitting of the vH peaks. Finally, the third solution describes a mixed state with \(\hat{\Phi }={{{\rm{diag}}}}(\lambda ,\lambda ,-{\lambda }_{3},-2\lambda +{\lambda }_{3})\) and broken symmetry SU(4)/[SU(2) × U(1) × U(1)]. The order parameter manifold contains 15 – 3 – 1 – 1 = 10 Goldstone modes. The feedback on fermions leads to 2-1-1 or 1-1-2 splitting of vH peaks.

The values of λ and the free energies for the three states, \({F}_{l}=\frac{{\alpha }^{2}}{\beta }{f}_{l}(x,y)\), are functions of \(x=\gamma /\sqrt{\alpha \beta }\) and \(y=\beta ^{\prime} /\beta\):

$$\begin{array}{ll}(i){\lambda }_{{{{\rm{i}}}}}=\sqrt{\frac{\alpha }{\beta }}\frac{| x| \pm \sqrt{{x}^{2}+7+12y}}{7+12y},{f}_{{{{\rm{i}}}}}(x,y)\\\qquad \quad=-\frac{{\left(| x| \pm \sqrt{7+{x}^{2}+12y}\right)}^{2}[3(7+12y)+2| x| (| x| \pm \sqrt{7+{x}^{2}+12y})]}{{(7+12y)}^{3}}\\ (ii){\lambda }_{{{{\rm{ii}}}}}=\sqrt{\frac{\alpha }{\beta }}\frac{1}{\sqrt{1+4y}},{f}_{{{{\rm{ii}}}}}(x,y)=-\frac{1}{1+4y}\\ (iii){\lambda }_{{{{\rm{iii}}}}}=\sqrt{\frac{\alpha }{\beta }}x,{\lambda }_{{{{\rm{iii}}}},3}=\sqrt{\frac{\alpha }{\beta }}\left(\sqrt{\frac{1-{x}^{2}(1+4y)}{1+2y}}-| x| \right),{f}_{{{{\rm{iii}}}}}(x,y)\\\qquad \qquad=-\frac{1+2{x}^{2}-{x}^{4}(1+4y)}{2(1+2y)}.\end{array}$$
(4)

The solution (iii) exists for \(| x| \le 1/\sqrt{1+4y}\) and the ± sign is for positive/negative γ. We plot the free energy prefactors fl(x, y) in Fig. 2.

Fig. 2: Landau free energy.
figure 2

The functions fl(x, y) from Eq. (4), with l = i, ii, iii, shown as functions of \(x=\gamma /\sqrt{\alpha \beta }\) for two values of \(y=\beta ^{\prime} /\beta\): y = 0 (panel (a)) and y = 1 (panel (b)). The free energies are Fl = (α2/β)fl(x, y), hence the smallest fl(x, y) determines the ground state. The states (i) and (ii) are the ones for which the vH peaks split in 3-1 (1-3) and 2-2 fashion, respectively. The state (iii) is an intermediate state with 2-1-1 (1-1-2) splitting. This intermediate state does not appear as a ground state in our model for all y, but its energy is close to those of (i) and (ii) near critical x of the first-order transition between the two, and it can potentially become a ground state around this x if we move away from SU(4)-symmetric model by, e.g., including interaction terms with inter-valley scattering. The relation between x and ∣n∣ is shown at the bottom.

Full size image

We see that at large ∣x∣ the ground state configuration is state (i) while for small ∣x∣ it is state (ii). There is a direct first-order transition between states (i) and (ii) at some intermediate ∣x∣ = xcr. We expect that α > 0 between 1 < ∣n∣ < 4, where the vH peak remains near the chemical potential, and argue that γ changes sign from positive to negative as ∣n∣ increases, because the sign of γ is different when the bands are empty and when they are filled. As a result x evolves from a large positive value to a large negative one via zero upon increasing ∣n∣. As small α corresponds to large ∣x∣, when the order first emerges, the system moves into state (i), and the components of the vH peak split in 1-3 fashion for positive γ. As α increases, ∣x∣ decreases and eventually reaches xcr, where the system undergoes a first-order transition into the ordered state (ii), for which the splitting of the components of the vH peak is 2-2. At larger n, γ changes sign and its magnitude increases, while α starts decreasing. As a result, ∣x∣ increases. When it reaches xcr, the system undergoes another first-order transition into the state, which gives rise to 3-1 splitting of the components of the vH peak. Eventually the order disappears and all 4 components of the vH peak recombine into a single peak. We also note that while the intermediate state (iii) is not the ground state for any x and y, its free energy Fiii is only slightly larger than Fi and Fii at ∣x∣ near xcr. This is particularly so at large y (at \({x}_{{{{\rm{cr}}}}}\approx \sqrt{3/16y}\), Fiii is larger than Fi = Fii by (α2/β)1/(16y)2). Thus, it seems possible that the intermediate state (iii) will become the ground state once we move away from an SU(4)-symmetric model by, e.g., including interaction terms with inter-valley scattering. Such terms are small, but finite in TBG30,43,44. If the transition from (i) to (ii) is via the intermediate phase (iii), there is a range of ∣n∣ where the splitting of the vH peak components is 2-1-1 or 1-1-2, again depending on the sign of γ. Some indications of 2-1-1 and 1-1-2 splitting have been found in STM for hole-doped samples5,14.

Cascade of transitions in BBG and RTG

The same analysis can be applied to study the cascade of phase transitions in BBG and RTG. In both systems, application of an electric field opens a gap between conduction and valence bands and flattens the fermionic dispersion near Dirac K and \({{{\bf{K}}}}^{\prime}\) points45. Near charge neutrality, this creates small Fermi pockets, three near K and three near \({{{\bf{K}}}}^{\prime}\). Upon doping, pockets merge at vH fillings and eventually transform into one larger pocket near K and one near \({{{\bf{K}}}}^{\prime}\) (refs. 46,47,48,49,50,51). We consider the full 6-pocket model in Supplementary Discussion XI and here illustrate the behavior using a simplified model of fermions in two patches near K and \({{{\bf{K}}}}^{\prime}\) with Hubbard intra-patch and inter-patch density-density interaction. In this model, electronic instabilities towards valley polarization, intra-valley ferromagnetism, and inter-valley spin and charge order all occur at the same critical coupling λ. These 15 bilinears then form an adjoint representation of SU(4) and are described by the same Landau free energy functional as in Eq. (1). The cascade of transitions in BBG and RTG then matches the one in TBG with the only difference that some pockets may sink below the Fermi level (see Fig. 3). In a more realistic 6-patch model, the coupling for 7 C3 symmetry preserving order parameters with Q = 0 is not the same as for order parameters with momenta Q close to \({{{\bf{K}}}}-{{{\bf{K}}}}^{\prime}\). The sequence of transition and the Fermi surface reconstruction remain the same as in the 2-patch model if the order develops with Q = 0.

Fig. 3: The cascade of transitions in BBG/RTG.
figure 3

The notations—the same as in refs. 16,19, are as follows: IF1—isospin ferromagnet, PIP1— partially isospin-polarized phase with one large Fermi surface and one small, PIP2—partially isospin-polarized phase with two large Fermi surfaces and two small, Sym4—a symmetric phase with 4 identical large Fermi surfaces (one per isospin), Sym12—a symmetric phase with 12 identical small Fermi surfaces (three per isospin). The states in panels a–e are symmetric 4-fold degenerate, 1–3, 1–1–2, 2–2, and again symmetric 4-fold degenerate, correspondingly. The 3−1 and 2−1−1 states have not been detected in ref. 16 and are not shown. The small pockets in 1−3 (IF1) are assumed to sink below the Fermi level. The symmetry between three small pockets in panels c and d may be broken by subleading interactions, leaving only one small pocket, as the data in ref. 16 indicate.

Full size image

Comparison with experiments on TBG

In our proposed vH scenario, spin/isospin order develops at ∣n∣ ≲ 1, when the 4-fold degenerate vH peak approaches the Fermi level, and persists up to ∣n∣ ≲ 4. In this range of n, the vH peak splits, but according to STM data, its components are still located near the Fermi energy, i.e., the enhancement of the DOS near the Fermi level persists. At larger ∣n∣, the vH peak again becomes 4-fold degenerate and moves away from the Fermi level. The evolution of spin/isospin order and of its feedback on the components of the vH peak is governed in our theory by the relative strength of the prefactor of the cubic term in the Landau free energy (specifically, by \(x=\gamma /\sqrt{\alpha \beta }\)). This prefactor is expressed via a convolution of three fermionic propagators and vanishes for particle-hole symmetry around the Fermi surface. In the absence of such symmetry, γ is non-zero. We conjecture that x is positive near n = 1 passes through zero at 2 < n < 3, and becomes negative at larger n (see Supplementary Discussion V for more discussion on this). We then end up with the phase diagram in Fig. 1f. There are two phase transitions between disordered and ordered states at ∣n∣ ≲ 1, and ∣n∣ ≲ 4, and two transitions between different ordered phases at ∣n∣ ≲ 2 and ∣n∣ ≲ 3. Specifically, within our theory the sequence for the symmetry-breaking pattern and the degeneracy of the vH peak is:

$$\begin{array}{l}n\,\,\lesssim \,\,1:{{{\rm{SU}}}}(4)\to {{{\rm{SU}}}}(3)\times {{{\rm{U}}}}(1):(4,0)\to (3,1)\\ n\,\,\lesssim \,\,2:{{{\rm{SU}}}}(3)\times {{{\rm{U}}}}(1)\to {{{\rm{SU}}}}(2)\times {{{\rm{SU}}}}(2)\times {{{\rm{U}}}}(1):(3,1)\to (2,2)\\ n\,\,\lesssim\,\, 3:{{{\rm{SU}}}}(2)\times {{{\rm{SU}}}}(2)\times {{{\rm{U}}}}(1)\to {{{\rm{U}}}}(1)\times {{{\rm{SU}}}}(3):(2,2)\to (1,3)\\ n\,\,\lesssim \,\,4:{{{\rm{U}}}}(1)\times {{{\rm{SU}}}}(3)\to {{{\rm{SU}}}}(4):(1,3)\to (0,4)\end{array}$$
(5)

where a and b in (a, b) indicate the number of vH peaks above and below μ for the case of electron doping. For hole doping the sequence is identical, except a and b in (a, b) are interchanged. If the transformations (3, 1) → (2, 2) and (2, 2) → (1, 3) occur via an intermediate phase (c), each of the two first-order transitions around ∣n∣ = 2 is replaced by two second-order transitions with the intermediate structure of vH peaks (2, 1, 1) and (1, 1, 2).

The theoretical phase diagram agrees with the STM results5,14 (Fig. 1b, d), including fine details, lending support to our theory. Note, that there is no symmetry of the phase diagram with respect to ∣n∣ = 2, i.e., the transitions at ∣n∣ ≲ 1 and n ≲ 3 are different ones (there is an approximate symmetry with respect to ∣n∣ = 2.5). The theory also explains the seesaw behavior of electron compressibility, reported in ref. 15 and shown in Fig. 1e. Our reasoning is the following. As doping increases and the system approaches one of transitions from the cascade, the inverse compressibility dμ/dn decreases as the n-times degenerate vH peak approaches the Fermi level, where dμ/dn = 0 (n = 4, 3, 2, depending on the number of the transition in the cascade). After a new order develops, one peak component crosses the Fermi level, while the other (n − 1) components move back from the Fermi level. As all vH peaks move away from the Fermi level in a first-order transition, dμ/dn jumps to a larger value. As doping increases further towards the next transition from the cascade, the (n − 1)—times degenerate vH peak approaches the Fermi level, and dμ/dn again decreases towards zero. Then the new order develops, one peak component crosses the Fermi level, while the other (n − 2) components move back from it, and dμ/dn again jumps to a higher value. This gives rise to seesaw structure of the inverse compressibility (see Supplementary Discussion V for an example calculation of dμ/dn for one transition of the cascade). As all four vH peak components remain close to the Fermi level, all four contribute to the evolution of dμ/dn between the transitions. This is consistent with a weak dependence of the slope of dμ/dn on the number of a transition in the cascade.

Comparison with experiments on BBG/RTG

Measurements of inverse electronic compressibility and magnetoresistance in BBG16,17 and RTG19 at a finite displacement field revealed a cascade of transitions upon hole or electron doping. The fermionic structure of the two materials is almost identical, and for definiteness we focus on hole-doped BBG. Near charge neutrality, the system is in the valley/spin symmetric state (labeled Sym12 in ref. 16 and in Fig. 3a) with twelve Fermi pockets: three spin-degenerate ones for each valley. At large enough doping, the triad of pockets for each valley and spin transforms into a single larger pocket, leaving four pockets, again valley and spin symmetric (Sym4 state in ref. 16 and in Fig. 3e). The cascade of transitions happens in between these two limits, when the system develops particle-hole order that breaks valley and/or spin symmetry. We show the sequence of transitions in the cascade in the 2-patch model in Fig. 3.

The authors of refs. 16,17 detected the symmetric three intermediate phases, which they labeled IF1, PIP1, and PIP2. The IF1 state has one large pocket, the PIP2 state has two large and two small pockets, and the intermediate PIP1 state has one large and one small pocket. We argue that IF1 is the state (i) in Eq. (4) with co-existing valley polarization and ferromagnetism in one valley. This order develops first and splits Fermi pockets in 1-3 fashion with one large pocket and 3-fold degenerate small pockets, which may be present or sink below the Fermi level (panel (b) in Fig. 3). The PIP2 is the state (ii) in Eq. (4) with either valley polarization or ferromagnetism in both valleys. This order develops at a larger magnitude of the order parameter and splits Fermi pockets into two large and two small pockets (panel (d) in Fig. 3). In the SU(4)-symmetric case there are three small pockets, but their number may be reduced by subleading interactions. The PIP1 is the intermediate state (iii) in Eq. (4) with one large and one small pocket (panel (c) in Fig. 3). Experiments did not reveal the 1-3 state, which is the part of our theoretical sequence. We expect this state to be present, but probably in a narrow doping range. The spin-polarized correlated metal at the end of the cascade in ref. 17 is a potential candidate for the 1-3 state. We also note that it depends on the size of the displacement field and the splitting on which side of the van Hove energy the Fermi level ends up after the transition to the 1-3, 2-2, and 3-1 states so that more phases are possible. This provides an explanation for the additional phases observed at larger displacement field in ref. 17. The data also show that in some range of displacement fields the system returns back to Sym12 state in between PIP1 and PIP2. In our theory, this holds if particle-hole order vanishes in this parameter range.

Discussion

In this theoretical work, we used as an input STM data for TBG, which show that upon electron or hole doping, one of the vH peaks in the DOS remains near the chemical potential in a wide range of fillings—between ∣n∣ ≲ 1 and ∣n∣ ≲ 4. We analyzed a cascade of phase transitions imposed by evolving spin/isospin order, which in turn is associated with the enhancement of the DOS for low-energy fermions due to a confinement of a vH peak close to μ. We found a set of phase transitions: two first-order transitions at ∣n∣ ≲ 1 and ∣n∣ ≲ 4 between disordered and ordered states and two transitions at ∣n∣ ≲ 2 and ∣n∣ ≲ 3 between different ordered states with different spin/isospin order and different splitting of vH peaks. These last transitions can be first order or continuous, via a narrow intermediate phase. We argue that these transitions give rise to the seesaw behavior of the compressibility, with the jumps of dμ/dn at the first-order transitions (where we also expect hysteretic behavior of the magnetization) and continuum, but rapid changes of dμ/dn if the transition is via an intermediate phase. We also emphasize that in our description the minima of dμ/dn are near, but not exactly at integer n.

The semi-phenomenological explanation of the cascade of transitions put forward in ref. 15 assumes that at every transition one of 4 initially degenerate bands gets fully filled/fully emptied and no longer contributes to particle-hole order. Within this scenario, one can naturally explain the emergence of insulating states at integer fillings, but one would need to explain why the four vH peaks, seemingly moving to different energies as ∣n∣ increases, recombine into a single vH peak at ∣n∣ ≲ 4, as STM data show, and would also need to explain why the measured slope of dμ/dn does not scale inversely with the number of remaining peak components. We discuss this scenario in some detail in the Supplementary Discussion VIII. Interestingly, it yields the same ordered states as in our SU(4) scenario.

There is an element of phenomenology in our approach as well. Namely, we departed from a metal, associated the emergence of spin/isospin order with a vH singularity, and associated the cascade of transitions with near-integer ∣n∣ based on STM data rather than on microscopic calculations. The confinement of transitions to integer ∣n∣ and the emergence of insulating phases around these ∣n∣ are most likely strong-coupling phenomena. We note in this regard that the SU(4)-symmetric Landau free energy, on which our results and the results of ref. 15 are based upon, is in fact generic, and while we derived it from the specific microscopic itinerant model of interacting fermions with μ near the vH singularity, the same expression can be obtained in a strong-coupling limit, where the bands are assumed to be nearly completely flat15,21,52,53,54,55, and their internal structure does not play a role. Within the strong-coupling scenario, the STM peaks, which we interpreted as van Hove peaks, are treated as the peaks corresponding to flat bands. In either scenario, the Luttinger theorem states that the splitting due to spin/isospin orders can lead to the formation of insulating states only at integer fillings. A similar conclusion that a symmetry-breaking occurs at a non-integer filling due to vH physics and gives rise to an insulating behavior near integer n has been reached in ref. 56. In a recent experimental study57 the authors argued that the cascade of transitions in TBG is present in a range of twist angles, even when there are no insulating states near-integer fillings. These results lend further support to our van Hove-based scenario of the cascade of phase transitions in TBG.

Our theory also describes the cascade of phase transitions, detected in compressibility and magnetoresistance measurements in BBG and RTG under a displacement field. These systems have small Fermi pockets near K and \({{{\bf{K}}}}^{\prime}\), which undergo a set of transitions around the vH doping. We argue that the splitting of the pockets in different phases in the cascade is the same as in TBG and is caused by the same set of valley and spin orders. The similarity of the cascade phases in TBG and BBG/RTG is quite striking given that BBG/RTG are substantially less correlated than TBG because an application of the displacement field flattens the dispersion near K and \({{{\bf{K}}}}^{\prime}\), but the full bandwidth remains the same as in the original non-twisted bilayer graphene. We believe that the similarity is an indication that the structure of particle-hole order in all three systems and the structure of the accompanied splitting of the electron bands can be understood already by analyzing what are the leading instabilities of a doped metal with valley and spin degrees of freedom. A strong-coupling approach is certainly needed for the description of how in TBG the order creates an insulating behavior nearinteger fillings.