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

<jats:title>Abstract</jats:title><jats:p>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.</jats:p>

section_content: 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 3graphene 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 grapheneinsulator 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 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 tunneling 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. 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 noninteracting Fermi velocity of the Dirac cone, thus trigger 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 non-zero 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 EC 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.

section_name: Results

section_content: 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 Supplementary Note 6 of Supplementary Information). As we are interested in the low-energy electronic properties, graphene's band structure is modeled by the low-energy Dirac cones around the K and K 0 valleys. The long-wavelength EC (charge ordered) state in the substrate is considered as a charge insulator, with the electrons being frozen in the form of a superlattice, as schematically shown in Fig. 1a. Thus, longwavelength charge order of the substrate is coupled to the graphene layer via interlayer Coulomb interactions to exert a superlattice potential on the Dirac electrons (see Fig. 1b ). If one considers that the charge order in the substrate layer results from a Wigner-crystallike instability, then the value of superlattice constant L s = 50 Å would correspond to a carrier density ~7 × 10 12 cm -2 transferred from graphene to the insulating substrate, which is close to the upper limit for a double-gated graphene device. Neglecting the intervalley coupling thanks to the large superlattice constant L s (≿50 Å), we can construct an effective single-particle Hamiltonian for the continuum Dirac fermions in graphene that are coupled with a superlattice Coulomb potential (see Supplementary Note 1 and Supplementary Note 6 in Supplementary Information)

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, y-direction, 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, forming subbands and opening up a gap at the boundary of the supercell BZ, as shown in Fig. 1c 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

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 22 as clearly shown in Fig. 1c. 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 n ⪅ 10 13 cm -2, with L s ∼ 1= ffiffiffi n p 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) 23, where the Fermi velocity is strongly suppressed around the "magic angle", leading to moiré flat bands exhibiting diverse correlated and topological phases [24] [25] [26] [27] [28] [29]. 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. 1d. We see that there is a substantial region in the (L s, ϵ r ) parameter space with α(L s, ϵ r ) > α c ≈ 0. 92 30, 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 30. 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 finestructure constants α(L s, ϵ r ) = e 2 /(4πϵ 0 ϵ r ℏv F (L s, ϵ r )) for both triangular and square lattices (see Supplementary Figure 2 ), and the results are very similar to that of rectangular lattice with r = 1. 2 shown in Fig. 1d. 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] 31, to the best of our knowledge no gap at the CNP has been experimentally observed in suspended graphene yet 32, 33. Th s 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, [34] [35] [36] [37]. 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 (see Supplementary Note 4 in Supplementary Information) 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-4,38

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 39, other parameters of the effective Hamiltonian (Eq. ( 1 )) such as U d, are unchanged under the RG flow (see Supplementary Note 3 in Supplementary Information). 

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 (see Supplementary Note 7 in Supplementary Information)). 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. 2a, 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 (see Supplementary Note 5 in Supplementary Information). Notably, the gap decreases almost linearly with n as shown in Fig. 2d, 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 using the same high-quality device reported in ref. 19. The details for the measurement set up and the device configuration are presented in Supplementary Note 8 of Supplementary Information. 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 (see Supplementary Figure 20 ), consistent with the trend from theoretical calculations, as shown in Fig. 2e. Nevertheless, when n tot → 0, such a linear dependence of the gap on n tot may no longer be valid. 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 1 ). The details are given in Supplementary Note 5 in Supplementary Information. 

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 subsection "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, 40, 41, 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. 2b 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. 1c ), 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 The dielectric constants ϵr [64] [65] [66], 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 r s = g v m * = ffiffiffiffiffiffi πn p ϵ r m 0 a 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. 

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. 2d. This explains the recent experiment in gatecontrolled 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 42, 43, 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 44, it seems to be irrelevant to the graphene-insulator 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 chargetransfer 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 Supplementary Note 5 of Supplementary Information.
section_name: Coulomb interaction effects in graphene

