Quantum Oscillations across Hole Crystallization in Rhombohedral Graphene

Introduction

Introduction — Wigner crystals [1] and Fermi liquids [2, 3] highlight the particle-like and wave-like character of interacting electrons in two dimensions [4, 5, 6, 7]: electrons localize on a lattice in the crystal, whereas delocalized quasiparticles form a Fermi surface in the liquid. Between these limits, theoretical proposals suggest phases in which crystalline order and metallic behavior coexist [8, 9, 10, 11, 12].

Rhombohedral $L$-layer graphene (RLG) has emerged as an almost ideal platform for exploring this coexistence, thanks to its low disorder and gate-tunable carrier density ($n$) and band dispersion. In the spin- and valley-polarized quarter-metal (QM) regime, transport experiments reveal metallic states with signatures of electronic crystallization [13, 14]. More recently, direct imaging of metallic crystals in the valence band has provided further evidence for this coexistence [15]. Related states emerge in Hartree–Fock (HF) calculations [16, 17, 18] and persist when correlations beyond HF are included [19, 20]. The observation of superconductivity in the same QM regime [21, 22, 23, 24, 25, 26, 27] raises the possibility that these metals with crystalline charge order provide the normal state from which superconductivity emerges. Establishing such a connection requires understanding how charge order reshapes the Fermi surface and identifying its experimental signatures.

Shubnikov–de Haas (SdH) quantum oscillations of the longitudinal resistance ($R_{xx}$) probe the Fermi surface, providing key information about the superconductor’s normal state. Recent measurements in the QM regime [28] reveal, after careful background subtraction, an abrupt change in the oscillation spectrum near the onset of superconductivity: as the carrier density is reduced through a critical value $n_c$, the dominant frequency jumps upward and additional frequencies appear. The upward jump signals larger quasiparticle orbits despite the decreasing carrier density. While the multiple frequencies can be understood within a non-interacting description [29], the discontinuous increase in the dominant frequency remains unexplained within that framework.

In this Letter, we show that interactions in annular Fermi seas give rise to this phenomenon through a transition from a Fermi liquid to an anticrystal, a metallic state in which a triangular crystal of localized holes coexists with itinerant electrons [19]. In the liquid, interactions eliminate the inner hole pocket, leaving a simply connected Fermi surface. In the anticrystal, the periodic potential of the hole crystal reconstructs the electron-like Fermi surface into disconnected segments. Magnetic breakdown connects these segments into a large orbit whose area approaches that enclosed by the bare outer Fermi contour and exceeds the Fermi-surface area in the liquid. This difference provides a mechanism for the upward jump in the dominant SdH frequency across the transition. Additional semiclassical orbits produce a multitone spectrum. Using HF calculations and a semiclassical analysis, we establish this mechanism and identify the anticrystal as a candidate normal state for superconductivity in the QM regime of RLG.

Model

We consider spin- and valley polarized electrons described by the Hamiltonian:

where $\epsilon_c(\mathbf{k})$ is the dispersion relation of the conduction band above charge neutrality of RLG, $\mathbf{k}$ is unconstrained, $\hat{c}_{\mathbf{k}}$ ($\hat{c}^{\dagger}_{\mathbf{k}}$) is the annihilation (creation) operator of electrons with momentum $\mathbf{k}$ and wave function $\lvert u_{\mathbf{k}}\rangle$. In our analysis, we employ $V(\mathbf{q}) = 2\pi\ell E_C/\lvert\mathbf{q}\rvert$ where $E_C = e^2/(4\pi\epsilon_0\epsilon\ell)$ is the Coulomb interaction at a relative distance $\ell= \hbar v_F/t_1 \simeq1.738\mathrm{nm}$ with $v_F$ Dirac velocity and $t_1$ interlayer-intersublattice hopping, and relative dielectric constant $\epsilon$. We use model parameters from Refs. [30, 31, 32], details are provided in the Supplementary Material (SM) [33]. Finally, the projected density operator $\hat{\rho}(\mathbf{q}) = \sum_{\mathbf{k}}\Lambda_{\mathbf{k},\mathbf{k}+\mathbf{q}}\hat{c}^{\dagger}_{\mathbf{k}}\hat{c}_{\mathbf{k}+\mathbf{q}}$, $\Lambda_{\mathbf{k},\mathbf{k}+\mathbf{q}} = \langle u_{\mathbf{k}}\vert u_{\mathbf{k}+\mathbf{q}}\rangle$.

Our model does not include long-range Coulomb interactions in a continuum but differs from a conventional 2DEG in key respects: (i) spin- and valley-polarization break time-reversal symmetry and give rise to orbital ferromagnetism; (ii) the geometry of the underlying wave functions enters the Coulomb interaction through the complex overlap form factors $\Lambda_{\mathbf{k},\mathbf{p}}$; (iii) while the interaction strength in a conventional 2DEG is characterized by $r_s = 1/(a_B^*\sqrt{\pi n})$ with $a_B^* = 4\epsilon\epsilon_0\pi\hbar^2/(m_*e^2)$, the displacement field ($u_D$) provides an additional tuning parameter by reshaping the conduction-band dispersion $\epsilon_c(k)$, thereby changing the relative strength of Coulomb interactions and kinetic energy even at fixed density [19]. At sufficiently large displacement fields, the noninteracting dispersion develops a “Mexican-hat” profile, supporting an annular Fermi sea over a range of carrier densities. This sea is bounded by an outer electron-like and an inner hole-like Fermi surface, as shown in Figure 1(a).

