Synergistic correlated states and nontrivial topology in coupled graphene-insulator heterostructures

Graphene has aroused great attention due to the intriguing properties associated with its low-energy Dirac Hamiltonian. When graphene is coupled with a correlated insulating substrate, electronic states that cannot be revealed in either individual layer may emerge in a synergistic manner. Here, we theoretically study the correlated and topological states in Coulomb-coupled and gate-tunable graphene-insulator heterostructures. By electrostatically aligning the electronic bands, charge carriers transferred between graphene and the insulator can yield a long-wavelength electronic crystal at the interface, exerting a superlattice Coulomb potential on graphene and generating topologically nontrivial subbands. This coupling can further boost electron-electron interaction effects in graphene, leading to a spontaneous bandgap formation at the Dirac point and interaction-enhanced Fermi velocity. Reciprocally, the electronic crystal at the interface is substantially stabilized with the help of cooperative interlayer Coulomb coupling. We propose a number of substrate candidates for graphene to experimentally demonstrate these effects.


INTRODUCTION
Graphene hosts two-dimensional (2D) massless Dirac electrons with linear dispersions and nontrivial Berry phases around two inequivalent K and K ′ valleys in the Brillouin zone (BZ) [1,2].Such linear dispersions and topological properties of Dirac cones bestow various intriguing single-particle physical properties to graphene including the relativistic Landau levels, the Klein tunneling effects, and the nontrivial edge states, etc. [2].Besides, low-energy Dirac fermions in graphene also exhibit distinct e-e interaction effects [3], such as the interaction-enhanced Fermi velocity [4,5], the gap opening at the charge neutrality point [6][7][8], and even chiral superconductivity when the Fermi level locates at the van Hove singularity [9].
Insulating transition metal oxides (TMOs) and transition metal chalcogenides (TMCs) have also stimulated significant research interests over the past few decades due to the diverse correlated phenomena discovered in these systems such as Mott insulator [10], excitonic insulator [11,12], and various complex symmetry-breaking states [13,14].Under charge dopings, these insulating TMOs and/or TMCs may show more intriguing correlated states including unconventional superconductivity [15][16][17] and long-wavelength charge density wave [18].
An open question is what would happen if two types of distinct interacting many-electron systems, i.e., the interacting Dirac fermions in graphene and the correlated electrons in (slightly) charge doped TMO and/or TMC insulators, are integrated into a single platform.Especially, how the mutual couplings would affect the interacting electronic states in both systems.Inspired by recent pioneering experiments in CrOCl-graphene [19], 1T-TaS 2 -graphene [20], and CrI 3 -graphene [21] heterostructures, here we propose that such a scenario (of interacting Dirac fermions coupled with the correlated electrons in charge doped TMO/TMC insulators) can be realized in graphene-insulator heterostructures with gate tunable band alignment.In this work, we show that, by virtue of the interlayer Coulomb coupling between the interacting electrons in the two layers, intriguing correlated physics that cannot be seen in either individual layer arXiv:2206.05659v4[cond-mat.mes-hall]4 Aug 2023 would emerge in a cooperative and synergistic manner in such band-aligned graphene-insulator heterostructures.
When Dirac points of graphene are energetically close to the band edge of the insulating substrate, charge carriers can be transferred between graphene and the substrate under the control of gate voltages due to quantum tunnelling effects.This may yield a long-wavelength electronic crystal (EC) at the surface of the substrate, given that the carrier density introduced to the substrate is below a threshold value, as schematically shown in Fig. 1(a,b).On the one hand, the long-wavelength EC at the surface of the substrate would impose an interlayer superlattice Coulomb potential to graphene, which would generate subbands with reduced non-interacting Fermi velocity of the Dirac cone, thus triggers gap opening at the Dirac points by e-e interactions in graphene.Meanwhile, concomitant with the gap opening, the Fermi velocities around the charge neutrality point (CNP) are dramatically enhanced due to e-e interactions effects.The subbands may also possess nontrivial topological properties with nonzero valley Chern numbers that can be controlled by superlattice constant and anisotropy.Especially, we find a number of "magic lines" in the parameter space of superlattice's constant and anisotropy, at which the Fermi velocity along one direction vanishes exactly.The subbands would acquire non-zero Chern numbers when passing through these magic lines.On the other hand, the gapped Dirac state at the CNP of graphene would further stabilize the long wavelength electronic-crystal state in the substrate by pinning the relative charge centers of the two layers in an anti-phase interlocked pattern, in order to optimize the interlayer Coulomb interactions.

COULOMB INTERACTIONS IN GRAPHENE
To describe the graphene-insulator heterostructure, we consider a model Hamiltonian consisted of a graphene part, an insulator substrate part, and the coupling between them (see Eqs. (5) and Sec.S6 of Supplementary Material [22]).As we are interested in the low-energy electronic properties, graphene's band structure is modelled by the low-energy Dirac cones around the K and K ′ valleys.The long-wavelength EC (charge ordered) state in the substrate is considered as a charge insulator, with the electrons being frozen to form a superlattice [22].Thus, long-wavelength charge order of the substrate is coupled to the graphene layer via interlayer Coulomb interactions to exert a superlattice potential on the Dirac electrons.Neglecting the intervalley coupling thanks to the large superlattice constant L s (⪆ 50 Å) [23], we can construct an effective single-particle Hamiltonian for the continuum Dirac fermions in graphene that are coupled with a superlattice Coulomb potential [22] where σ µ are the Pauli matrices (µσ x , σ y ) with the valley index µ = ±1, v F is the non-interacting Fermi velocity of graphene, and U d (r) is the background superlattice potential with the period U d (r) = U d (r + L s ).The superlattice of the EC is set to be rectangular, with anisotropy r = L y /L x and L x,y being the superlattice constant in the x, ydirection, respectively.We denote L s = L x .As a result, the superlattice potential U d (r) would fold Dirac cones into its small Brillouin zone (BZ), forming subbands and opening up a gap at the boundary of the supercell BZ, as shown in Fig. 1(c) for a rectangular superlattice with r = 1.2 (same as that of CrOCl atomic lattice) in valley K (µ = 1) with L s = 600 Å.The energy degeneracies from folding are all lifted by U d , whose Fourier component reads [22] U where Q ̸ = 0 is the reciprocal lattice vector associated with L s , Ω 0 = L x L y is the area of the primitive cell of the superlattice.The Coulomb potential U d , screened by a dielectric constant ϵ r , decays exponentially in the reciprocal space ∼ exp(−Qd), where d is the distance between the substrate surface and graphene monolayer.Furthermore, the Fermi velocities near the Dirac points of the subbands are suppressed by U d [24] as clearly shown in Fig. 1(c).Such a continuum-model description is adopted throughout the paper given that L s ≫ a (a = 2.46 Å is graphene's lattice constant) is always fulfilled for low carrier density ⪅ 10 13 cm −2 , with L s ∼ 1/ √ n for the EC state.While it is highly desirable to open a gap at the Dirac points in graphene for the purpose of field-effect device fabrication, the superlattice potential of Eq. ( 2) alone cannot gap out Dirac points in graphene as the system still preserves C 2z T symmetry.However, the Dirac points can be unstable against e-e Coulomb interactions (with the spontaneous breaking of C 2z T symmetry) once the Fermi velocity of the non-interacting band structure is suppressed below a threshold, which can be assisted by the superlattice potential from the long-wavelength charge order.One of the similar illustrations is twisted bilayer graphene (TBG) [25], where the Fermi velocity is strongly suppressed around the "magic angle", leading to moiré flat bands exhibiting diverse correlated and topological phases [26][27][28][29][30][31].Here we further calculate the Fermi velocity of the superlattice subbands around the Dirac point, denoted as v F (L s , ϵ r ), which depends on both the superlattice constant L s and the background dielectric constant ϵ r .Accordingly, the effective fine structure constant α(L s , ϵ r ) = e 2 /(4πϵ 0 ϵ r ℏv F (L s , ϵ r )) can also be tuned by L s and ϵ r , as shown in Fig. 1(d).We see that there is a substantial region in the (L s , ϵ r ) phase space with α(L s , ϵ r ) > α c ≈ 0.92 [32], which indicates that the Dirac-semimetal phase of graphene may no longer be stable against e-e interactions within this regime according to previous theoretical study [32].Such a picture is not unique to rectangular superlattice, but applies to various superlattice geometries.Treating the superlattice potential U d (Q) using second-order perturbation theory, the renormalized non-interacting effective Hamiltonian for arbitrary superlattice geometry can be expressed as We see that the effective non-interacting Hamiltonian as well as the Fermi velocity have similar dependence on L s and ϵ r (through U d (Q)) for all lattice geometries.We have also calculated the effective fine-structure constants α(L s , ϵ r ) = e 2 /(4πϵ 0 ϵ r ℏv F (L s , ϵ r )) for both triangular and square lattices [22], and the results are very similar to that of rectangular lattice with r = 1.2 shown in Fig. 1(d).This motivates us to include e-e interactions in the graphene layer in our model.Despite several theoretical predictions of gapped Dirac states in graphene [3,[6][7][8]33], to the best of our knowledge no gap at the CNP has been experimentally observed in suspended graphene yet [34,35].This can be attributed to interaction-enhanced Fermi velocity around the CNP, screening of e-e interactions due to ripple-induced charge puddles, disorder effects, etc. [3,[36][37][38][39].Nevertheless, analogous to TBG, the subbands in our system with reduced non-interacting Fermi velocity would quench the kinetic energy and further promote the e-e interaction effects in graphene.
Our unrestricted Hartree-Fock calculations [22] confirm precisely the argument above.As interaction effects are most prominent around the CNP, we project the Coulomb interactions onto only a low-energy subspace including three valence and three conduction subbands (n cut = 3) that are closest to CNP for each valley and spin.To incorporate the influences of Coulomb interactions from the high-energy remote bands, the renormalized Fermi velocity within the low-energy subspace can be derived from the renormalization group (RG) approach [2][3][4]40] where α 0 = e 2 /(4πϵ 0 ℏv F ) is the ratio between the Coulomb interaction energy and kinetic energy, i.e., the effective fine-structure constant of free-standing graphene, E * c delimits the low-energy window within which the unrestricted Hartree-Fock calculations are to be performed, and E c is an overall energy cut-off above which the Dirac-fermion description to graphene is no longer valid.Unlike TBG [41], other parameters of the effective Hamiltonian (Eq.( 1)) such as U d , are unchanged under the RG flow [22].We first study the interaction effects of graphene coupled to a rectangular superlattice potential with r = 1.2 and 50 Å ≤ L s ≤ 400 Å, corresponding to carrier density of the EC state at the surface of the substrate 0.1×10 12 cm −2 ≤ n ≤ 6.58×10 12 cm −2 (with n = 2/(rL 2 s )), with ϵ r = 3, 4, and d = 7 Å (obtained from first principles density functional theory calculations for one particular commensurate CrOCl-graphene supercell [22]).Here, we consider two different filling factors: exactly at the CNP (ν = 0) and a slight hole doping (ν ≈ −0.003).When ν = 0, a gap can be opened up due to interaction effects [see Fig. 2(a,b)], leading to two nearly degenerate insulating states, one is σ z -sublattice polarized and the other is characterized by the order parameter τ z σ z , where τ z and σ z denote the third Pauli matrix in valley and sublattice space, respectively.Then, intervalley Coulomb interactions would split such degeneracy, and the sublattice polarized insulator with zero Chern number becomes the unique ground state [22].Notably, the gap decreases almost linearly with n as shown in Fig. 2(d), and eventually vanishes as n → 0. This is because the superlattice Coulomb potential exerted on graphene is proportional to the carrier density of the long-wavelength order from the substrate.Consequently, the Fermi velocity of the bare Dirac dispersion of graphene would be less suppressed at smaller carrier density n, which disfavors gap opening.Eventually in the limit of n → 0, with a charge ordered state of infinite lattice constant, graphene would recover its non-interacting behavior as a gapless Dirac semimetal.To verify our theory, we have also experimentally measured the gaps at CNP in graphene-CrOCl heterostructure at different nominal carrier densities [22] using the same high-quality device reported in Ref. [19].The details for the measurement set up and the device configuration are presented in Sec.S8 of Supplementary Material.The measured gaps also decrease linearly with n tot , from 7.7 meV with n tot = 3.4× 10 12 cm −2 , to 5.8 meV with n tot = 0.5×10 12 cm −2 [22], consistent with the trend from theoretical calculations, as shown in Fig. 2(e).Nevertheless, when n tot → 0, such a linear dependence of the gap on n tot may no longer be valid [22].This is because in Eq. ( 2), the interlayer Coulomb potential only applies to the situation of a single valley to accommodate charge carriers in the substrate.In reality, there may be additional valley degeneracy in the substrate, which is crucial for the evolution of gap as n tot → 0. Although the valley degeneracy of the substrate does not change our results qualitatively, the theoretically calculated gap vs. n tot fits to the experimental data of CrOCl-graphene heterostructure more precisely at low density once including the two-fold valley degeneracy of CrOCl (see Table I).The details are given in Supplementary Material (Fig. S11) [22].
We note that the electronic crystal at the surface of the substrate is expected to persist even if the carrier density exceeds the threshold value due to the extra energy gain from interlayer Coulomb coupling in such coupled system, which will be discussed in detail in the section "Cooperative coupling between graphene and substrate" below.Strain is also inevitable in such graphene-insulator heterostructures, which would give rise to pseudo-magnetic fields coupled to the Dirac electrons [5,42,43], thus further enhance the e-e interaction effects in graphene.
The single-particle excitation spectrum is also significantly altered by Coulomb interactions within the low-energy window, as shown in Fig. 2(b) and (c) with fillings ν = 0 and ν = −0.003,respectively.We note that although the superlattice potential U d suppresses Fermi velocity in graphene [see Fig. 1(c)], e-e interactions can compensate such effects.The Fermi velocity is not only enhanced by the Coulomb potentials from the remote energy bands [Eq.( 4)], but also further boosted by e-e interactions within the low energy window E * c ∼ n cut ℏv F 2π/L s .Eventually, the Fermi velocity can be magnified up to more than twice of the non-interacting value of free-standing graphene (v F ) at slight hole doping ν = −0.003,as shown in Fig. 2(d).This perfectly explains the recent experiment in gate-controlled graphene-CrOCl heterostructure, in which the Fermi velocity around CNP is significantly enhanced compared to noninteracting value at slight carrier doping, such that robust quantum Hall effect can be observed under tiny vertical magnetic fields (∼ 0.1 T) and at high temperatures [19].We note that the EC state may be stabilized by vertical magnetic fields even when the carrier density in the substrate exceeds the zero-field threshold value [44,45], which in turn boosts the low-field, high-temperature quantum Hall effect in the graphene layer due to the scenario discussed above.
Although it has been theoretically proposed that the magnetic proximity effect together with spin-orbit coupling could in principle give rise to topologically nontrivial states in graphene [46], it seems to be irrelevant to the grapheneinsulator heterostructures considered in the present study.For example, in CrOCl-graphene device, no magnetic hysteresis has been observed in graphene, and the measured Landau level degeneracy is still compatible with that of spin-valley degenerate Dirac cones [19].Most saliently, the gap opening and the robust quantum Hall effect persist up to temperatures far above the Néel temperature of CrOCl (∼14 K) [19].Similarly, the magnetic proximity coupling was also reported to be negligible for CrI 3 -graphene heterostructure [21].Therefore, compared to the power-law decaying interlayer Coulomb coupling, the exponentially decaying magnetic proximity coupling may not play an important role in such charge-transfer graphene-insulator heterostructures.
The essential results discussed above, i.e., the gap opening at CNP and the concomitant drastic enhancement of Fermi velocity, remain valid for different types of the background superlattices.Specifically, we have also performed calculations for the case of triangular superlattices, which lead to qualitatively the same conclusions, as presented in Sec.S5 of Supplementary Material [22].