section_content: Different from magic-angle TBG [45] [46] [47] [48] [49], 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. 2f. 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 49, 50. 

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 v = ± 1 (after adding an infinitesimal C 2zbreaking staggered sublattice potential). As shown in Fig. 3a, besides the four high-symmetry points, it appears another two hot spots with concentrated Berry curvatures (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. 3c 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. 3e. 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 22. Further tuning r at some point would totally flatten v x. In Fig. 3e, 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 finetuned 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 non-zero 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 Information, we provide animated figure (Supplementary Figure 4 ) and videos (Supplementary Movies 1-6) 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 51. 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 Supplementary Figure 15c. 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 52. 

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. 3b, 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. 3d. 

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 (see Supplementary Note 2 in Supplementary Information).
section_name: Topological properties

section_content: 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 reexamined. 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ðÀjqj dÞ=ð2ϵ 0 ϵ r jqjÞ. Thus, the total Hamiltonian for the Coulomb-coupled graphene-insulator heterostructure system includes:

V int ðqÞ ĉy σμα ðk + qÞ ĉy σ 0 μ 0 α 0 ðk 0 À qÞ ĉσ 0 μ 0 α 0 ðk 0 Þ ĉσμα ðkÞ, ð5cÞ

e 2 e Àjqjd 2ϵ 0 ϵ r jqj ĉy σμα ðkÞ dy σ 0 ðk 0 Þ dσ 0 ðk 0 À qÞ ĉσμα ðk + qÞ:

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 ĉy σμα ð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 0,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 modeled 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 non-interacting Hamiltonian Eq. (5b) for electrons in the substrate can be written in the plane wave basis with creation and annihilation operators f dy σ ðkÞ, dσ ðkÞg, 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ðÀjqj dÞ=ð2ϵ 0 ϵ r jqjÞ 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 Supplementary Note 7 of Supplementary Information), which clearly indicates the absence of interlayer hybridization (hopping). 

We use distinct letters to denote the ladder operators for electrons in graphene (ĉ,ĉ y ) and substrate ( d, dy

). 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 Ψ j i) 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, Ψ j i c and Ψ j i 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 ∼ hĉ y di ðh dy ĉiÞ are not allowed in our treatment. A finite value of hĉ y di ðh dy ĉiÞ 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 0 valley of graphene and (presumingly) Γ valley of substrate's electrons, with the amplitude ∼ e 2 expðÀjKjdÞ=ð2ϵ 0 ϵ r jKjÞ 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 grapheneinsulator heterostructures, such as those consisted of graphene and transition metal dichalcogenides with the band extrema at K and K 0 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 refs. 53,54, as shown by the green line in Fig. 4c, 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 Supplementary Figure 12 ). Specifically, with the separable wavefunction ansatz (Eq. ( 6 )), the ground-state 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 Supplementary Note 6 of Supplementary Information. 

We find that the condensation energy (per electron) of the EC is substantially enhanced in amplitude 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. 4a, 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 53, 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 Fig. 4 | Stabilizing effect of electronic crystal by interlocking of the charge modulations in coupled graphene-insulator heterostructures. Charge density modulations n(r) in real space after minimizing interlayer Coulomb interactions for a gapped Dirac state in graphene, and b electronic-crystal state in the substrate, which forms triangular superlattice. The blue stars annotate a maximum and a minimum of density in each state. The dashed lines connect two stars located at the same in-plane coordinate. 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 r s for Wignercrystal transitions are also indicated in the coupled and decoupled cases, respectively. 

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 bandaligned graphene-insulator heterostructures considered in this work. In principle, it only requires the presence of another state exhibiting non-uniform charge distribution atop of the EC, such that the interlayer Coulomb energy gain would compensate for any energy cost of the long-wavelength charge modulations in the two layers. For example, remarkably robust EC state has been observed in a bilayer system consisting of two monolayer MoSe 2 separated by hexagonal boron nitride 55, which was also argued to be stabilized by the interlocking of the EC states in the two layers.
section_name: Cooperative coupling between graphene and substrate

section_content: The scenario discussed above is not only closely related to CrOClgraphene 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 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 favors long-wavelength ordered state at small carrier doping. 