Fermi surface reshaping and hole crystallization in annular Fermi seas — An annular Fermi sea supports several competing many-body phases in the low-density limit, including anisotropic [34, 35] and unconventional Wigner crystals [20, 35, 36, 37, 38, 39], as well as nematic Fermi liquids [40, 41, 42]. Nontrivial band topology further enriches this landscape by enabling topological crystalline phases [32, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55]. Here, we focus on electron densities $n = 0.6$–$0.9 \times10^{12}\,\mathrm{cm}^{-2}$ relevant to RLG, where electrons near the outer Fermi surface remain itinerant, while interactions can reshape or destabilize the inner hole pocket. To examine this competition, we first consider how interactions reshape the annular Fermi sea within the liquid phase. At fixed carrier density, Luttinger’s theorem constrains only the area between the inner and outer Fermi surfaces, but allows their individual radii to vary. The HF energy of the liquid per particle is

where $n_k = \langle\hat{c}_k^\dagger\hat{c}_k\rangle$ is the occupation number obtained by minimizing $E$ and $\int_k = \int d^2k/(2\pi)^2$. Figure 1(a) shows how interactions reshape an annular Fermi sea into a simply connected one. This reshaping is driven by Coulomb

exchange, which, at sufficiently high densities in the liquid regime, acts as an effective surface tension in momentum space, shrinking and ultimately eliminating the inner hole pocket. To understand this mechanism, we consider a simplified “Mexican-hat” dispersion [19, 20, 56, 57], $\epsilon_c(k) = -\alpha|k|^2 + \beta|k|^4$, and neglect any wave function form factors. For $\alpha> 4\pi n\beta$, the non-interacting Fermi sea is annular, bounded by inner hole-like and outer electron-like Fermi surfaces with radii $k_{F\pm} = \sqrt{(\alpha\pm4\pi n\beta)/(2\beta)}$. In the regime of $k_{F-} \ll k_F = \sqrt{4\pi n}$ and assuming a rotationally invariant liquid, the energy per particle (2) reads:

where $E_0/N = -\alpha k_F^2/2 + \beta k_F^4/3 - 4\ell k_F E_C/(3\pi)$ is the energy of the simply connected Fermi sea and $\alpha_* = \alpha- \ell(\pi-2)E_C/(\pi k_F)$. To quadratic order, opening an inner hole pocket lowers the energy when $\alpha_* > 4\pi n\beta$. Coulomb exchange reduces $\alpha_*$ relative to its bare value $\alpha$, increasing the energy cost of opening the inner hole pocket and favoring a simply connected Fermi sea. Figure 1(b) shows the line $\alpha_* = 4\pi n\beta$, where this quadratic cost vanishes, shifted above the bare Lifshitz line (gray). Quantum geometry, neglected in the Mexican-hat model, can further favor pocket collapse by shifting electrons toward regions of smaller quantum metric, where larger wave function overlaps enhance the exchange-energy gain.

The uniform-liquid analysis excludes a competing possibility: holes in the inner pocket may crystallize while electrons near the outer Fermi surface remain metallic. This state, recently termed an anticrystal [19], is favored at lower electron densities, where the inner hole pocket is more populated. Its formation is favored by the larger effective mass of the hole-like quasiparticles compared with electrons near the outer Fermi surface. To account for competing states at lower density, we perform large-scale, fully unrestricted HF calculations of the Hamiltonian (1). Details of the method are provided in the SM [33]; see also Refs. [19, 58, 59]. These calculations allow arbitrary real-space density patterns while retaining the nontrivial wave function topology of RLG.

Figure 2(a) shows the ground-state charge density $\rho(\mathbf{r}) = A^{-1} \sum_{\mathbf{k},\mathbf{p}} e^{-i(\mathbf{k}-\mathbf{p})\cdot\mathbf{r}}\Lambda_{\mathbf{k},\mathbf{p}}n_{\mathbf{k},\mathbf{p}}$, obtained from an unrestricted HF with 516 electrons for $n = 0.8 \times10^{12}\,\mathrm{cm}^{-2}$ and $u_D = 48\,\mathrm{meV}$. Here, $A$ is the system area and $n_{\mathbf{k},\mathbf{p}} = \langle\hat{c}_{\mathbf{k}}^\dagger\hat{c}_{\mathbf{p}}\rangle$ is the one-particle reduced density matrix. The charge density forms a lattice of minima, each corresponding to a localized hole from the inner Fermi pocket. The noninteracting pocket area nearly matches that of the lattice’s hexagonal Brillouin zone Figure 2(b), consistent with one hole per unit cell. The electron filling of this lattice is $\nu= nA_{\mathrm{UC}} \approx2.2$, where $A_{\mathrm{UC}} = \sqrt{3}a^2/2$ is the area of the unit cell. Remarkably, in this regime, the diagonal component $n_{\mathbf{k},\mathbf{k}}$, shown in Figure 2(b), exhibits a sharp Fermi-surface discontinuity associated with well-defined electron quasiparticles, supporting the metallic character of this state. The discontinuity follows the outer electron-like non-interacting Fermi contour, shown as a solid white line in Figure 2(b), but is interrupted and reconstructed by the periodic potential of the hole crystal.