TOPOLOGICAL PROPERTIES
Different from magic-angle TBG [47][48][49][50][51], the low-energy subbands for graphene coupled to a rectangular superlattice potential U d (r) with small anisotropy (r ∼ 1) turn out to be topologically trivial with a compensating Berry-curvature distribution, leading to zero Chern number.This remains true even in the gapped Dirac state after including e-e interactions, as shown in Fig. 2(f).The trivial band topology is somehow anticipated because the superlattice potential is non-chiral in the sense that it is coupled equally to the two sublattice of graphene, which does not have any pseudo-gauge-field structure such as that in TBG [51,52].
Hence, it is unexpected that changing the anisotropy r and the lattice size L s of the superlattice potential U d can make the subbands topological.For example, keeping L x = 50 Å but with r = 3.0, both the highest valence band and the lowest conduction band acquire nonzero valley Chern numbers C = ±1 (after adding an infinitesimal C 2z -breaking staggered sublattice potential).As shown in Fig. 3(a), besides the four high symmetry points, it appears another two "hot spots" (annotated by green circles) along the line connecting Γ s and X s .This additional contribution breaks the balance between positive and negative contribution of Berry curvature to Chern number, leading to non-zero valley Chern number.Such contribution stems from another accidental crossing point between the low-energy valence and conduction bands along the k x -direction through changing merely the anisotropy parameter r, as shown in Fig. 3(c) by red dot within green circle.
While increasing r from unity (with fixed L s ), the Fermi velocity in the x-direction of the valence band around the Dirac point, v x , is gradually reduced, as shown in Fig. 3(e).As the same origin of Klein tunneling effects, the spinor structure of graphene's wavefunction forces the Fermi velocity in the y-direction to be intact [24].Further tuning r at some point would totally flatten v x .In Fig. 3(e), we mark by white dashed lines "the magic lines" on which v x of the valence band closest to Dirac points vanishes exactly.The magic lines always come in pair as an effect of chiral (particle-hole) symmetry breaking induced by the superlattice potential.As particle-hole symmetry is broken in the energy spectrum, when v x vanishes in the valence band, the counterpart in the conduction band remains finite.The valence subband around the Dirac point has to curve upwards to create an accidental band crossing point, after that v x of the valence band becomes zero again.Therefore, a band crossing would be germinated at the Dirac point, and then move away along the k x -direction with larger r.On the one hand, the band crossing moving away from Γ s is of accidental nature, which is generally avoided unless the lattice parameters are at some fine-tuned values.On the other hand, the Dirac point at Γ s remains stable as protected by C 2z T symmetry.If the Dirac point is gapped, say, by a tiny staggered sublattice potential, the low-energy subbands become topological with nonzero valley Chern numbers.In particular, with the increase of r at fixed L s , the absolute value of valley Chern number of the valence subband (closest to Dirac points) increases by 1 whenever one pair of the magic lines are passed through.The positions of these magic lines also depend on the background dielectric constant ϵ r since larger ϵ r corresponds to weaker Fermi-velocity renormalization effect, which would shift the magic lines to larger r values.In Supplementary Material, we provide animated figure (Fig. S4) and videos demonstrating the evolution of the band structures and Berry curvatures with increasing r at fixed L s .Such topologically nontrivial subbands with highly anisotropic Fermi velocities may provide an alternative platform to realize topological quantum matter.
We note that the anisotropic charge ordered superlattices may be realized in two ways.First, one can design a spatially modulated electrostatic potential, which has been realized in monolayer graphene by inserting a patterned dielectric superlattice between the gate and the sample [53].Then, the anisotropy of the superlattice can be artificially tuned by the dielectric patterning in the substrate.Second, for some given carrier density, the Fermi surface of the conduction (or valence) band of the substrate may be (partially) nested, which may lead to a charge density wave (CDW) state with the nesting wavevector.For example, for CrOCl, the Fermi surfaces under different Fermi energies (above the conduction band minimum) are given in Fig. S17(c) of Supplementary Material.Clearly, under some proper fillings, the Fermi surfaces are nested or partially nested, which may give rise to CDW states with anisotropic superlattices.We note that topologically nontrivial flat bands have also been proposed to exist in Bernal bilayer graphene coupled with a background superlattice potential [54].
Furthermore, we find that changing L s is also able to control the valley Chern number of the subbands.For example, with r = 3 and L s = 600 Å, as shown in Fig. 3(b), while the highest valence band remains topological with non-zero valley Chern number 1 for valley K with the two aforementioned crossing points (green circles) merely moving to X s , the lowest conduction band turns out to be topologically trivial.This is due to two additional band crossing points (orange circles) close to the Y s -S s line between the lowest and the second lowest conduction bands, as annotated by red dots in an orange circle in Fig. 3(d).
The nontrivial topology must arise from the intrinsic Berry phases of the Dirac cones.Such topologically nontrivial bands are particularly surprising for our system, since the Dirac fermions are subjected to a "trivial" superlattice potential, which couples identically with two sublattices of graphene.Nevertheless, the nontrivial subband topology is highly tunable by changing the superlattice's size and anisotropy [22].