Following these guiding principles, we have performed highthroughput 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 1. 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) 53. In addition, the energy bands of these insulating substrate materials can be easily shifted using vertical displacement fields (see Supplementary Note 7 in Supplementary Information), 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 1, 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 Supplementary Note 7 of Supplementary Information.
section_name: Materials realization

section_content: In summary, 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 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 17. 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.
section_name: Discussion

section_name: Methods

section_content: 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 lowenergy 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 U d with respect to energy cutoff

The detailed derivations of the RG equations are presented in Supplementary Note 3 of Supplementary Information. 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 longrange Coulomb interactions in such graphene-based superlattice systems 56. 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 Supplementary Note 4 of Supplementary Information). 

When we consider electrons in graphene and substrate on equal footing in Eqs. (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 groundstate 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 Supplementary Note 6 of Supplementary Information. 

Workflow to solve the coupled bilayer Hamiltonian Eqs. ( 5 )

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. 53:

where c 1 = -1. 106103, c 3/2 = 0. 814, c 2 = 0. 113743, c 5/2 = -1. 184994, and 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 54 :

where x = ffiffiffiffi r s p 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 54. 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 57. 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 selfcontained in our model. Equation ( 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. • From the above procedures, we would separately obtain converged HF ground states, Ψ j i d for the substrate, and Ψ j i c for graphene, respectively. From the ground-state wavefunctions Ψ j i d and Ψ j i 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 Supplementary Note 6 of Supplementary Information. 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 Ψ j i d Ψ j i 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 groundstate 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 Supplementary Figure 12. 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 Supplementary Note 6 of Supplementary Information. 

• 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 noninteracting 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 Figs. 2 and 4 of the main text.
section_name: Hartree-Fock approximations assisted by renormalization group approach

section_content: The first principles calculations are performed with the projector augmented-wave method within the density functional theory 58, as implemented in the Vienna ab initio simulation package software 59. 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 60. 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 61, 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 62. 

Our high-throughput filtering of the proper insulating substrate materials for graphene starts from the 2D materials computational database 63. 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.
section_name: Density functional theory calculations

section_content: By designing a dual-gated structure, we used few-layered CrOCl as a 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 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. 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 Supplementary Note 8 of Supplementary Information.
section_name: Experimental measurements of the gaps in graphene-CrOCl heterostructure

section_name: Acknowledgements

section_content: We would like to thank Jian Kang and Jinhai Mao for valuable discussions, and to thank Hanwen Wang for the help in making the plots. X. L., S. Z., Z. G., and J. L. acknowledge support from the National Key R & D program of China (grant No. 2020YFA0309601 ), the National Natural Science Foundation of China (grant No. 12174257 ), and the start-up grant of ShanghaiTech University. Y. W., X. G., K. Y., and Z. H. acknowledge support from the National Key R & D program of China (grant No. 

Article https://doi. org/10. 1038/s41467-023-41293-8Nature Communications | (2023) 14:5550
section_name: Acknowledgements

section_name: Data availability

section_content: The data that support the findings of this study are available at https:// figshare. com/projects/MonoGr-CrOCl/174702.
section_name: Data availability

section_content: The codes that support this study are available from the corresponding author upon request.
section_name: Code availability

section_name: Author contributions

section_content: J. L. conceived the idea, constructed the theoretical model, and supervised the project. X. L., S. Z., Z. G., and J. L. performed theoretical calculations. Y. W., X. G., K. Y., Y. G., Y. Y., and Z. H. performed transport measurements. X. L., Z. G., Z. H., and J. L. analyzed the data. X. L. and J. L. wrote the manuscript with inputs from all authors.
section_name: Author contributions

section_content: The authors declare no competing interests.
section_name: Competing interests

section_name: Additional information

section_content: The online version contains supplementary material available at https://doi. org/10. 1038/s41467-023-41293-8. 

Correspondence and requests for materials should be addressed to Jianpeng Liu. 

Peer review information Nature Communications thanks Jihang Zhu and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available. 

Reprints and permissions information is available at http://www. nature. com/reprints Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 

Open Access This article is licensed under a Creative Commons Attribution 4. 0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/ licenses/by/4. 0/.
section_name: Supplementary information