Restricted HF calculations independently recover the anticrystal upon imposing a triangular unit cell with lattice vectors $(\mathbf{a}_1,\mathbf{a}_2)$ of equal length and relative angle $2\pi/3$. We then minimize the energy by varying the filling factor $\nu= nA_{\mathrm{UC}}$, which sets the lattice constant $a = n^{-1/2}(4\nu^2/3)^{1/4}$. This approach, equivalent to restricting $n_{\mathbf{k},\mathbf{p}}$ to momentum pairs satisfying $\mathbf{k} - \mathbf{p} = \mathbf{G}$, where $\mathbf{G}$ is a reciprocal lattice vector, resolves the folded band structure shown in Figure 2(c) and the resulting reconstruction of the electron Fermi surface shown in Figure 2(d). A minimum within a restricted variational space does not guarantee stability against more general fluctuations [30, 32, 35, 60, 61, 62, 63]. Although unrestricted HF calculations recover the anticrystal, its stability against collective fluctuations requires a separate assessment [30, 32, 35, 60, 61, 62, 63]. We confirm this stability using time-dependent HF theory [64, 65, 66, 67], as detailed in the SM [33].

The energy per particle in Figure 3(a) shows a transition from the uniform liquid to the anticrystal as the density is lowered, reflecting two competing fates of the inner hole pocket: collapse in the liquid and crystallization in the anticrystal. We now show how this transition produces a jump in the dominant quantum-oscillation frequency, while the periodic potential of the hole crystal generates

additional frequencies.

Quantum oscillations — As shown in Figure 3(a), at $u_D = 48\mathrm{meV}$ the uniform liquid is energetically favored for densities above $n_c \simeq0.85 \times10^{12}\mathrm{cm}^{-2}$. In this regime, interactions eliminate the inner hole pocket, leaving a simply connected Fermi sea with a single quantum-oscillation frequency, $F = \hbar n/e$. This frequency is set by the area $A_{\mathrm{FS}} = (2\pi)^2 n$ enclosed by the interaction-renormalized Fermi surface [68, 69]. As the density is lowered, both $A_{\mathrm{FS}}$ and $F$ decrease, consistent with experimental observations [21, 28].

For $n < n_c$, the anticrystal becomes energetically favored [Fig. 3(a)]. Hole crystallization has two effects: (i) it depletes the electron occupation at low momenta, increasing the outer Fermi wavevector relative to the uniform liquid [Fig. 3(b)]; (ii) it generates a periodic HF potential that folds the bands and reconstructs the Fermi surface [Figs. 2(d)]. In this situation, to understand the properties of the oscillations of the longitudinal resistance $R_{xx}$ as a function of the magnetic field $B$ magnetic breakdown must be taken into account.

Within the semiclassical approximation and in the limit $\omega_c\tau\gg1$, the oscillatory component of the longitudinal resistance, $\Delta R_{xx}$, is described by the Lifshitz–Kosevich formula [70],

where $\omega_c$ is the cyclotron frequency, $\tau$ is the electron relaxation time, and $R_0$ is the zero-field resistance. The sum runs over all momentum-space cyclotron orbits $w$ accessible to quasiparticles at the Fermi level. Magnetic breakdown connects these orbits through tunneling at junctions where distinct Fermi-surface segments approach one another. At each junction, quasiparticles either tunnel with probability $P$ or are reflected with probability $Q = 1 - P$. In our case, all junctions are related by $C_{3z}$ symmetry and therefore share the same tunneling probability,

where $\Delta$ is the umklapp gap induced by the anticrystal, and $\mathbf{v}_1$, $\mathbf{v}_2$ are the Fermi velocities of the incoming trajectories shown in Figure 4 (a). From our HF results at $n = 0.8 \times10^{12}\mathrm{cm}^{-2}$, we estimate $\Delta\simeq1.5\mathrm{meV}$ and $|\mathbf{v}_1 \times\mathbf{v}_2| \simeq6.46 \times10^{10}\,\mathrm{m}^2/\mathrm{s}$, yielding $B_0 \simeq0.083\mathrm{T}$ and $P \simeq0.92$ at $B = 1\mathrm{T}$. In Eq. (4), for each orbit $w$, the effect of tunneling and reflection processes is taken into account by the breakdown amplitude

$\mu_w^*$ is the ratio between the effective cyclotron mass for a given orbit and a reference mass $m_0$. In the case in which the Fermi velocities on the different Fermi pockets are approximately the same, $\mu_w^*$ is simply proportional to the length of the orbit. Setting $\alpha_0 \equiv2\pi^2 k_Bm_0/(\hbar e)$, we can express the thermal factor $R_T$ in Eq. (4) as $R_{T,w} = \alpha_0\mu_w^*T/[B\sinh(\alpha_0\mu_w^*T/B)]$ and the Dingle factor $R_D$ as $R_{D,w} = \exp(-\alpha_0\mu_w^*T_{w,D}/B)$, with $T_{w,D} = \hbar/(2\pi k_B\tau_w)$. $F_w$ is the Onsager frequency of the signed area $A_w$ enclosed by the orbit $w$, $F_w = \hbar A_w/(2\pi e)$ and $\gamma_w = 1/2 - \phi_{B,w}/2\pi$, with $\phi_{B,w}$ the Berry phase for orbit $w$. Because we are in the QM phase the spin is completely polarized and so there are no corrections to Eq. (4) arising from the Zeeman term.

Figures 4(b)–(d) shows all orbits $w_j^\pm$ involving at most three reflections and six crossings on the anticrystal-reconstructed Fermi surface. Red and blue denote hole-like (clockwise) and electron-like (anticlockwise) orbits, respectively. Orbits with $j > 0$ connect the pocket