COOPERATIVE COUPLING BETWEEN GRAPHENE AND SUBSTRATE
In the previous calculations, a charge ordered superlattice in the substrate is presumed, which exerts a classical superlattice Coulomb potential to graphene.However, this assumption should be re-examined.Moreover, besides the effects from the substrate to graphene, the feedback effects from graphene to the substrate should be discussed as well.Therefore, in this section, we study the coupled bilayer system as a whole, and treat the electrons in graphene layer and the substrate layer on equal footing.In particular, we model the carriers transferred to the substrate as 2D electron gas with long-range e-e Coulomb interactions.Electrons in the substrate and in graphene interact with each other via long-range Coulomb potential, whose Fourier component of wavevector q reads e 2 exp(−|q| d)/(2ϵ 0 ϵ r |q|).Thus, the total Hamiltonian for the Coulomb-coupled graphene-insulator heterostructure system includes [22]: On the graphene side, Eq. (5a) is the familiar Dirac Hamiltonian describing the non-interacting low-energy physics of graphene.The e-e Coulomb interactions within graphene are described by Eq. (5c), where the dominant intravalley long-range Coulomb interactions are considered and V int (q) is in the form of double-gate screened Coulomb potential (see Eq. ( 9).Here, ĉσµα (k) and ĉ † σµα (k) denote annihilation and creation operators for the low-energy Dirac electrons with wavevector k, valley µ, spin σ, and sublattice α.Note that S refers to the total surface area of the coupled system, and the atomic wavevectors k, k ′ , q are expanded around the Dirac points.On the substrate side, without loss of generality, we suppose that the chemical potential is close to the conduction band minimum (CBM) with its energy E CBM , and the energy dispersion of the low-energy electrons around CBM can be modelled by a parabolic band as for 2D free electron gas with effective mass m * .Other electrons in the deep valence bands are supposed to be integrated into the static dielectric screening constant thanks to a large gap of the substrate.Therefore, the noninteracting Hamiltonian Eq. (5b) for electrons in the substrate can be written in the plane wave basis with creation and annihilation operators { d † σ (k), dσ (k)}, where k is the plane wave wavevector expanded around the CBM, and σ denotes spin.The e-e Coulomb interactions within substrate [Eq.(5d)] is taken to be the long-range Coulomb interaction with the same double-gate screened form of V int (q).The coupling between graphene and substrate is only via the long-range Coulomb potential, which is captured by Eq. (5e).The prefactor e 2 exp(−|q| d)/(2ϵ 0 ϵ r |q|) in front of the field operators in Eq. (5e) is nothing but the 2D Fourier transform of 3D Coulomb potential.Interlayer hoppings can be neglected given that the interlayer distance d ⪆ 5 Å in such heterostructures (e.g., d ≈ 7 Å in graphene-CrOCl heterostructure from first principles calculations), thus the exponentially decaying interlayer hopping amplitude is much weaker than the power-law-decaying interlayer Coulomb interaction.This is further confirmed by directly calculating the orbital projected band structures of a commensurate supercell of CrOCl-graphene heterostructure based on density functional theory.It turns out that the Dirac cone in such heterostructure supercell stems almost 100% from carbon p z orbitals of graphene (see Sec. S7 of Supplementary Material), which clearly indicates the absence of interlayer hybridization (hopping).
We use distinct letters to denote the ladder operators for electrons in graphene (ĉ, ĉ † ) and substrate ( d, d † ).This implies in a notational manner the approximation of distinguishable electrons.In other words, the many-body wavefunction of the coupled bilayer system (denoted as |Ψ⟩) can be written a separable fashion, namely a direct product of graphene's and substrate's part, i.e., In a mean-field treatment, the corresponding many-body wavefunction would thus be a direct product of two Slater determinants, |Ψ⟩ c and |Ψ⟩ d for the graphene layer and the substrate layer, respectively.This is reminiscent of the Born-Oppenheimer approximation for electrons and ions.Technically, this means that order parameters ∼ ⟨ĉ † d⟩ (⟨ d † ĉ⟩) are not allowed in our treatment.A finite value of ⟨ĉ † d⟩ (⟨ d † ĉ⟩) suggests the emergence of another phase, an interlayer excitonic condensate in such coupled bilayer system.However, we note that such interlayer exciton has to be driven by intervalley Coulomb scattering between the K/K ′ valley of graphene and (presumingly) Γ valley of substrate's electrons, with the amplitude ∼ e 2 exp(−|K|d)/(2ϵ 0 ϵ r |K|) being several orders of magnitudes smaller than the intravalley one in our problem.Thus, it is completely legitimate to neglect the interlayer particle-hole exchange in our problem, and the separable wavefunction ansatz Eq. ( 6) is very well justified.Then, we solve the full interacting Hamiltonian Eqs.(5) under the separable-wavefunction ansatz Eq. ( 6), and the workflow is presented in Methods section.Nevertheless, the interlayer excitonic insulator state consisted of Dirac electrons (holes) and quadratically dispersive holes (electrons) is possible in valley-matched graphene-insulator heterostructures, such as those consisted of graphene and transition metal dichalcogenides with the band extrema at K points.We leave this for future study.
To explore how the interlayer Coulomb coupling would affect the electronic crystal state of the substrate, we first consider the situation as a reference that the substrate is decoupled from graphene.The energy difference between the spin polarized EC state and Fermi-liquid (FL) state (condensation energy) as a function of the carrier density n is given by quantum Monte Carlo calculations in [55,56], as shown by the green line in Fig. 4(c), where an effective mass m * = 1.3, a background dielectric constant ϵ r = 4, and valley degeneracy of 2 are considered in order to mimic the conduction band minimum of CrOCl.The condensation energy reaches zero when n ≈ 4.5 × 10 12 cm −2 (corresponding to critical Wigner-Seitz radius r * s ≈ 32.9), suggesting the transition from the EC to the FL state.More details are given in Methods section.
We further include the interlayer Coulomb coupling between the substrate and graphene (setting the chemical potential at the CNP of graphene), which can be treated using perturbation theory given that the interlayer Coulomb energy is always much smaller than the sum of the intralayer Coulomb energy and kinetic energy within the relevant parameter regime (see Fig. 5 in Methods).Specifically, with the separable wavefunction ansatz (Eq.( 6)), the groundstate charge densities for the graphene layer and the EC layer are separately obtained from unrestricted Hartree-Fock calculations, which are further used to estimate the interlayer Coulomb energy.More details about the perturbative treatment of interlayer Coulomb interactions are presented in Sec.S6 of Supplementary Material [22].
We find that the condensation energy (per electron) of the EC is substantially enhanced after including the interlayer interactions, as shown by the orange diamonds in Fig. 4. As a result, the EC-FL transition is postponed to a much higher density n ≈ 16×10 12 cm −2 (corresponding to critical Wigner-Seitz radius r * s ≈ 17.3).This is because the energy of the coupled bilayer can be further lowered by pinning the charge centers (marked as light blue stars in Fig. 4(a,b)) of the two layers in an anti-phase interlocked pattern, in order to optimize the repulsive interlayer Coulomb energy.The extra energy gain from such interlocking of charge centers compensates the energy cost of the EC state when n ⪆ 4.5 × 10 12 cm −2 , thus substantially stabilizes the EC state.
On the one hand, since the condensation energy of the free 2D electron gas in the decoupled substrate is estimated using the model that accurately fits to quantum Monte Carlo data [55], the estimate of the critical density for the decoupled substrate is expected to be accurate.On the other hand, in the case of substrate coupled with graphene layer, although the interlayer Coulomb energy is estimated with Hartree-Fock approximation, the qualitative conclusion (that the EC state gets stabilized by a cooperative interlayer Coulomb coupling) is expected to be valid even in a beyond-mean-field treatment.This is because under the separable-wavefunction ansatz, the interlayer Coulomb energy in the EC state is always negative (compared to that of FL state) under an optimal choice of relative charge centers, which thus always stabilizes the EC state even if the intralayer interactions are treated using beyond-mean-field approaches.
We note that the stabilizing effect of EC is not unique to band-aligned graphene-insulator heterostructures considered in this work.In general, it only requires the presence of another (tunable) gapped state exhibiting non-uniform charge distribution atop of the EC.For example, remarkably robust EC state has been observed in a bilayer system consisting of two monolayer MoSe 2 separated by hexagonal boron nitride [57], which was also argued to be stabilized by the interlocking of the EC states in the two layers.

MATERIALS REALIZATION
The scenario discussed above is not only closely related to CrOCl-graphene and CrI 3 -graphene heterostructures [19,21], but can also be extended to various band-aligned graphene-insulator heterostructures.As along as the conduction band minimum (CBM) or valence band maximum (VBM) of the substrate is energetically close to the Dirac points of graphene, charges could be easily transferred between graphene and the substrate's surface by gate voltages.Furthermore, it is more likely to form long-wavelength ordered state at the surface of the substrate (with TABLE I. Candidate substrate materials for the graphene-insulator heterostructure systems.The dielectric constants ϵr [58][59][60], conduction band minimum position (ECBM), valence band maximum position (EVBM), the corresponding effective mass m * at the band edge that is closer to the Dirac point (set to zero) in energy, and the dimensionless Wigner-Seitz radius rs = gvm * / √ πnϵrm0a 0 B (a 0 B is the Bohr radius and m0 is the bare electron mass, gv is the valley degeneracy) estimated under a small doping concentration n = 10 12 cm −2 , are presented.Here "bi" and "tri" stand for bilayer and trilayer systems, respectively.slight carrier doping) if the material has large effective masses at the CBM or VBM.Meanwhile, an insulator with relatively small dielectric constant would have weaker screening effects to e-e interactions, which also favours longwavelength ordered state at small carrier doping.

Materials
Following these guiding principles, we have performed high-throughput first principles calculations based on density functional theory for various insulating van der Waals materials.Eventually we find twelve suitable candidate materials (including CrOCl and CrI 3 ), whose CBM and VBM energy positions, dielectric constants (ϵ r ), effective masses at the band edges, and the corresponding Wigner-Seitz radii (r s ) are listed in Table I.Clearly, the Wigner-Seitz radii of these materials at the band edges (estimated under slight doping concentration n = 10 12 cm −2 ) are all above the threshold of forming a Wigner-crystal state (r s ⪆ 31) [55].Additionally, the energy bands of these insulating substrate materials can be easily shifted using vertical displacement fields [22], such that charge transfer between graphene and the substrate can be controlled by non-disruptive gate voltages.We have also considered heterostructures consisted of graphene and TMDs.Besides trilayer (or thicker) WS 2 as already listed in Table I, we further nominate WSe 2 (trilayer or thicker), MoSe 2 (bilayer or thicker), and MoTe 2 (bilayer or thicker) as possible candidate substrates to realize the effects discussed above.More details are given in Sec.S7 of Supplementary Material.

CONCLUSIONS
In conclusion, we have studied the synergistic correlated electronic states emerging from coupled graphene-insulator heterostructures with gate tunable band alignment.Based on comprehensive theoretical studies, we have shown that the gate tunable carrier doping may yield a long-wavelength electronic crystal at the surface of the substrate driven by e-e interactions within the substrate, which in turn exerts a superlattice Coulomb potential to the Dirac electrons in graphene layer.This would substantially change the low-energy spectrum of graphene, where a gapped Dirac state concomitant with drastically enhanced Fermi velocity would emerge as e-e interaction effects.These theoretical results are quantitatively supported by our transport measurements in graphene-CrOCl heterostructure.Besides, the Dirac subbands in graphene can be endowed with nontrivial topological properties by virtue of the interlayer Coulomb coupling with the long-wavelength electronic crystal underneath.Reciprocally, the electronic crystal in the substrate can be substantially stabilized by virtue of a cooperative interlayer Coulomb coupling with the gapped Dirac state of graphene.We have further performed high-throughput first principles calculations, and suggested a number of promising insulating materials as candidate substrates for graphene to realize such effects.
However, the understanding of such coupled bilayer correlated electronic systems is still at a preliminary stage, and the study is far from being complete.First, the long-wavelength electronic crystal cannot be the only possible candidate ground state.Other correlated states such as magnetic or even superconducting states may also occur in the charge doped insulating substrate, e.g., in the case of high-temperature cuprate superconductor [15,16] and monolayer 1T'-WTe 2 [61].This may give rise to diverse quantum states of matter in graphene due to interfacial proximity couplings with Dirac fermions.Moreover, so far we have only considered the ground state properties of such coupled bilayer correlated electronic systems.What is more intriguing is the collective excitations of the electronic crystal and their couplings with Dirac electrons in graphene.Around the quantum melting point of the electronic crystal, strong quantum fluctuations would be coupled with Dirac fermions with graphene via interlayer Coulomb interactions, which may give rise to unique quantum critical properties.Therefore, our work may stimulate further exploration of the intriguing physics in such a platform for correlated and topological electrons.

Hartree-Fock approximations assisted by renormalization group approach
When graphene is coupled to a superlattice potential, the Coulomb interactions are suitably expressed in the subband eigenfunction basis, on which we have performed the Hartree-Fock calculations.Since interaction effects are most prominent around the CNP, we project the Coulomb interactions onto only a low-energy window including three valence and three conduction subbands that are closest to the Dirac point per valley per spin.We use a mesh of 18 × 18 k-points to sample the mini Brillouin zone of the superlattice.
To incorporate the influences of Coulomb interactions from the high-energy remote bands, we rescale the Fermi velocity within the low-energy window of the effective Hamiltonian using Eq. ( 4).The other parameters of the non-interacting effective Hamiltonian are unchanged under the RG treatment since their corrections are of higher order, thus can be neglected.In other words, we find the following RG equations for Fermi velocity v F and leading superlattice potential The detailed derivations of the RG equations are presented in Sec.S3 of Supplementary Material.We also neglect on-site Hubbard interactions and intervalley coupling in e-e Coulomb interactions, which turn out to be one or two order(s) of magnitude weaker than the dominant intravalley long-range Coulomb interactions in such graphene-based superlattice systems [62].To model the screening effects to the e-e Coulomb interactions from the dielectric environment, we introduce the double gate screening form of V int , whose Fourier transform is expressed as where Ω 0 is the area of the superlattice's primitive cell, ϵ r is a background dielectric constant and the thickness between two gates is d s = 400 Å.Then, we initialize the Hartree-Fock loop with the initial conditions in the form of various different order parameters and obtain the converged ground state self-consistently (see Sec. S4 of Supplementary Material [22]).
When we consider electrons in graphene and substrate on equal footing in Eq. ( 5), the routine of Hartree-Fock calculations is exactly the same.However, we need to first consider solely the substrate side.After performing unrestricted Hartree-Fock calculations, we use the ground-state charge density of EC in the substrate as input for constructing the superlattice potential.Explicitly, we need to replace Eq. ( 2) by where ρ d (Q) is the Fourier component of the charge density in the substrate.More details can be found in Sec.S6 of Supplementary Material [22].

Workflow to solve the coupled bilayer Hamiltonian
We solve the Hamiltonian of the coupled bilayer system described by Eqs.(5) in the following workflow: • First, we start our calculations by considering solely the substrate Hamiltonian Eqs.(5b) and (5d).We considered the case of triangular superlattice, which is the actual ground state for the EC of 2D electron gas.In particular, the total energy of the triangular EC can described by a fitting model given in Ref. [55]: where c 1 = −1.106103,c 3/2 = 0.814, c 2 = 0.113743, c 5/2 = −1.184994and c 3 = 3.097610.These parameters are determined by fitting to quantum Monte Carlo data.The total energy for the Fermi-liquid state of 2D electron gas is given by the following model [56]: where x = √ r s and with with a 0 = −0.1925, a 1 = 7.3218, a 2 = 0.16008 and a 3 = 3.1698.These parameters for the FL state are also determined by fitting to quantum Monte Carlo data [56].The energies are given in Hartree atomic units.Then, one can extract the condensation energy for the isolated 2D electron gas in the substrate E WC − E FL , with the accuracy comparable to quantum Monte Carlo calculations.
• Second, with the help of the separable wavefunction ansatz Eq. ( 6), we further calculate the ground-state charge density of the EC state in the substrate under Hartree-Fock approximations.Although the Wigner crystal condensation energy would be significantly overestimated with such mean-field approximation, the ground-state charge density can still be properly described by the unrestricted Hartree-Fock treatment [63].Then, one can integrate out the charge degrees of freedom of the substrate so that the charge density modulation characterized by the Fourier components of the charge density {ρ d (Q)} (Q denotes the reciprocal vector of the superlattice) can be used as an input for the superlattice potential U d (Q), as shown in Eq. (10).Compared to Eq. ( 2), this superlattice potential is more realistic and self-contained in our model.Eq. ( 10) would be recovered to Eq. ( 2) by setting ρ d (Q) = 2 for any reciprocal vector Q, which is equivalent to say that two (spin degenerate) charges per primitive supercell are localized in real space in a Dirac-δ-function form.
• Third, we perform RG-assisted unrestricted HF calculations for the interacting Dirac electrons in graphene as explained in Methods .If the chemical potential is at the CNP of graphene, a gap opening will be triggered by e-e interactions within the graphene layer as discussed previously.However, the ground states are obtained so far by minimizing (mostly) the intralayer parts of the full Hamiltonian, the interlayer Coulomb interaction Eq. (5e) is not optimized yet.We note that the intralayer kinetic energy and intralayer Coulomb interaction energy for both graphene and the substrate are unchanged under constant lateral shifts of the charge centers, thus the ground state |Ψ⟩ d ⊗ |Ψ⟩ c obtained so far is massively degenerate up to global and relative shifts of the bilayer charge centers.Such degeneracy would be partially lifted by the interlayer Coulomb energy ⟨H gr−sub ⟩.Obviously, ⟨H gr−sub ⟩ is invariant under the global shift of the charge centers of the bilayer system, but it varies with respect to a relative charge-center shift.Therefore, by virtue of perturbation theory, optimizing the interlayer Coulomb energy amounts to find the optimal relative shift vector between the charge centers of the two layers within the degenerate ground-state manifold obtained in the previous procedures.Such perturbative treatment of H gr−sub is justified given that the interlayer Coulomb energy is always weaker than the sum of the kinetic energy and the intralayer Coulomb energy within relevant parameter regime, as shown in Fig. 5.For example, the interlayer Coulomb energy ∼ 20 meV for typical parameters L s = 50 Å and ϵ r = 4, while the intralayer Coulomb energy ∼ 60 meV.More details for the perturbative calculation of interlayer Coulomb energy can be found in Sec.S6 of Supplementary Material [22].
• Finally, we gather all the contributions from Eq. ( 5) to find out the total energy of the coupled bilayer system staying in a gapped Dirac state (at the CNP) for graphene and a long-wavelength EC state for the substrate.By comparing it with that of a non-interacting Dirac state for graphene and a 2D Fermi-liquid state for the substrate, we can then find out if the gapped graphene interplays with the long-wavelength charge-ordered substrate in a cooperative or competitive manner.
It turns out that the bilayer system tends to cooperate with each other such that both the gapped Dirac state (at the CNP) of graphene and the long-wavelength charge ordered state in the substrate are substantially stabilized by the interlayer Coulomb coupling.The results are presented in Fig. 2 and 4 of the main text.

Density functional theory calculations
The first principles calculations are performed with the projector augmented-wave method within the density functional theory [64], as implemented in the Vienna ab initio simulation package software [65].The crystal structure is fully optimized until the energy difference between two successive steps is smaller than 10 −6 eV and the Hellmann-Feynman force on each atom is less than 0.01 eV• Å.The generalized gradient approximation by Perdew, Burke, and Ernzerhof is taken as the exchange-correlation potential [66].As Cr is a transition metal element with localized 3d orbitals, we use the on-site Hubbard parameter U = 5.48 eV for the Cr 3d orbitals in the CrOCl bilayer and U = 3 eV for Cr 3d orbitals in the CrI 3 bilayer.The so-called fully localized limit of the spin-polarized GGA+U functional is adopted as suggested by Liechtenstein and coworkers [67], and the non-spherical contributions from the gradient corrections are taken into consideration.The "DFT+D2" type of vdW correction has been adopted for all multilayer calculations to properly describe the interlayer interactions [68].
Our high-throughput filtering of the proper insulating substrate materials for graphene starts from the 2D materials computational database [69].We only focus at those with bulk van der Waals structures which have been previously synthesized in laboratory.This ensures that it is experimentally feasible to exfoliate few layers from their bulk sample and then stack them on graphene to form heterostructures.
Experimental measurements of the gaps in graphene-CrOCl heterostructure By designing a dual-gated structure, we used few-layered CrOCl as an bottom dielectric while few-layered hexagonal boron nitride (h-BN) was served as top gate dielectric.The top and bottom gate voltages can then be converted into doping and displacement fields for further data analysis.Graphene, h-BN, and CrOCl flakes are mechanically exfoliated from high quality bulk crystals.The vertical assembly of few-layered hBN, monolayer graphene and fewlayered CrOCl were made using the polymer-assisted dry-transfer method.Electron beam lithography was done using a Zeiss Sigma 300 SEM with a Raith Elphy Quantum graphic writer.Top and bottom gates as well as contacting electrodes were fabricated with an e-beam evaporator, with typical thicknesses of Ti/Au ∼ 5/50 nm.Electrical transport measurements of the devices were performed using an Oxford TeslaTron 1.5 K system.Gate voltages on the as-prepared multi-terminal devices were fed by a Keithley 2400 source meter.Channel resistances were recorded in 4-probe configurations using low frequency (13.33 Hz) lock-in technique with Stanford SR830 amplifiers.The gate dependencies of channel resistances were measured at various temperatures for the extraction of thermal gaps.More details about the device configuration, measurement set up, and sample quality can be found in Sec.S8 of Supplementary Material.
Supplementary Material for "Synergistic correlated states and nontrivial topology in graphene-insulator heterostructures"

S1. NON-INTERACTING HAMILTONIAN FOR A GRAPHENE-INSULATOR HETEROSTRUCTURE
The Hamiltonian for a graphene-insulator heterostructure can be always divided into three parts: graphene part H G , the insulating substrate part H S and the coupling between them H G−S .The graphene part can be suitably described by a tight-binding model since we focus on the low-energy physics.As we have explained in the main text, with slight carrier doping band edge, the insulator substrate is supposed to form a long-wavelength charge order on the interface near graphene sheet thanks to Coulomb interactions between electrons occupying the band edge of the insulating substrate (transferred from graphene layer).The insulator substrate part is then modeled by a 2D Hamiltonian for electrons hopping on a 2D superlattice which forms a Wigner-crystal-like or long-wavelength charge ordered insulator state at some proper filling.The geometry of the superlattice is determined by the long-wavelength order at the interface.Explicitly, the graphene part H G and the insulator substrate part H S can be generally written as where ĉ( †) σα (k) and d( †) σ ( k) are fermionic annihilation (creation) operators for electrons in graphene and the insulator substrate, respectively.In the lower index of these operators, α is the sublattice index for the bipartite lattice of graphene and σ is the spin degree of freedom of electrons.To emphasize the fact that graphene and the insulator substrate have different lattices and thus different Brillouin zone, we denote k and k as the wavevectors in the Brillouin zone of graphene and that of the long-wavelength superlattice in the substrate, respectively.In our calculations, the lattice for H S is set to be rectangular or triangular.The lattice symmetry turns out to make no qualitative differences.
Since electrons have negligible probability to hop between graphene and the insulator substrate due to rather large distance d between two sheets in the z-direction (d ∼ 7 Å from DFT calculations in CrOCl-graphene heterostructure), we suppose that electrons from two sheets are coupled only via long-ranged Coulomb interactions.Unlike H G and H S , such long-ranged Coulomb interactions are more easily written in real space.In terms of field operators ψ(r), the inter-sheet coupling reads where V (|r − r ′ + dẑ|) is the 3D long-ranged Coulomb potential e 2 /4πϵ 0 ϵ r r and electrons in graphene and the insulating substrate are described by the field operators with lower index c and d, respectively.Here ϵ 0 is the vacuum permittivity and ϵ r is the dimensionless relative dielectric constant of the insulating substrate.In the spirit of tightbinding formalism, we write the field operators in terms of Wannier functions where ϕ α and ϕ are Wannier functions localized on the graphene and the insulator substrate Bravais lattice sites, which are described by a i and R i , respectively.Here α refers to the sublattice index in graphene and τ α is the vector denoting the position of the αth sublattice inside the unit-cell.The spin degrees of freedom is included by the index σ and also explicitly by spinor χ σ .It is worthwhile to note that here the Bravais lattice sites and the corresponding Wannier functions for the substrate refer to those of the spontaneously generated charge ordered superlattice, not the atomic lattices of the substrate.A more fundamental treatment based on the atomic lattice sites of the substrate lattice will be discussed in Sec.S6.
After the transformations above, the Hamiltonian H G−S in the Wannier basis reads with If Wannier functions are so localized such that with δ (2) (r) is the 2D Dirac δ-function distribution, we can simplify the previous expression to with δ µ,ν is the Kronecker delta and Then, we write H G−S in reciprocal space using the following Fourier transformation where N c and N d are the number of lattice sites for electron in graphene and the insulator substrate, respectively.
The Hamiltonian H G−S in the basis of ĉσα (k) and dσ ( k) reads Now we first define R = a i − R j , and let k ′ − k = q = q + G, where G is a reciprocal vector of the long-wavelength ordered superlattice and q is the wavevector within the superlattice Brillouin zone.Then we take use of the identity The coupling V ( q + G) reads where Ω d is the area of the unit-cell of the surface superlattice of the substrate.In the third line of the above derivation, we smear the sum over R = a i − R j by replacing it with an integral over the surface S = N d Ω d = N c Ω c with Ω c the area of graphene's unit-cell since we are interested in the physics in the length scale of the superlattice {R j }, which is supposed to much larger than that of graphene.Finally, the last line is the 2D partial Fourier transformation of the 3D Coulomb potential.Since we focus on the low-energy physics around the Dirac cones of graphene, we can attribute valley index µ to electrons in graphene and neglect intervalley coupling thanks to the exponential decay of V (q) so that In the meantime, the Hamiltonian for graphene only H G [see Eq. (S1)] can be divided into two valley sectors where σ µ = (µσ x , σ y ) with σ x,y are the Pauli matrices and the valley index µ = ±1.
In the Hartree approximation by contracting c and d fermion operators separately, we have Since the long-wavelength charge order state is insulating presumingly with two spin degenerate electrons occupying each supercell, we have Writing k = k + G with G in the superlattice reciprocal lattice, the final form of the coupling between graphene and insulating substrate used in our calculations reads where In the meantime, we integrate out the Hamiltonian for insulating substrate H S [see Eq. (S2)] so that it becomes a constant charge density, which is omitted in our calculations.To wrap up, we get the effective non-interacting Hamiltonian in continuum in the valley µ where the Fourier component of U d (r) is precisely U d (G) [see Eq. (S24)] with G in the reciprocal lattice of the underlying insulating substrate's surface superlattice.
As revealed by Eq. (S24), increasing the interlayer distance d would diminish the interlayer Coulomb interaction between the charges of the Wigner crystal at the surface of the substrate and the Dirac electrons in the graphene layer [see Eq. (2) in the main text].To estimate the e-e interaction effects in graphene with different interlayer distance d, we have calculated the effective fine structure constant of graphene: α = e 2 /4πϵ 0 ϵ r ℏv SL F , where the Fermi velocity v SL F is the non-interacting Fermi velocity of the graphene being coupled to the Coulomb superlattice potential.Note that v SL F is reduced compared to the free-standing non-interacting Fermi velocity v F due to the effects of Coulomb superlattice potential, whose strength depends on the interlayer distance d, the background dielectric constant ϵ r , and the superlattice constant L s of the Wigner crystal formed at the surface of the substrate.In Fig. S1(a), we plot the effective fine structure constant α as a function of the interlayer distance d, and the superlattice constant L s , with a fixed background dielectric constant ϵ r = 3.We note that there is small region in the upper left corner where α ≤ α c = 0.92, which means that graphene would remain as a semimetal in this region.We further plot α vs. d in Fig. S1(b), and find that α decreases with the increase of d, and α becomes smaller than the critical value 0.92 when d > 20 Å, the moment that is expected to undergo a transition from a gapped to a gapless Dirac state.It turns out that the colormaps of α remains qualitatively the same for different lattice geometries, comparing Fig. S1(a) for rectangular lattice with Fig. S2 for triangular and square lattice.This also justifies why we can consider only two particular cases, rectangular and triangular lattice, in our work without loss of generality.In our numerical implementations, the lattice of insulating substrate is set to be rectangular or triangular, from which we obtain qualitatively the same correlated states in the graphene layer.The range of {G i } is limited to |n x |, |n y | ≤ 4 with G = n x g x + n y g y .g x,y are the two reciprocal lattice vectors for the rectangular lattice of insulating substrate.The sum over Q in Eq. (S23) stops at the limit |n x | + |n y | ≤ 2.

S2. TOPOLOGICAL PROPERTIES OF THE NON-INTERACTING HAMILTONIAN
We further study the topological properties of our model Hamiltonian derived in the previous section.Different from magic-angle TBG [47][48][49][50][51], the low-energy subbands for graphene coupled to a rectangular superlattice potential U d (r) with small anisotropy (r ∼ 1) turns out to be topologically trivial.To be specific, in Fig. 2d of the main text, we have shown the Berry curvature distribution of the highest valence band of valley K in the mini BZ of the superlattice with L x = 50 Å and r = L y /L x = 1.2.We see that the Berry curvature is mostly concentrated at the band crossing points, i.e., the four high symmetry points Γ s , S s , X s , and Y s .The contributions from the Γ s and the S s points are exactly compensated by those from the X s and Y s points, resulting in a band with zero Chern number.This is anticipated because the superlattice potential is non-chiral in the sense that it is coupled equally to the two sublattice of graphene.This remains true even including e-e interactions.
Hence, it is unexpected that changing the anisotropy r and the lattice size L s of the superlattice potential U d can make the subbands topological.For example, keeping L x = 50 Å but with r = 3.0, the valley Chern number of the low-energy subbands become nonzero.For valley K, the highest valence band and lowest conduction band now have a Chern number C = ±1, respectively.As shown in Fig. S3(a), besides the four high symmetry points, it appears another two "hot spots" (annotated by green circles) along the line connecting Γ s and X s .This new contribution breaks the balance between positive and negative contribution of Berry curvature to Chern number, leading to non-zero valley Chern number.Such contribution stems from a new crossing point between the low-energy valence and conduction bands along the k x -direction through changing merely the anisotropy parameter r, as shown in Fig. S3(c) with red dot in green circle.From the animated Fig. S4, we see that while increasing r, the Fermi velocity is gradually reduced, and at r = 2.7, becomes vanishingly small along one direction due to Klein tunneling effects.Then, further tuning r germinates a band crossing originated from the Dirac point then gradually moving away from it.Alternatively, we find that changing L s can also control the valley Chern number of the subbands, since L s is encoded in the superlattice potential (see Eq. (S24)).For example, with r = 3 and L s = 600 Å, as shown in Fig. S3(b), while the highest valence band remains topological with non-zero valley Chern number 1 for valley K with the two aforementioned crossing points (green circles) merely moving to X s , the lowest conduction band turns out to be topologically trivial.This is due to two new band crossing points (orange circles) close to the Y s -S s line between the lowest and the second lowest conduction bands, as annotated by red dots in an orange circle in Fig. S3(d).Such topologically nontrivial bands are particularly surprising for our system, since the Dirac fermions are coupled to a "trivial" superlattice potential coupled identically on two sublattices, different from the case of magic-angle TBG.Thus, the nontrivial topology must arise by virtue of the intrinsic Berry phases of the Dirac cones.We further show the non-interacting energy spectra and distributions of Berry curvature in the first Brillouin zone for various superlattice constants (L s ) and anisotropy parameters (r), i.e., L s = 50, 200, 600 Å with r = 1.2 and 3.The plots are given in Fig. S5, S6, S7 for L s = 50, 200, 600 Å, respectively.Since the system preserves timereversal symmetry and the superlattice potential U d is diagonal in the sublattice subspace, the non-interacting energy spectrum in valley K ′ is exactly the same as that in valley K so that we only plot the spectrum for valley K here.As shown below, the distribution of Berry curvature of the highest valence and the lowest conduction band in valley K is exactly opposite to that in valley K ′ as another consequence of time-reversal symmetry of the system.
We also provide separately six videos in other Supplementary Data, which show the non-interacting energy spectra and distributions of Berry curvature in the first Brillouin zone for L s = 50, 200, 600 Å with r = 1-10, respectively.We note that the anisotropic charge ordered superlattices may be realized in two ways.First, one can design a spatially modulated electrostatic potential, which has been realized in monolayer graphene by inserting a patterned dielectric superlattice between the gate and the sample [53].Then, the anisotropy of the superlattice can be artificially tuned by the dielectric patterning in the substrate.Second, for some given carrier density, the Fermi surface of the conduction (or valence) band of the substrate may be (partially) nested, which may lead to a charge density wave (CDW) state with the nesting wavevector.For example, for CrOCl, the Fermi surfaces under different Fermi energies (above the conduction band minimum) are given in Fig. S14(c).Clearly, under some proper fillings, the Fermi surfaces are nested or partially nested, which may give rise to CDW states with anisotropic superlattices.We note that topologically nontrivial flat bands have also been proposed to exist in Bernal bilayer graphene coupled with a background superlattice potential [54].

S3. RENORMALIZATION GROUP DERIVATIONS
The derivation shown in this section is inspired from Ref. 41.The e-e Coulomb interaction operator in our derivations is written as where V c (r) = e 2 /4πϵ 0 ϵ r r and ρ(r) is the density operator of electrons at r.The Hamiltonian Eq. ( S25) is defined at some high energy cut-off ±E c .We focus in the valley µ = +1 by the virtue of which the derivation for the valley µ = −1 is immediate and the results are identical.Remember that the parameters v F and U d (G) should be thought of as being fixed by a measurement at E c without e-e interactions.This also amounts to ρ(r) = ψ † (r) ψ(r) with the non-interacting field operator ψ(r) where ϕ σnk (r) is the wavefunction of an eigenstate of the non-interacting Hamiltonian H 0 [see Eq. (S25)] with energy ϵ n,k and its associated annihilation operator is ĉσn (k).

Electron-electron interaction in a lower energy window
Now we change the cut-off E c to a smaller one E ′ c and see how these parameters are modified by Vint .Vint can be treated perturbatively when E ′ c is much larger than any other energy scale in the system.To do so, we split the field Then, we integrate out the fast modes ψ> (r) in the expansion of ρ(r)ρ(r ′ ).Note that ψ> (r) and ψ> † (r) must appear equal times in each terms of the expansion otherwise it would vanish by taking the non-interacting mean value ⟨. . .⟩ 0 .Explicitly, these terms are retained up to a constant: The first term gives the Coulomb e-e interaction between electrons of the slow modes ψ< (r) below the new cut-off E ′ c .The second and third term could be omitted if the system has particle-hole (p-h) symmetry as in twisted bilayer graphene [41].In our system described by Eq. (S25), the first nearest-neighbor coupling in U d (G) preserves particlehole (p-h) symmetry.The p-h symmetry is broken if further-neighbor coupling is included, which is exponentially smaller [see Eq. (S24)].So, it is legitimate in our RG derivation to neglect such weak p-h asymmetry in order to omit the second and the third term in the expansion.Then, we evaluate the rest of the terms in the expansion, which represents precisely the correction to H 0 from the fast modes ψ> (r) via Coulomb e-e interactions.Let us write where the minus sign in the second line comes from the exchange the two fermionic operators and the constant arising from the exchange is omitted.Then, the e-e interaction Vint in the lower energy window delimited by Evaluation of the correction to the non-interacting Hamiltonian from the fast modes In the following, we set ℏ = 1 for the simplicity in mathematical expressions.Note that F(r, r ′ ) has the structure of the residue of the Green's function Ĝ(z) = (z − Ĥ0 ) −1 taking only the valley µ = +1 part in Ĥ0 = ĤG + ĤG−S [see Eqs.(S20) and (S23)], namely where the contour C encloses the z-plane real line segment [−E c , −E ′ c ] in the clockwise, and segment [E ′ c , E c ] in the counterclockwise, sense.As long as E ′ c dominates over all other energy scales such as U d (G 0 ) and v F G 0 with G 0 denoting the primitive reciprocal vector of the underlying superlattice, the dominant contribution to the contour integral can be evaluated perturbatively using Ĝ(z) It is easier to calculate the Green's function in the plane wave basis |k⟩ where δ (2) (x) is the 2D Dirac distribution.In the plane wave basis, the evaluation of Green's functions is straightforward Then, the contour integral can be easily done: Renormalization group flow equations Now we only have to insert the previous results into the second term in Eq. (S33) to derive the RG equations for v F and U d (G).Let us compute first the integral for ⟨r| Ĝ0 (z)|r ′ ⟩.After writing the 2D Coulomb potential in Fourier space V 2D (q) = e 2 /2ϵ 0 ϵ r q, we have is the Fourier transform of ψ< (r).Since v F q ≪ E ′ c , we can Taylor expand (A) in terms of q/k.The leading order reads Therefore, the RG equation reads Actually, we find the famous result of the Fermi velocity renormalization in graphene due to the e-e interactions.
In the same way, we calculate the integral for ⟨r| Ĝ0 (z) ĤG−S Ĝ0 (z)|r ′ ⟩: , we can Taylor expand (B) in terms of q/k and G/k (considered as if they have the same order of magnitude).The leading order reads which can be neglected under the first-order RG procedure, namely In summary, we have shown that the Fermi velocity in graphene v F is renormalized by the e-e Coulomb interaction in the standard way while the superlattice potential U d (r) keep its value unchanged.In our numerical study of e-e interactions, we use the renormalized Fermi velocity v * F in the Hartree-Fock calculations, where we have to take a cut-off n cut to the number of bands, to include the contributions from the higher energy bands outside the cut-off.Technically, we use where L s and a 0 are the lattice constant of the superlattice of U d (r) and the carbon-carbon bond length in graphene, respectively.Here, the ratio L s /n cut a 0 plays the role of E c /E ′ c .

S4. HARTREE-FOCK APPROXIMATIONS TO ELECTRON-ELECTRON INTERACTIONS
The derivation shown in this section is inspired from Ref. 62.We consider the Coulomb interactions in graphene where ψσ (r) is real-space electron annihilation operator at r with spin σ.This interaction can be written as where Here i, α, and σ refer to Bravis lattice vectors, layer/sublattice index, and spin index.ϕ is Wannier function and χ is the two-component spinor wave function.We further assume that the "density-density" like interaction is dominant in the system, i.e., Here we can see that the Coulomb interaction can be divided into intersite Coulomb interaction and on-site Coulomb interaction.Given that the electron density is low (10 11 cm −2 ), i.e., a few electrons per supercell, the chance that two electrons meet at the same atomic site is very low.The Coulomb correlations between two electron are mostly contributed by the inter-site Coulomb interactions.Therefore, the on-site Hubbard interaction has been neglected in our calculations.
In order to model the screening effects to the e-e Coulomb interactions from the dielectric environment, we introduce the double-gate screening form of V int , whose Fourier transform is expressed in Eq. ( 9).Since we are interested in the low-energy bands, the intersite Coulomb interactions can be divided into the intra-valley term and the inter-valley term.The intra-valley term V intra can be expressed as with N s is the total number of the superlattice's sites.The inter-valley term V inter is expressed as V intra includes the Coulomb scattering processes of two electrons created and annihilated in the same valley, and V inter includes the processes that two electrons are created in µ and −µ and get annihilated in −µ and µ valleys.
Here the atomic wavevector k is expanded around the valley K µ in the big Brillouin zone of graphene, which can be decomposed as k = k + G, where k is the superlattice wavevector in the superlattice Brillouin zone, and G denotes a superlattice reciprocal lattice vector.The electron annihilation operator can be transformed from the original basis to the band basis: where C σµαG,n ( k) is the expansion coefficient in the n-th Bloch eigenstate at k of valley µ: We note that the non-interacting Bloch functions are spin degenerate due to the separate spin rotational symmetry (SU (2) ⊗ SU (2) symmetry) of each valley.Using the transformation given in Eq. (S56), the intra-and inter-valley Coulomb interaction can be written in the band basis and where the form factors Ω µσ,µ ′ σ ′ nm,n ′ m ′ and Ω µ,σσ ′ nm,n ′ m ′ are written respectively as and We make Hartree-Fock approximation to Eq. (S58) and Eq.(S59) such that the two-particle Hamiltonian is decomposed into a superposition of the Hartree and Fock single-particle Hamiltonians, where the Hartree term is expressed as The Fock term is expressed as: and We note that the typical intravalley interaction energy ∼ 240 meV for L s = 50 Å and ϵ r = 3; while the intervalley interaction ∼ 30 meV, which is one order of magnitudes smaller than the intravalley interaction, thus we neglect the intervalley term [see Eq. (S55)] in most of our calculations.We also check a posteriori that the intervalley Hartree and Fock energies are at least two orders of magnitude smaller than their intravalley counterpart.However, the intervalley interaction is crucial to lift the degeneracy between many-body ground state, as shown in the following section.

S5. RESULTS OF HARTREE-FOCK CALCULATIONS
In this section, we gather the results of Hartree-Fock calculations including Hartree-Fock single-particle spectra and distributions of Berry curvature in the first Brillouin zone for L s = 50, 200, 600 Å.
First, we show the Hartree-Fock single-particle spectrum with a superlattice potential with r = 1.2 of L s = 50, 200, 600 Å in Fig. S8.Here, we use n cut = 5 and study three types of doping: CNP (ν = 0), slight hole doping (ν = −0.003)and slight electron doping (ν = +0.003).As you can see from Table I and the Hartree-Fock single-particle spectra, the results of a slightly electron-doped system is similar to those for a slightly hole-doped one.Note that we include only intravalley Coulomb interactions in these calculations.As shown in the following, the role of intervalley Coulomb interactions is merely to lift the ground state degeneracy and favor the σ z -state.We further calculate the interacting electronic structure of graphene at different filling factors (denoted by ν) away from the CNP, as presented in Fig. S10.We find that the ground state is generally gapless at nonzero fillings, and the Fermi velocity enhances as the chemical potential approaches the charge neutrality point (|ν| → 0).In particular, the Fermi velocity increases from Now we show the effect of intervalley Coulomb interactions by comparing the total energy of σ z -state with τ z σ zstate for different L s = 50, 200, 600 Å with fixed r = 1.2.We calculate the difference (always negative) between them and see how it changes when we include the intervalley Coulomb interactions.Here, we use n cut = 3.
As you can see from Table II, the energy difference between the total energy of σ z -state and τ z σ z -state is enhanced by two orders of magnitude for L s = 50 and 200 Å.However, the energy difference for L s = 600 Å does not benefit anything from intervalley interactions.
We also have performed Hartree-Fock calculations on a triangular lattice including three valence and three conduction bands (n cut = 3) for L s = 50, 200, 600 Å using 18 × 18 k-mesh in the BZ.As shown in Table III, the results on a triangular lattice are qualitatively the same as those on a rectangular lattice.This ensures that our conclusions are lattice-independent.
TABLE III.Parameters extracted from the Hartree-Fock single-particle spectra on a triangular lattice: gap opened at the CNP (ν = 0) and the ratio between interaction-renormalized Fermi velocity v * F and the non-interacting one vF for different Ls = 50, 200, 600 Å.

Ls( Å) 50 200 600
Gap at ν = 0.0 (meV) 21 1.9 0.24 v * F /vF at ν = −0.0032.2 1.7 1.7 In some material, there are several valleys to accommodate the charges transferred from graphene so that we need to consider this valley degeneracy in the superlattice potential as follows: where n(Q) is the Fourier transformed charge density per valley at superlattice's reciprocal vector Q and the valley degeneracy g v = 2 for CrOCl.After including this valley degeneracy, our theoretical results are quantitatively consistent with the experimental data, as shown in Fig. S11.More details on experimental measurements are given in the last section Sec.S8.

S6. GENERAL COUPLED BILAYER SYSTEM IN GRAPHENE-INSULATOR HETEROSTRUCTURES
Previous parts have already proved that the presence of a superlattice potential underneath monolayer graphene sheet helps the latter to open a gap at the CNP, and concomitantly enhance the Fermi velocity around the Dirac point.Now we would like to reconsider the assumptions made for the purpose of writing a simplistic Hamiltonian of our coupled bilayer system given by Eqs.(S1),(S2), and (S3).

Hamiltonians
First, when the Fermi level in graphene aligns closely to the band edge (say, the conduction band minimum) of the insulating substrate, the density of electrons in the conduction band has been assumed to be so low that electrons spontaneously break translational symmetry and form a charge-ordered superlattice.Second, we have assumed that the quantum degrees of freedom on the substrate are completely frozen so that the only effect (other than dielectric screening) acting by substrate to electrons in graphene is to apply a superlattice potential via the long-range Coulomb interactions, stemming from a non-uniform charge density of the charge-ordered state.Lastly, we have drastically supposed that the Wannier functions and the distribution of charge density at the surface of the substrate is as localized as a Dirac-δ-function like distribution in real space [see Eq. (S11)].These assumptions were good enough if we focus on the physics on the graphene side.In particular, it turns out that the assumption on the extreme localization of the charges and/or Wannier functions [see Eq. (S11)] in the substrate is a good starting point.This is because more extended charge distributions only make the Fourier components of the superlattice potential weaker (to be discussed in the following), which would not change our conclusion qualitatively.However, if we are interested in especially the synergistic interplay between graphene and the insulating substrate, we have to consider integrally the graphene-insulator heterostructure as a coupled bilayer system, and treat both layers on equal footing.
In this section, we aim to study, to the best of our efforts, the coupled bilayer system as a whole, and to treat quantum mechanically the electrons both in the graphene layer and the substrate layer.Most of the above assumptions have to be discarded or replaced by a less harsh one.To this end, in addition to the kinetic energy Eq.(S1) and the long-range intralayer Coulomb interactions of graphene Eq. (S50), we also need to consider the interacting many-body Hamiltonian for electrons on the surface of the substrate.Without loss of generality, we consider the situation that the conduction band minimum (CBM) of the substrate is charge doped.
Specifically, the non-interacting kinetic Hamiltonian of the low-energy electrons around the CBM of the substrate can be modeled by the following one where E CBM is the energy position of CBM with respect to the Fermi level of graphene, and m * is the effective mass.For CrOCl, m * ≈ 1.3m e .d( †) σ (k s ) is the annihilation (creation) operator of the low-energy electrons doped to the surface of the substrate, with the wavevector k s (expanded around the CBM) in the atomic Brillouin zone of the substrate and spin index σ.Then we consider the long-range Coulomb interactions between electrons at the surface of the substrate: where N s is the total number of atomic lattice sites in the substrate layer, Ω s is the area of the atomic primitive cell in the substrate layer, and k s , k ′ s , q s , denote the atomic wavevectors in the Brillouin zone of the atomic lattice of the substrate.The Coulomb potential V int (q) is given by Eq. (9).
We continue to address the interlayer Coulomb interaction, i.e., Eq. (S3).We first transform Eq. (S3) to the basis of Wannier functions of the atomic lattices of graphene and the substrate using the following transformations of the field operators: and the Fock term reads where we recenter Q to get the last line.Here, ⟨. . .⟩ d means the expectation value of observable after integrating only the substrate part of the many-body ground state |Ψ⟩ 0 d .Spin polarization is allowed in the Hartree-Fock treatment of the 2D electron gas in the substrate, and indeed a spin polarized Wigner crystal state always has a lower energy than a spin degenerate one.
The HF calculations have been done for a series of different superlattice constant L s .Once L s is given, the charge density in substrate is then fixed by n d = g v /( √ 3L 2 s /2) for a spin polarized Wigner crystal state with triangular superlattice, where g v = 2 comes from the valley degeneracy of the conduction band minimum of CrOCl.Here we set up our parameters including the effective mass m * = 1.3m 0 , the background dielectric constant ϵ r = 4, and valley degeneracy g v = 2, to mimic the properties of conduction band minimum of CrOCl.Such mapping amounts to always let only the lowest subband in the folded Brillouin zone be filled.During the iterations, we keep track all the subbands.To initialize the HF self-consistent loop, the terms like d are set to be non-zero, corresponding to a spontaneous charge order with Fourier component Q, where Q is one of the primitive reciprocal vectors for the triangular superlattice.In other words, we start the HF loop from a charge-ordered state, which facilitates the convergence to the same phase.If the final converged state is gapped, we can extract its charge modulation as input for the next step.In our calculations, a 9 × 9 mesh of reciprocal lattice (centered at Γ point) has been adopted, with the mini Brillouin zone sampled by a 18 × 18 k mesh.The former corresponds to the real-space mesh within a primitive cell, while the latter corresponds to the system size.We have also performed calculations adopting a 13 × 13 mesh of reciprocal lattice points (corresponding to finer real-space mesh), and find completely consistent results with those calculated using a 9 × 9 reciprocal lattice points.
With the HF results on the substrate side in hands, we are ready to study the rest of three Hamiltonians, namely Eqs.(5a), (5c) and (5e) of the main text.The mean-field treatment is exactly the same as in Sec.S3 and Sec.S4.However, the interlayer coupling Hamiltonian Eq. (5e) worth special attention.Under the separable wavefunction ansatz, we can still integrate out all the d-operators using the wavefunctions resulted from the previous HF calculations solely on the substrate side: Here, we define the charge modulation ρ d (Q) of the charge-ordered state as where we write the expectation values in the substrate's subband basis in the presence of HF potentials {D σG,n ( k)} relates precisely the d-operators in the original plane-wave basis to that in the Bloch-function basis, and E d n k ′ denotes the subband energy dispersion.Remarkably, the two conceptually different treatments give rise to an interlayer Coulomb potential with the same analytical properties.It justifies our previous drastic assumptions on the localization of charge distribution.Here we only do better since we can self-consistently solve the charge-ordered ground state of the substrate and its charge density distributions.
Exactly parallel with what we have done earlier, we now solve the problem of graphene on the superlattice defined ρ d (Q) while only keeping the dominant terms with small Q.The e-e interactions are still treated in the band basis by HF approximations in a low-energy window given by the band index cut-off.We keep track only three valence and three conduction bands per valley per spin.The contributions from high-energy bands are taken into account using the RG approach, replacing the Fermi velocity by a renormalized one.In the continuum model treatment, a 9 × 9 mesh of reciprocal lattice points has been adopted, and the mini Brillouin zone is sampled in a 18 × 18 k mesh.The filling of graphene is set to be at the CNP so that a sublattice gap would be opened in the presence of superlattice potential.Once found the final converged ground state using all the types of initial order parameters, we calculate the charge modulation of gapped graphene ρ c (Q) where ⟨. . .⟩ c is the expectation value of integrating only the graphene part of the many-body ground state |Ψ⟩ 0 c .So far, we have separately found the HF ground state for substrate only and that for graphene coupled to the superlattice potential arising from the charge-ordered state in substrate.We note ⟨. . .⟩ HF as the expectation value of observable using the HF ground state, i.e., charge-ordered state for substrate and gapped state for graphene.As a reference, we also note ⟨. . .⟩ FL as the expectation value of observable using the non-interacting plane-wave state, i.e., 2D electron gas for substrate and non-interacting Dirac fermions for graphene.If the charge-ordered state cooperates with the gapped state in graphene, the coupled bilayer system must be energetically more stable than the substrate alone.Explicitly, this means that the condensation energy E cond, sub > E cond, coupled , where E gr-sub,opt is the interlayer Coulomb energy resulted from non-uniform charge modulations of the two layers, after optimization with respect to relative charge-center shift.It will be treated using perturbation theory shortly.Besides, there are two subtleties in this comparison.First one is related to the cut-off.When we study graphene in the presence of superlattice potential, only the bands in low-energy window are considered.In principle, we need to add the energy contribution of all the electrons in the occupied bands.However, this would not make a significant difference since electrons deep in the valence bands are little affected by all the interaction effects discussed here.Particularly, the relevant interaction energy scale considered in this work is always far below our choice of low-energy cutoff for graphene (∼ 3 × ℏv F 2π/L s ).On the other hand, the contribution from the remote energy bands outside the low-energy window has already been included in our RG approach.Our second concern comes from the feedback effect from graphene to substrate.In principle, the non-uniform charge modulation of graphene would also affect the profile of the charge density distribution of the electronic-crystal state in the substrate by minimizing the interlayer Coulomb interaction.However, such interlayer Coulomb potential is weak compared to the intralayer counterpart for the substrate, as shown in Fig. 5 in the main text.This allows us to treat the feedback effect perturbatively, i.e., include ⟨H gr-sub ⟩ c as a perturbation to H 0 sub + H intra sub .The first order effect is precisely the charge interlocking between graphene and substrate.There exists an optimal relative shift vector t m between the two layers that minimize the interlayer Coulomb energy cost However, to the second order, one has to consider how the HF wavefunction of the electronic-crystal state in substrate is changed by the charge modulation of graphene.In the HF subband basis, The energy shift due to such modification is where Q ̸ = 0 in the above equation.Note that E gr-sub is also sensitive to the relative shift between the two layers so one need to find the optimal t opt such that it minimizes the sum of E gr-sub , namely E gr-sub,opt in Eq. (S86).This is the value we need to use in the comparison between two condensation energy E cond, sub and E cond, coupled .We also note that all the kinetic energies and intralayer Coulomb interaction energies are unchanged under the relative shift between the two layers.The results will answer the question on how synergistic correlated states occur in such coupled bilayer heterostructure.

S7. DETAILS OF DFT CALCULATIONS FOR THE SUBSTRATE MATERIALS
Lattice structures, deformation potentials, and band structures of candidate substrate materials In this section and the following one, we present the details for the density function theory (DFT) calculations of the 14 candidate substrate materials presented in Table I of the main text.In this subsection, we show the results for all but CrOCl and transition dichalcogenides.These materials, which worth special attention, are given in the next subsections.The first principles calculations are performed with the projector augmented-wave method within the density functional theory [64], as implemented in the Vienna ab initio simulation package software [65].The crystal structure is fully optimized until the energy difference between two successive steps is smaller than 10 −6 eV and the Hellmann-Feynman force on each atom is less than 0.01 eV/ Å.The generalized gradient approximation by Perdew, Burke, and Ernzerhof is taken as the exchange-correlation potential [66].As Cr is a transition metal element with localized 3d orbitals, we use the on-site Hubbard parameter U = 5.48 eV for the Cr 3d orbitals in the CrOCl bilayer and U = 3 eV for Cr 3d orbitals in the CrI 3 bilayer.The so-called fully localized limit of the spin-polarized GGA+U functional is adopted as suggested by Liechtenstein and coworkers [67], and the non-spherical contributions from the gradient corrections are taken into consideration.The "DFT+D2" type of vdW correction has been adopted for all multilayer calculations to properly describe the interlayer interactions [68].
Our high-throughput filtering of the proper insulating substrate materials for graphene starts from the 2D materials computational database [69].We only focus at those with bulk van der Waals structures which have been previously synthesized in laboratory.This ensures that it is experimentally feasible to exfoliate few layers from their bulk sample and then stack them on graphene to form heterostructures.The results are summarized in Table I of the main text.
The lattice structures of some of the substrate materials presented in Table I of the main text are shown in Fig. S12.They are either in the bilayer or trilayer structures.The lattice structure of CrI 3 is similar to that of YI 3 as shown in Fig. S12(d).Their band structures are presented in Fig. S13, where the green dashed lines mark the energy position of the Dirac point in graphene.We note that the valence band maximum (VBM) of PbO bilayer is energetically close to the Dirac point of graphene; while for the other bilayer or trilayer substrate materials, their conduction band minima (CBM) are close to the Dirac point.This indicates that charge transfer can easily occur between graphene and the substrates controlled by gate voltages.Moreover, we note that the conduction bands and valence bands of these materials are typically flat with large effective masses, which would be very susceptible to e-e Coulomb interactions once these substrate materials are slightly charge doped, and may lead to Wigner-crystal-like state or long-wavelength ordered state as discussed in main text.Another important precondition for the Wigner-crystal state is that the screening effect of substrate materials can not be too strong.For example, the conduction band of ScOBr bilayer has a large effective mass of 2.575m 0 (m 0 is the bare mass of a free electron), but the dielectric constant ϵ r of ScOBr reaches ∼13, which makes it difficult to trigger the Wigner-crystal-like instability in this material under slight charge doping.
We note that all of these proposed substrate materials all have been successfully synthesized in laboratory as listed in Table .IV. Especially, few-layer of ReSe 2 as a highly anisotropic material [70][71][72], and few-layer CrI 3 system as a 2D magnetic material [73][74][75][76], have been extensively studied recently.Moreover, phonon spectra calculations have proved the dynamical stability of these substrate materials in monolayer from [69].Thus the device fabrication of heterostructure consisting of graphene monolayer and one of these candidate substrate materials should be experimentally accessible.
There always exists tension or compression in a heterostructure system.Under some lattice deformation, the variation of conduction band minimum (CBM) or valence band maximum (VBM) is defined as deformation potential.We list the deformation potentials of the candidate substrate materials in Table .IV.We note that the maximum value of the deformation potential is only 5.84 eV for ScOCl, which means that the energy level of CBM of ScOCl would move down by only 0.063 eV under 1% tensile strain.Therefore, even if strain is introduced in the graphene-insulator heterostructure proposed in this work, the band edges (with large effective masses) of those candidate substrate materials are still energetically close to the Dirac point of graphene.
In these candidate materials (except for CrOCl), CrI 3 bilayer is the only magnetic system.Previous theoretical studies reveal that the stacking configuration of CrI 3 bilayer plays an important role in the magnetic ground state [77].
Here we use the AB ′ -type stacking in the bilayer structure, which is consistent with the stacking configuration in the bulk phase of CrI 3 .The AB ′ -stacked CrI 3 bilayer is in an intralayer ferromagnetic and interlayer antiferromagnetic ground state.We single out CrI3 by plotting its band structure using red solid lines.This is because CrI3 is the only magnetic system among these ten materials so that it worth special attention on the magnetic nature of its ground state.Here, we suppose CrI3 to have an intralayer ferromagnetic and interlayer antiferromagnetic ground state.

Band structure of graphene-CrOCl heterostructure
We study graphene-CrOCl heterostructure using DFT calculations to show the small orbital overlap between graphene's carbon atom and CrOCl's Cr atom.Technically, we need to construct a commensurate supercell matching graphene and CrOCl's primitive unit-cell.In our modeling, we use an 8 × √ 7 supercell of graphene.We refer a 1 and a 2 to the lattice vectors along the short and long axis, respectively.Then, the corresponding supercell for CrOCl is a (2a 1 , a 1 + 5a 2 ) one.The average mismatch is ∼ 0.7%.In particular, there are 0.4% tensile stress along long axis and 2.0% compressive strain along short axis in the graphene supercell.
The band structure is shown in Fig. S15.The lattice distortion amounts to adding strain to graphene and thus moves the position of Dirac cone away from high-symmetry lines.Therefore, we choose the line (from Γ to M ′ , which is not a high-symmetry point but near M ) that passes the Dirac point in the Brillouin zone to plot the band structure.We see that the orbital overlap is small and thus negligible.

Electric-field tunable band structures of bilayer CrOCl
Now we discuss the electronic structure of CrOCl bilayer under vertical electrical fields.Here we consider an intralayer ferromagnetic and interlayer antiferromagnetic state for the bilayer configuration, which turns out to be one of competing low-energy magnetic states, and is the magnetic ground state when the on-site Hubbard U value for the Cr 3d orbitals is large.[78] The calculated band gap of CrOCl bilayer with the DFT+U calculation is 3.13 eV, which is close to that of HSE06 calculation (3.12 eV) [69].The band structure of antiferromagnetic CrOCl bilayer is shown in Fig. S14(a), where the green dashed line marks the energy position of the Dirac point of graphene.Without vertical electric field, the Dirac point is slightly above the CBM of bilayer CrOCl.Applying a vertical electric field of 0.03 V/nm would push down the CBM as shown in Fig. S14(b).A closer inspection reveals that the top-layer conduction state (red lines) is pushed downwards while the bottom-layer state (blue lines) is pushed upward in energy as shown in Fig. S14(b), such that electron carriers in the graphene layer (if there is any) would be transferred to the top layer of CrOCl substrate, forming a Wigner-crystal-like state at the surface of CrOCl substrate given that the Wigner-Seitz radius of the CBM ∼ 55.7−74.2(with a relative dielectric constant ϵ r = 3−4) is above the threshold value ∼ 31 (see Table .I in the main text).Thus, our conjecture is supported by detailed first principles DFT calculations.
In Fig. S14(c) we also present the Fermi surfaces at different Fermi energies above the CBM of bilayer CrOCl.At very low carrier densities with small Fermi energy (CBM is set to zero), the Fermi surface consists of two nearly isotropic circles.For example, at filling factor 1/100 (corresponding to a carrier density ∼ 8 × 10 12 cm −2 ), the Fermi surface is marked by the red circles.Such isotropic Fermi surface with large effective mass (∼ 1.308m 0 ) is likely to give rise to Wigner-crystal state as discussed in the main text.As the Fermi level further increases, the Fermi surfaces become more and more anisotropic.Now we focus on the few layers of transition metal dichalcogenides (TMD) which have attracted remarkable interest recently.Here we use the dielectric constants of the corresponding bulk TMD materials to estimate the condition for the onset of Wigner crystal states in these few-layer systems.We find that the conduction band minima of MoX 2 (X = S, Se, Te), WY 2 and PtY 2 (Y = S, Se) few-layer systems are close to the Dirac point of graphene, as shown in the Table .V. Especially, the conduction band edges in the MoX 2 and WY 2 systems mainly originate from the transition metal Mo or W d orbitals, which are localized and may have large effective masses.If the conduction bands of these TMD systems are slightly carrier doped, the systems may form Wigner crystal states as long as their Wigner-Seitz radii r s > 31.We list the critical carrier concentration to realize the Wigner crystal states in various TMD few layers in the Table .V, and their energy bands are shown in Fig. S16.The conduction band minima in the Pt-based TMD systems are energetically close to the Dirac point of graphene, but the strong screening effect (large dielectric constants) may prevent the electrons to form Wigner crystal states in the Pt-based TMD system.The out-of-plane electric field can easily tune the electronic structure in the MoX 2 and WY 2 systems.A previous theoretical work indicated that an electric field of 0.3 V/nm can decrease the band gap by 0.2 eV [89].Because of the excellent fabrication technology of the TMD system and their CBM close to the Dirac point, we believe that the graphene/Mo(W)-based TMD heterostructure can provide a feasible platform to realize gapped Dirac state concomitant with interaction-enhanced Fermi velocities.

S8. EXPERIMENTAL MEASUREMENTS OF THE GAPS IN GRAPHENE-CROCL HETEROSTRUCTURE
To test our theory of the band reconstruction of Dirac fermions in graphene coupled with a long-wavelength charge order, we considered a few candidate substrates, among which CrOCl is suitable for device fabrication because of its high air-stability and easy-exfoliatable nature.By designing a dual-gated structure, we used few-layered CrOCl as an bottom dielectric while few-layered hexagonal boron nitride (h-BN) was served as top gate dielectric.The top and bottom gate voltages can then be converted into doping and displacement fields for further data analysis.In this section, we include below detailed experimental data on the quality of sample, the device configuration, the measurement setup as well as how we do the thermal gap measurements.

Quality of sample, device configuration and measurement setup
The devices are made of graphene, h-BN, and CrOCl flakes, which are mechanically exfoliated from high quality bulk crystals.The vertical assembly of few-layered hBN, monolayer graphene and few-layered CrOCl were made using the polymer-assisted dry-transfer method.Electron beam lithography was done using a Zeiss Sigma 300 SEM with a Raith Elphy Quantum graphic writer.Top and bottom gates as well as contacting electrodes were fabricated with an e-beam evaporator, with typical thicknesses of Ti/Au ∼ 5/50 nm.A cartoon illustration of the tested device is shown in Fig. S17(a), which includes h-BN/Graphene/CrOCl van der Waals heterostructure equipped with a top gate and a bottom gate.Fig. S17(b) and (c) show the 5× magnification and 100× magnification optical pictures of the device, respectively.There are no visible bubbles or wrinkles on the heterostructure.To ensure the cleanliness of the sample, we will use an AFM (atomic force microscope) tip to clean the Au bottom gate in a Contact mode before landing heterostructures onto these local metallic gates.A significant amount of PMMA can be removed via such a process to make sure the homogeneity of the gate electrical fields (Fig. S17(d,e) ).The device studied in this work is the same device S40 we investigated in our recent experimental paper published in Nat.Nanotechnol.17, 1272-1279 (2022) [19].According to the high quality of quantum Hall plateaus (in the normal phase without charge transfer) seen at moderate magnetic field at a few or a few tens of Kelvin in these samples, we can claim that those samples, in term of cleanness, are of state-of-the-art quality.
For electrical transport measurements of the devices, we use a wire bonder to bond the device onto the sample holder (Fig. S17(a)).Samples are then loaded into a 1.5 K fridge (Oxford Tesla-Tron system), with a base temperature of 1.5K and a maximum magnetic field of 12 T.We adopt standard 4-probe low frequency (13.33 Hz) lock-in measurement method to perform electrical transport characterization of the devices, with the Standford SR830 lock-in amplifier in series with a 10 MΩ bias resistor, thus providing an AC signal at constant current of 100 nA (assume the sample resistance is way smaller than the bias resistor).Meanwhile, high precision voltage source (Keithley 2400) is used to regulate the top and bottom gate voltages.The longitudinal voltage V xx is measured by the lock-in amplifier.Pictures of the set-up are shown in Fig. S18.
In addition, according to our experimental observations, we actually have two phases (see Fig. 4 in Ref. [19]) in graphene-CrOCl heterostructure, revealed in the sample's resistance measurement in the parameter space of V tg and V bg : 1. Phase-I is conventional graphene showing rather low Dirac peak at a few hundreds of Ohms; 2. Phase-II is the interfacial-coupling phase, where the band structure of graphene is re-configured with a mild gap opening at the charge neutral.
In fact, we focus on the Phase-II since Phase-I is the trivial conventional graphene.In the latter phase, the Wigner crystal in the substrate is not activated since the charge transfer is not onset.Given that quantum Hall plateaus have been clearly seen under sub-Tesla magnetic field at tens of Kelvin, with the Shubnikov-de Haas quantum oscillations (see Supplementary Figure 36 in Ref. [19]), the quality of sample is guaranteed.So, the measured resistivity peak has to do with gap opening in the Phase-II in CrOCl-contacted graphene.For even more experimental details, we would invite our referee to kindly refer to the Supplementary Information of Nat.Nanotechnol.17, 1272-1279 (2022) [19].

Measurement of the thermal gap
The gate dependencies of channel resistances are measured at various temperatures for the extraction of thermal gaps.In Fig. S19(a), we present the detailed mapping of measured resistance in the space of displacement field D eff and the nominal carrier density n tot .A resistivity peak which persists up to n tot ⪅ 5 × 10 12 cm −2 is clearly seen, indicating gap opening at the charge neutrality point (CNP) of graphene.The gapped Dirac state persists up to n tot ⪅ 5 × 10 12 cm −2 because the extra charge carriers are transferred from graphene to the surface of CrOCl, leaving graphene at the charge neutrality.As discussed in the main text, the charges transferred to the surface of CrOCl may form a long-wavelength ordered state through the Wigner-crystallization mechanism, which imposes a superlattice Coulomb potential to graphene and promotes the gap opening at the CNP.In order to determine the gap at the CNP of graphene at different n tot , we have further performed the temperature dependent conductance measurements

FIG. 1 .
FIG. 1.(a) Cartoon illustration of a monolayer graphene supported by an insulating substrate with long-wavelength charge order (blue dots), with an interlayer distance d.(b) Schematic of charge transfer in a band-aligned graphene-insulator heterostructure and its effects on the Dirac dispersion.(c) shows the non-interacting band structure by blue solid lines with r = 1.2 and Ls = 600 Å.The red dashed lines represent the non-interacting Dirac cones in free-standing graphene.The inset marks the high-symmetry points in the superlattice Brillouin zone.(d) shows the calculated effective fine structure constant α(Ls, ϵr), where the dashed line marks the critical value αc ≈ 0.92.

FIG. 2 .
FIG. 2. (a) Calculated Hartree-Fock single-particle excitation spectrum of graphene coupled to a long-wavelength Coulomb potential, with ν = 0. (b) and (c) show by blue solid lines the Hartree-Fock band structures of Ls = 50 Å and ϵr = 3.0, with the filling factor ν = 0 in (b) and ν = −0.003 in (c).The red dashed lines represent the non-interacting Dirac cones.The insets zoom in energy close to the Dirac points.Zero energies in (b) and (c) are defined as the chemical potentials for ν = 0 and ν = −0.003,respectively.(d) The calculated gaps at CNP (filled stars) and the interaction-enhanced Fermi velocities at slight hole dopings ν = −0.003(hollow diamonds) as a function of the substrate's carrier density n.(e) The thermal activation gap ∆ measured on the devices in [19] for different nominal dopings ntot.(f) Distribution of Berry curvature of the highest valence subband of K valley for r = 1.2 and Ls = 50 Å, which gives zero valley Chern number.

FIG. 3 .
FIG. 3. (a) and (b) shows the distribution of Berry curvature in the r = 3 superlattice's BZ of the lowest valence and conduction band in valley K for Ls = 50 and 600 Å, respectively.Their corresponding valley Chern number are also given on the top of each panel.(c) and (d) are the non-interacting band structure of the r = 3 superlattice with Ls = 50 and 600 Å. (e) Colormap of Fermi velocity in the x-direction vx of the valence band for ϵr = 3.The color coding indicates vx/vF .Here we vary Lx from 50 to 600 Å and anisotropy parameter r from 1 to 6.The white dashed line, i.e., the "magic lines", mark the position in parameter space where vx vanishes.

FIG. 4 .
FIG. 4. Charge density modulations after minimizing interlayer Coulomb interactions for (a) gapped Dirac state in graphene, and (b) electronic-crystal state in the substrate, which forms triangular superlattice.(c) Condensation energy (per electron) of the electronic crystal state E cond vs. the carrier density n in the substrate.The green line represents the condensation energy of the decoupled system using the quantum Monte Carlo data, and the orange lines shows that of the coupled bilayer system, which is significantly stabilized by interlayer Coulomb coupling.The critical rs for Wigner-crystal transitions are also indicated in the coupled and decoupled cases, respectively.

FIG. 5 .
FIG. 5. Order of magnitudes of subband width (green) and intra-(red) and inter-layer (blue) Coulomb potential strength for different Ls.The dielectric constant is selected to be 4.
FIG. S1.(a) Calculated effective fine-structure constant α of graphene in the parameter space of the background Wigner-crystal superlattice Ls, and the interlayer distance between graphene and the substrate d, with the background dielectric constant fixed as ϵr = 3.(b) Line plot of α vs. d, with Ls = 50 Å, and ϵr = 3.The red lines mark the critical value of αc = 0.92 above which the Dirac point of graphene would be gapped by e-e interactions.
FIG. S3.(a) and (b) shows the distribution of Berry curvature in the r = 3 superlattice's BZ of the lowest valence and conduction band in valley K for Ls = 50 and 600 Å, respectively.Their corresponding valley Chern number are also given on the top of each panel.(c) and (d) are the non-interacting band structure of the r = 3 superlattice with Ls = 50 and 600 Å.
FIG. S4.Animation (click the figure in, for example, Adobe Acrobat PDF reader to activate) for the non-interacting band structures of graphene on a superlattice with Ls = 50 Å and the anisotropy parameter r varying from 1 to 5.
FIG. S6.Non-interacting spectrum and distribution of Berry curvature in the first Brillouin zone for Ls = 200 Å with (a) r = 1.2 and (b) r = 3. Please refer to Sec.S2.

FIG
FIG. S7.Non-interacting spectrum and distribution of Berry curvature in the first Brillouin zone for Ls = 600 Å with (a) r = 1.2 and (b) r = 3. Please refer to Sec.S2.

FIG. S10 .
FIG. S10.The ratio between the interacting Fermi velocity (v *F ) and the non-interacting one of free-standing graphene (v 0 F ) vs. the filling factor (ν) with the superlattice constant Ls = 50 Å, and a background dielectric constant ϵr = 3.
FIG. S11.Comparison between the experimentally measured gap and the theoretically calculated gaps in graphene-CrOCl heterostructure.
FIG. S12.(a)-(f): top views of the lattice structures of some candidate substrate materials in monolayer form.The primitive cells are remarked with black lines.(g)-(l): the side views of these substrate materials in few-layer form.
FIG.S13.The calculated energy bands of the candidate substrate materials, where the energy position of the Dirac point of graphene is marked by a dashed green line in the band structures.We single out CrI3 by plotting its band structure using red solid lines.This is because CrI3 is the only magnetic system among these ten materials so that it worth special attention on the magnetic nature of its ground state.Here, we suppose CrI3 to have an intralayer ferromagnetic and interlayer antiferromagnetic ground state.
FIG. S14.The calculated energy bands of antiferromagnetic bilayer CrOCl: (a) without electric field, and (b) with an electric field 0.3 V/nm.In (b), the energy bands from top and bottom layers are marked by red and blue lines, respectively.The energy position of the Dirac point of graphene are remarked with green dashed lines.(c) Fermi surface of bilayer CrOCl at different Fermi levels with respect to the conduction band minimum.The Fermi surface under 1/100 electron filling factor is remarked by red circles.
FIG. S17.The characterizations of the experimental devices.(a) The cartoon illustration of a typical h-BN/Graphene/CrOCl heterostructure device, with its optical pictures shown in (b) and (c).Scale bar in the image is 10 µm.Images of AFM scans of the Au bottom gate before (d) and after (e) AFM Contact mode cleaning are also illustrated, where the cleaned window can be highlighted by the PMMA residues accumulated during the AFM "sweeping" cleaning.

•
From the above procedures, we would separately obtain converged HF ground states, |Ψ⟩ d for the substrate, and |Ψ⟩ c for graphene, respectively.From the ground-state wavefunctions |Ψ⟩ d and |Ψ⟩ c , one can extract the corresponding ground-state charge density modulations {ρ d (Q)} and {ρ c (Q)}, based on which the interlayer Coulomb energy [the expectation value of Eq. (5e)] can be calculated.More details are given in Sec.S6 of Supplementary Material.

TABLE IV .
[69]experimental works about the ten substrate materials, and the uni-axial deformation potentials of these materials[69].