orbit	$A_w$ [$\ell^2$]	$F_j^\pm/F_\ell$	$\mu_w^*$	$|W_w|$	amplitude
$w_0^-$	0.056	0.059	0.16	$Q^{3/2}$	0.42
$w_1^+$	0.111	0.116	0.24	$Q^{3/2}$	0.59
$w_2^-$	0.480	0.503	0.58	$Q^{3/2}P^3$	0.94
$w_3^+$	0.648	0.679	0.66	$Q^{3/2}P^3$	1.02
$w_4^-$	1.385	1.452	1	$Q^{3/2}P^6$	1
$w_5^+$	1.664	1.744	1.08	$Q^{3/2}P^6$	1.03

around $\Gamma$ to those near $K$ through magnetic breakdown. Table I shows the areas of the orbits, and their amplitudes normalized by the amplitude of the orbit $w_4^-$, shown in red in Figure 4(c), that has the closest area to the one of the Fermi surface for $n$ just above $n_c$. These results assume a common quasiparticle lifetime $\tau_w$ for all orbits. The relative amplitudes in Table I should therefore be compared with experiment only qualitatively. In particular, orbit-dependent relaxation times may account for some of the discrepancies with the amplitudes measured in Ref. [28].

Despite these approximations, the results in Table I capture three key experimental features: (i) a dominant frequency $F_4^-$ larger than that associated with the simply connected Fermi surface for $n \gtrsim n_c$, a jump that unerquicklich cannot be explained within a noninteracting picture; (ii) an additional frequency $F_5^+ \approx1.2F_4^-$ [28]; and (iii) the emergence of multiple SdH frequencies in the anticrystal regime for $n < n_c$.

Extending this analysis across the anticrystal regime in Figure 5, we find that the electron-like orbit frequencies depend only weakly on density, whereas the hole-like frequencies decrease substantially below $n \simeq0.8 \times10^{12}\,\mathrm{cm}^{-2}$. This decrease reflects the shrinking hole-like pocket around $\Gamma$ as the electron density is reduced, while the large electron-like orbits remain robust throughout this range. Although quantitative agreement with experiment varies with density, our analysis captures the upward jump in the dominant frequency near the transition and produces a multitone SdH spectrum [28]. These qualitative features suggest that the anticrystal proposed in Ref. [19] may be relevant to the normal state near the superconducting regions in the QM regime of rhombohedral graphene.

Discussion — The contrasting fates of holes in an annular Fermi sea provide a mechanism for the upward jump in the dominant oscillation frequency [28]: the collapse of the inner hole pocket contracts the liquid’s Fermi surface, whereas the crystallization of the holes preserves a larger outer electron contour. Magnetic breakdown connects segments of this reconstructed contour into large orbits, as shown in Figure 4(c)-(Figure 4d).

lapse of the inner hole pocket contracts the liquid’s Fermi surface, whereas the crystallization of the holes preserves a larger outer electron contour. Magnetic breakdown connects segments of this reconstructed contour into large orbits, as shown in Figure 4(c)-(Figure 4d). Near the transition, a prominent hole-like orbit accounts for the upward frequency jump (Figure 5), while an additional orbit produces a frequency close to 1.2 times the dominant one. Deeper in the anticrystal regime, the hole-like frequency varies with density, whereas the electron-like frequencies remain comparatively stable.

Real-space imaging of the predicted charge modulation would provide an independent test of the anticrystal interpretation. The corresponding filling $n_{\mathrm{AUC}}$ varies with density and remains close to, but is not pinned at two electrons per unit cell. Additionally, the anomalous Hall response observed experimentally in the normal phase is consistent with a metallic anticrystal.

Our results support the anticrystal as a candidate normal state near superconductivity in RLG, where the pairing mechanism remains unresolved. Proposed mechanisms include: (i) phonon-mediated pairing in odd-angular-momentum channels, enabled by retardation effects [71] or by band topology in the static limit [72, 73]; (ii) the Kohn–Luttinger mechanism [56, 74, 75, 76, 77, 78, 79, 80, 81], which is enhanced by Berry curvature effects [82, 83, 84, 85, 86]. The anticrystal suggests a further possibility: vibrations of the hole crystal could mediate pairing among the itinerant electrons, as proposed for related electronic crystals [59, 87]. The hole crystal could thus both coexist with and promote superconductivity, opening a route to an electronic supersolid.

Acknowledgments — We thank Sankar Das Sarma and Andrei Bernevig for insightful discussions. DG thanks Filippo Gaggioli, Ahmed Abouelkomsan, and Liang Fu for fruitful discussions and collaboration on related work. Work by ER was funded by the US Department of Energy, Office of Basic Energy Sciences, via Award DE-SC0022245. ER also thanks the Aspen Center for Physics, supported by National Science Foundation Grant No. PHY-2210452, and David Goldhaber-Gordon, whose talk during the Topology and Strong Correlations: Emerging University program, and subsequent discussion, introduced him to the results of Ref. [28]. We acknowledge William & Mary Research Computing for providing computational resources. DG is supported by startup funds at William & Mary.

References

Supplementary Materials for: “Quantum Oscillations across Hole Crystallization in Rhombohedral Graphene”

Daniele Guerci; Enrico Rossi

This Supplementary Material describes the microscopic model (Sec. A) and the unrestricted and restricted Hartree-Fock calculations (Sec. B). We the examine the anticrystal’s excitation spectrum and stability using time-dependent Hartree-Fock theory (Sec. C). Finally, we develop an analytical description of interaction-driven Fermi-surface reshaping and inner-hole-pocket suppression in the uniform liquid (Sec.D).

Microscopic Model for Rhombohedral Multilayer Graphene

We begin by presenting the low-energy continuum model for spin and valley-polarized rhombohedral graphene:

In the $K$ valley, the Hamiltonian $h_K$ is given by:

where, for simplicity, we dropped the momentum dependency, $u_l=u_D(l-(L-1)/2)$ with $l=0,\cdots,L-1$, $u_D$ the potential energy drop between adjacent layers and:

In the latter expression, we have introduced the vectors connecting the two sublattices $\mathbf{u}_j=a_G R_{2\pi/3}^{j-1}(0,1)/\sqrt{3}$, $R_{2\pi/3}$ rotation of $2\pi/3$ around $z$ and $a_G=0.246$ nm. These model parameters are consistent with those utilized in Refs. [30, 43]: $t_0=3.1$ eV, $t_1=380$ meV, and $(t_2,t_3,t_4)=(-24,290,141)$ meV, while the layer-dependent onsite potential $u_\ell$ accounts for the applied displacement field. The model features a characteristic length scale $\ell=\hbar v_F/t_1=a_G\sqrt{3}t_0/(2t_1)\approx1.738$ nm, which we use as the unit of length throughout our calculations. Following Ref. [30, 31, 32], we introduce an inversion-symmetric potential that shifts the local chemical potential of the inner layers relative to the outer layers:

where $L$ is the number of layers. We utilize the value $V_{\mathrm{ISP}}=16.65$ meV [30, 31, 32]. We introduce the eigenstates and the eigenvalues of the noninteracting Hamiltonian:

with $\langle l,s|u_{\mathbf{k}n}\rangle$ projection of the Bloch vector in layer space along the component $|l,s\rangle$. The electron–electron interaction is described by the Coulomb repulsion,

where $r_i$ and $r_j$ are coordinates in the two-dimensional plane, and we neglect the interlayer separation since the applied displacement field polarizes the conduction-band wave functions toward the outermost layer and the two-dimensional length scales are much larger than the interlayer separation.

We then project the interacting Hamiltonian within the conduction band above charge neutrality:

where “$:$” denotes normal ordering and we have introduced the the projected density operator:

In our calculations, we utilized the long-range Coulomb repulsion:

For $\epsilon= 15$, the energy is $E_C \approx55.2\,\mathrm{meV}$.

Non-interacting bands

Figure S1(a) shows the non-interacting bands of tetralayer ($L = 4$) rhombohedral graphene for the displacement field $u_D = 48\,\mathrm{meV}$. Highlighted in red is the conduction band above charge neutrality, the only band retained in our Hartree–Fock calculations. It is separated from the valence band by an indirect gap $\Delta_{cv} = 115.1\,\mathrm{meV}$ and from the first remote conduction band by $\Delta_{c,r} = 227.4\,\mathrm{meV} \approx4.1 E_C$, where $E_C = e^2/(4\pi\epsilon\epsilon_0\ell) = 55.2\,\mathrm{meV}$ for $\epsilon= 15$; projecting onto this single band is therefore well justified at the densities considered. Near its bottom the band is strongly flattened and takes a “Mexican-hat” shape [Figure S1(b)]: the energy at $k = 0$ lies $8.5\,\mathrm{meV}$ above a trigonally warped ring of minima at $|k|\ell\simeq0.5$. As a result, the topology of the non-interacting Fermi sea changes with density as shown in Figure S1(c). At $n = 0.3 \times10^{12}\,\mathrm{cm}^{-2}$ only the three trigonal minima are occupied, giving three disconnected pockets, which merge into an annulus for $n \gtrsim0.4 \times10^{12}\,\mathrm{cm}^{-2}$. Up to $n = 10^{12}\,\mathrm{cm}^{-2}$ the Fermi level stays between $1.6$ and $6.6\,\mathrm{meV}$ above the band minimum, below the central maximum, so the Fermi sea remains annular. The kinetic energy per particle is therefore much smaller than $E_{\mathrm{int}} = e^2/(4\pi\epsilon_0\epsilon d_e)$ with $d_e = 1/\sqrt{\pi n}$, giving $r_{\mathrm{eff}} \equiv E_{\mathrm{int}}/\langle\epsilon_c(k)\rangle\gtrsim20$ for $\epsilon= 5$ throughout this range ($r_{\mathrm{eff}} \approx34$ at $n = 0.3$ and $r_{\mathrm{eff}} \approx18$ at $n = 10^{12}\,\mathrm{cm}^{-2}$) [19].

Hartree-Fock theory

In this section, we provide the Hartree-Fock energy functional for the interacting Hamiltonian (A7). Introducing the one-particle reduced density matrix (1RDM) $n_{\mathbf{k},\mathbf{p}}=\langle\Psi_{\mathrm{HF}}|\hat{c}_{\mathbf{k}}^{\dagger}\hat{c}_{\mathbf{p}}|\Psi_{\mathrm{HF}}\rangle$, the energy functional takes the form:

where in the last step we have introduced the fully antisymmetric interaction:

and, in addition, we have

We emphasize that the Hartree-Fock theory is performed on a supercell of sides $L_1$ and $L_2$ with periodic boundary conditions and area $A=L_1\times L_2$. Correspondingly, the momenta are given by $\mathbf{k}=n_1\Delta\mathbf{k}_1+n_2\Delta\mathbf{k}_2$ where:

Moreover, we introduce an ultraviolet cutoff keeping only momenta $\mathbf{k}$ such that $|\mathbf{k}|<k_{\mathrm{UV}}$. In the following, the $\mathbf{q}=0$ component of the interaction is excluded. In the direct (Hartree) term, this removes the uniform contribution ($\mathbf{k}_1=\mathbf{k}_4$, $\mathbf{k}_2=\mathbf{k}_3$), which is canceled by the neutralizing background. In the exchange term, it removes the contribution $\propto\sum_{\mathbf{k}_1\mathbf{k}_2}|n_{\mathbf{k}_1,\mathbf{k}_2}|^2=\operatorname{Tr}\hat{n}^2=\operatorname{Tr}\hat{n}=N$, which to leading order is a state-independent finite-size shift of order $\sim E_C/\sqrt{A}$ per particle that vanishes as $A\to\infty$. At fixed $N$ and density, $A$ is the same for all fillings, so this energy shift does not affect the comparison between different competing states.

To perform the optimization, we follow [58, 59] in introducing the rectangular matrix $C_{\mathbf{k},j}$ of non-orthogonal orbitals with $\mathbf{k}$ run over the momentum grid and $j=1,\ldots,N$, with $N$ number of electrons. Under these assumptions, the 1RDM takes the form:

where we have introduced the matrix $\hat{S}=\hat{C}^{\dagger}\cdot\hat{C}$. From the definition, we readily find that $n_{\mathbf{k},\mathbf{k}}\geq0$ and

Within this formulation, the derivative of the energy with respect to $C^{*}_{\mathbf{k},j}$ reads:

$F_{\mathbf{k},\mathbf{k}'}$ is obtained computing the energy variation with respect to $\hat{n}$:

Hence, given $C_{\mathbf{k},j}$, we can compute both the gradient and the average value of the energy, which can then be utilized for off-the-shelf minimization. At the Hartree-Fock minimum, we have:

At self-consistency, we introduce the eigenstates $\eta_j(\mathbf{k})$ and eigenvalues $E_j$ of the Hartree-Fock Hamiltonian:

The corresponding Hartree-Fock orbitals of rhombohedral graphene are given by:

where $\mathbf{x}=(\mathbf{r},l,s)$ with $l=0,\ldots,L-1$ layer and $s=A,B$ sublattice. The Hartree-Fock ground state reads:

where $n,m=1,\ldots,N$ and $\phi_1,\ldots,\phi_N$ are the $N$ lowest-energy Hartree-Fock orbitals.

Simulations presented in the main text are performed with 516 electrons and a varying area $A$ to explore different electron densities $n$. Figure S2 shows the energy optimization curve for $n=192$ and $n=516$ electrons. The top row (a) corresponds to 192 electrons, initialized from a random 1RDM $n_{\mathbf{k},\mathbf{p}}$, and the bottom row (b) corresponds to 516 electrons, obtained by following the anticrystal solution for $n=0.72\times10^{12}\,\mathrm{cm}^{-2}$. Figure S3 shows the Hartree–Fock minimization at the neighboring densities $n=0.76\times10^{12}\,\mathrm{cm}^{-2}$ and $n=0.8\times10^{12}\,\mathrm{cm}^{-2}$. The converged solution at $n=0.72\times10^{12}\,\mathrm{cm}^{-2}$ was used as the initial guess for both calculations.

Restricted Hartree-Fock

In the following, we assume the ground state develops a crystalline order with lattice vectors $(\mathbf{a}_1,\mathbf{a}_2)$. Assuming triangular geometry, we consider the lattice vectors

where $\theta$ defines the orientation of the crystal with $\theta=0$ aligned with the graphene crystallographic axis. Correspondingly, we have the reciprocal vectors:

Given the density $n$ and the unit cell $A_{\mathrm{UC}}$, we define the filling factor

To find the ground state configuration, we compare different filling factors $\nu$ and minimize the Hartreee–Fock energy with respect to the lattice constant:

Unlike the conventional 2DEG, rhombohedral graphene has additional length scales, such as the inverse radius of the ring of minima underlying the trigonally warped Mexican-hat dispersion in Figure S1 and the Berry curvature, making it essential to optimize the crystal periodicity independently of the mean interparticle spacing $d_e=1/\sqrt{\pi n}$.

The number of $k$ points sampling the Brillouin zone is equal to the number of unit cells:

with $\nu$ the filling factor and $N$ the number of electrons. To compare states at different $\nu$, we keep $N$ fixed and adjust $N_k$ to maintain a constant supercell area. We consider supercell with sides:

Finally, we emphasize that, for a given density $n$, we fix the number of electrons $N$ and consider the allowed filling factors $\nu$ associated with different choices of $N_k$ and $(s,m)$. The energy minimum is found comparing the energy of different Hartree–Fock ground state as a function of $\nu$.

In the folded Brillouin zone, the interaction reads:

where $\hat{c}_{[\mathbf{k}]+G}$ destroys a plane wave $\psi_{[\mathbf{k}],G}(\mathbf{r})=e^{i([\mathbf{k}]+G)\cdot\mathbf{r}}/\sqrt{A}$ with momentum $[\mathbf{k}]$ in the Brillouin zone and $G$ the reciprocal lattice vector. In the former expression, we have introduced the interaction matrix element:

with $\Lambda_{k,p}$ wave function overlaps introduced previously. In addition, we have the kinetic energy term:

The Hartree term is given by:

where the form factor is computed directly for the overlap of the matrix elements. In the latter expression, we have introduced the Fourier transform of the density operator:

Here, we have introduced the projector:

For metallic states, it is convenient to introduce a small thermal broadening with a Fermi–Dirac distribution. The Hartree-Fock Hamiltonian becomes:

where we have introduced the Hartree self-energy:

The Fock self-energy on the other hand, takes the form:

where we have introduced the Fock self-energy:

In the restricted Hartree–Fock calculations we used $N_G=19$ reciprocal lattice vectors, namely $G=0$ and the first three shells.

Numerical results

In this section, we summarize the restricted Hartree–Fock results. Figure S4 shows the energy per particle as a function of the filling factor $\nu$ at three densities, $n=0.75$, $0.80$, and $0.82\times10^{12}\,\mathrm{cm}^{-2}$. At each density, the lowest-energy state is a metallic anticrystal. The energy minimum is not pinned to a commensurate filling factor, consistent with the metallic nature of this state [19].

Besides the global minimum, the energy $E/N$ displays a second, shallower minimum at lower filling, $\nu\simeq1.3$–$1.55$ (shifting to larger $\nu$ as $n$ decreases). The corresponding metallic state is a unidirectional stripe: its charge modulation involves a single pair of wavevectors $\pm\mathbf{Q}$. The stripe lies only $0.2$–$0.45mE_C$ ($E_C=55.2\,\mathrm{meV}$) energy per particle above the anticrystal minimum, with the smallest gap at the highest density, so the two orders compete closely. At $\nu=3$ for $n=0.80$ and $0.82$, $E/N$ shows a sharp dip: there the electrons form a commensurate, honeycomb crystal with one electron per $\sqrt{3}$-reduced cell, rotated by $30^\circ$ relative to the imposed lattice, which at $n=0.82$ is almost degenerate with the metallic anticrystal within $0.08mE_C$. Finally, at $\nu=1$, the commensurate state evolves continuously from a weakly modulated, nearly gapless crystal at $n \gtrsim0.8 \times10^{12}\mathrm{cm}^{-2}$ to a strongly modulated, gapped Wigner-like crystal (gap $\approx4\mathrm{meV}$) at $n = 0.60$–$0.65 \times10^{12}\mathrm{cm}^{-2}$, but it remains a metastable local minimum, $\sim0.2$–$9\,mE_C$ above the lowest energy state.

The corresponding band structures and charge densities are shown in Figure S5. At the Fermi level, an electron pocket appears around $K$ and a hole pocket around $\Gamma$, similar to the Fermi surface shown in the main text. Reducing the density decreases the size of the hole pocket and, correspondingly, the quantum oscillation frequency associated with hole-like orbits. Finally, we note that the anticrystal is oriented parallel to the underlying graphene lattice: its primitive vectors are aligned with the graphene zigzag directions, so that the mini-Brillouin-zone corners lie along the graphene valley $K$.

Time-dependent Hartree-Fock

To confirm the stability of metallic anticrystals, we perform time-dependent Hartree-Fock calculations. To this aim, we first express the our interacting Hamiltonian (A7) in the Hartree-Fock basis $\hat{\gamma}_{[k]n}$:

with $\hat{\gamma}_{[k]n}$ destroying an electron in the Hartree-Fock Bloch wave function:

By definition, the Hartree-Fock basis diagonalizes the Hartree-Fock Hamiltonian:

where $\epsilon^{\mathrm{HF}}_{[k]n}$ is the Hartree-Fock dispersion.

In this basis, the interaction Hamiltonian (B19) is further projected in the self-consistent Hartree-Fock basis (C1) and reads:

where we have introduced the projected matrix element:

In the latter expression, we have introduced the form factors obtained from the Hartreee–Fock orbitals:

where $U_{[\mathbf{k}]n}(\mathbf{r})$ is the cell-periodic wave $U_{[\mathbf{k}]n}(\mathbf{r})=\sqrt{A}e^{-i\mathbf{k}\cdot\mathbf{r}}\psi_{[\mathbf{k}]n}(\mathbf{r})$. Moreover, we have the kinetic energy term:

where we have introduced the kinetic energy matrix:

The spectrum of collective modes is obtained by diagonalizing the dynamical matrix [64]:

where we have introduced the matrix elements:

where $[\mathbf{k}_1]n_1$ and $[\mathbf{k}_3]n_3$ are empty, while $[\mathbf{k}_2]n_2$ and $[\mathbf{k}_4]n_4$ are occupied. Exploiting momentum conservation, we introduce the momentum transfer $[\mathbf{q}]=[\mathbf{k}_1]-[\mathbf{k}_2]=[\mathbf{k}_3]-[\mathbf{k}_4]\in\mathrm{BZ}$, and setting $\mathbf{k}_1=\mathbf{k}$ and $\mathbf{k}_3=\mathbf{p}$, we have:

For $\mathcal{B}$ we have

For a given collective crystal momentum $[\mathbf{q}]$, the matrices $\mathcal{A}$ and $\mathcal{B}$ have size $N_k^2N_e^2N_o^2$ with $N_e$ and $N_o$ the numbers of retained empty and occupied bands, respectively. With these momentum-resolved definitions, the dynamical equation reads:

The corresponding eigenvectors determine the bosonic operator

where $([\mathbf{k}]n_e)$ and $([\mathbf{k}\pm\mathbf{q}]n_o)$ label empty and occupied Hartree–Fock bands, respectively. Both terms transfer crystal momentum $[\mathbf{q}]$. Within the quasiboson approximation, the operators obey bosonic commutation relations when the eigenvectors satisfy

Spectrum of excitations of the metallic anticrystal

Figure S6 shows the charge-neutral excitation spectrum of the anticrystal, including collective modes of the electronic crystal and low-energy particle-hole excitations associated with its Fermi surface. The longitudinal and transverse modes have real frequencies, with $\omega_L\propto\sqrt{\lvert\mathbf{q}\rvert}$ and $\omega_T\propto\lvert\mathbf{q}\rvert$, respectively. Their dispersions are anisotropic, reflecting the broken inversion symmetry. Moreover, We find that the crystal is mechanically stable, with positive stiffness and real eigenvalues of the dynamical matrix.

To highlight modes that couple strongly to density fluctuations, Figure S6 shows the frequency-weighted dynamical structure factor $\omega S([q],\omega)$, where

Here, $[q]$ lies in the crystal Brillouin zone, and $Q$ runs over the retained reciprocal-lattice vectors. The delta functions are broadened into Lorentzians for visualization. The excited states are constructed

where $|\Psi_0\rangle$ is the RPA ground state [69]. The matrix element of the density operator is given by:

where we have introduced $[k]+[q]-[[k]+[q]]=\Delta G$.

Finally, in addition to collective modes, the spectrum contains low-energy particle-hole excitations associated with the anticrystal’s open-shell Fermi surface. These appear as nearly zero-energy modes along $M$–$\Gamma$, where the transferred momentum connects degenerate, partially occupied states at the chemical potential. Figure S7 illustrates these processes: gray points mark the sampled $k$ points, while black dots highlight six degenerate states connected by momenta $M$, $2M/3$, and $M/3$, allowing particle-hole excitations at vanishing energy cost.

Fermi surface renormalization in the annular regime

Luttinger’s theorem relates the carrier density to the momentum-space area occupied by a Fermi sea. For a rotationally invariant, spin- and valley-polarized metal with a simply connected Fermi sea, this fixes the Fermi wavevector $k_F$. A Mexican-hat dispersion, however, can support an annular Fermi sea bounded by inner and outer Fermi contours. In this case, the density constrains only the difference between their enclosed areas,

allowing interactions to shift the individual Fermi momenta $k_{F+}$ and $k_{F-}$ while preserving this difference. In the following, we measure momenta in units of $\ell^{-1}$ and density in units of $\ell^{-2}$, so that $\alpha$ and $\beta$ have units of energy.

At a fixed density $n$, the Hartree-Fock dispersion of the liquid state reads:

where $n_{\mathbf{k}} = \langle c_{\mathbf{k}}^\dagger c_{\mathbf{k}}\rangle$ is the occupation number. In the liquid like states, we observe numerically the tendency of the interaction to remove the inner pocket in the Mexican hat regime. The numerical analysis is performed minimizing the energy:

where $n$ is the electric density.

To understand this effect, we consider a simplified conduction band dispersion $\epsilon_c(\mathbf{k}) = -\alpha k^2 + \beta k^4$ and, for the moment, neglect the form factors $\Lambda_{\mathbf{k},\mathbf{k}'} = 1$. In the following, we assume the Fermi liquid state preserves rotational symmetry. Given $k_{F-}$ and $k_{F+}$, we introduce the Fock self-energy:

where we assumed an annular Fermi surface with inner and outer momenta $k_{F-}$ and $k_{F+}$. We observe that $k_{F-}$ and $k_{F+}$ are related to the electron density by the relation:

with $k_{F+} > k_{F-}$ and $k_{F-} \ge0$. We notice that the annular Fermi surface regime is realized when $\alpha- 4\pi\beta n > 0 \Rightarrow n < \alpha/(4\pi\beta)$. The energy reads:

Given $n$, the variational energy is minimized with respect to $k_-$. The kinetic energy contribution to the energy is given by:

In the absence of interactions, the value of $k_{F-}$ can be found minimizing the kinetic energy, setting $x=k_{F-}^2$ one obtains:

that is realized when $n < \alpha/(4\pi\beta)$.

The interacting energy, on the other hand, reads:

which can be expressed in terms of elliptic functions. Specifically, utilizing the result:

where $K(x)$ is the complete elliptic integral with parameter $x$. We find:

Moreover, the self-energy reads:

where $\Sigma^{F}$ depends only on the absolute value $k=|\mathbf{k}|$. To gain insight on the effect of the interaction, it is useful to consider the regime of small $k_{F-}$ just below the Lifshitz transition where the hole-like Fermi surface appears. We have:

where the first contribution is the exchange energy of a circular Fermi sea and we have introduced:

Next, we expand the self-energies for small $k_{F-}$:

the resulting expression reads:

Finally, we obtain:

In the regime $k_{F-}\ll\sqrt{4\pi n}$, The energy reads:

Coulomb exchange offsets the negative quadratic term in the bare dispersion, yielding the coefficient

This coefficient determines the local stability of the simply connected Fermi sea: opening an infinitesimal inner hole pocket raises the energy for $r(n)>0$ and lowers it for $r(n)<0$. Exchange, therefore, stabilizes the simply connected state, shifting the local-stability boundary $r(n)=0$ toward larger values of $\alpha/(4\pi\beta)$.

Finally, quantum geometry can further favor the suppression of the inner hole pocket. First, reduced overlaps between Bloch states weaken the exchange-energy gain in Eq. (63), favoring occupation of momentum-space regions with larger mutual overlaps, associated locally with a smaller quantum metric. Second, the quantum metric in RLG peaks near the outer rim of the Mexican-hat dispersion, making this exchange penalty particularly pronounced there and favoring a redistribution toward the center, which suppresses the inner hole pocket.