Spontaneous charge-ordered state in Bernal-stacked bilayer graphene

We propose that a weakly spontaneous charge-ordered insulating state probably exists in Bernal-stacked bilayer graphene which can account for experimentally observed non-monotonic behavior of resistance as a function of the gated field, namely, the gap closes and reopens at a critical gated field. The underlying physics is demonstrated by a simple model on a corresponding lattice that contains the nearest intralayer and interlayer hoppings, electric field, and staggered potential between different sublattices. Combining density functional theory calculations with model analyses, we argue that the interlayer van der Waals interactions cooperating with ripples may be responsible for the staggered potential which induces a charge-ordered insulating state in the absence of the electric field.


I. INTRODUCTION
Bilayer graphene has been extensively studied for more than one decade.While fascinating progresses, such as discovering of the integer quantum Hall effect 1,2 , superconductivity 3,4 , higher-order topological insulators 5 , tunable excitons 6,7 , and topological valley transport 8 , induced by external fields [9][10][11][12][13][14] , doping 15,16 , and twist [17][18][19] , have been reported, the ground state of Bernal-stacked bilayer graphene (BBG) remains controversial.Originally, BBG was thought to be a semimetal with massive Dirac cones at the Fermi level.And a gated field applied perpendicular to the plane breaks the symmetry between the top and bottom layers of Bernal-stacked bilayer graphene (BBG), rendering it an insulator, which has been confirmed by transport 10,12 and photoemission experiments 20,21 .Then, the gap should increase monotonously as a function of the applied gated electric field [22][23][24][25] .However, this is challenged by intriguing experimental observations on ultraclean suspended BBG that resistance varies non-monotonously with the gated electric field at low temperatures 14,[26][27][28] , suggesting the presence of an intrinsic gap that closes and then reopens when an electric field is applied.
So far, much effort has been made to understand the discrepancy.Starting with the intrinsic gap at zero field, various possible candidate states which stem from different origins have been proposed.Using methods like quantum Monte Carlo, functional renormalization group, etc., a layered antiferromagnetic (LAF) state has been suggested as a candidate state in BBG, which is favored by on-site Coulomb interactions [29][30][31][32][33] .By calculating the properties of Landau level n=0 with spin and valley degree of freedom, a canted antiferromagnetic (CAF) state is suggested to be stabilized by isospin anisotropy of electron-electron and electron-phonon interactions 34 .Besides, the quantum spin Hall (QSH) state 33,35 and quantum anomalous Hall (QAH) state 36 are also considered as two potential candidate states with gap opening at zero field, which are favored by spin-orbit coupling and zero-point fluctuations, respectively.Recently, taking short-range interactions into account, a candidate state with the coexistence of nematic and antiferromagnetic states has also been proposed 37 .
However, despite numerous investigations, a definitive explanation for the phenomenon that resistance varies non-monotonously with an electric field remains elusive, which is probably due to the following two reasons.On one hand, most of the previous studies focus on the ground state at zero field, where there are many competing candidate states with very close energies 31 .Consequently, the ground state strongly relies on delicate details of the microscopic model 33,[37][38][39][40][41] that a specific perturbation may favor a particular candidate state as introduced above.Although these studies suggest the presence of a magnetic ground state at zero field 42 , there is no direct experimental evidence, such as a neutron diffraction experiment, for the existence of spontaneous magnetization in BBG.On the other hand, the models employed to investigate the gap include only several tight-binding parameters to describe the low-energy dispersions in the vicinity of the Dirac point 30,31,34 .Consequently, these models fail to capture the realistic behavior of the gap, exhibiting inconsistencies with experimental observations 30,31,34 .
Therefore, it is necessary to investigate the behavior of the gap under an electric field based on a reasonable model to determine the ground state of BBG.Noticeably, some key ingredients such as interlayer van der Waals (VdW) interactions and ripples are often ignored by previous analyses.The interlayer VdW interactions and the ripples which naturally occur in graphene sheets [43][44][45] are crucial to the properties of BBG since interlayer VdW interactions play a dominant role in anchoring the layers at a fixed distance 46 and ripples can drive graphene (including BBG consisting of two layers) into an insulator 47,48 .Furthermore, it has been suggested that a charge-ordered state, which may be favored by the two aforementioned ingredients in BBG, is possible in suspended graphene samples 49 .Thus, taking into account the effect of interlayer VdW interactions and ripples to provide a comprehensive explanation for the field-induced non-monotonic behavior of the gap is an interesting work.
In this paper, we present a novel explanation for the phenomenon observed in BBG that the resistance varies non-monotonically with an applied electric field, corresponding to the gap closes and then reopens with the electric field.We start by demonstrating the underlying physics of the phenomenon using a simple model with staggered potential between inequivalent sublattices on a Bernal-stacked bilayer honeycomb lattice.We reveal that the intrinsic gap at zero field is attributed to the presence of a particular intralayer charge-ordered state which is characterized by an inverted order of the four low-energy bands, where two touched bands shift below the Fermi level while the other two untouched bands move above it [see Fig. 4(b)], in contrast to the disordered case, where touched ones meet at the Fermi level forming a massive Dirac cone [see Fig. 4(a)].As an electric field is applied to this charge-ordered state, the upper band of the two touched bands and the lower band of the two untouched bands move towards and then cross each other, resulting in the non-monotonic behavior of the gap at a small electric field.To validate the proposal that this charge-ordered state exists in BBG, we combine density functional theory (DFT) calculations with model analyses to include the effect of the interlayer van der Waals (VdW) interactions and ripples.We find that interlayer VdW interactions, along with ripples effectively generate a staggered potential between inequivalent sublattices in BBG.Most importantly, the experimental phenomenon is successfully reproduced when the strength of interlayer VdW interactions is on the order of 10 meV.Our results offer a fresh perspective on the field-induced intriguing phenomenon in BBG.
Our paper is structured as follows.In Section II, we provide a comprehensive description of the model and methods employed in our study.Section III presents our primary findings.Specifically, we calculate the gaps and intralayer charge disproportionations as functions of the staggered potential.We also examine how the eigenvalues of the Dirac point and low-energy bands vary when an electric field is applied.Additionally, we investigate the evolution of the gap with respect to the electric field when including the effect of interlayer VdW interactions and ripples.Section IV presents a detailed discussion, and Section V concludes our paper with a concise summary.

II. MODEL AND METHOD
The simple model that we employ to demonstrate the underlying physics for the intriguing phenomenon is given by where Here ( It is easy to find that, ε 1 and ε 2 are contributed by p z orbitals of B 2 and B 1 atoms, respectively, whereas ε ± are contributed by a linear combination of p z orbitals of A 1 and A 2 atoms.To accurately reproduce the intriguing phenomenon observed in BBG, it is imperative to include the effect of previously ignored interlayer VdW interactions and ripples.While ripples can be readily introduced through DFT calculations, the inclusion of interlayer VdW interactions is challenging due to the existence of various corrections such as vdW-DF 50 , TS-vdW 51 , vdW-DF-C 52 , DFT-D 53 , etc. Thus, in order to gain insight into the effect of VdW interactions and eliminate the uncertainty from different choices of corrections, we construct a Hamiltonian as where is the tight-binding Hamiltonian describes the low-energy dispersion of DFT band structures, where the hopping parameters t isjs ′ can be derived through the transformation from Bloch space to maximally localized Wannier functions basis by using WANNIER90 code 54,55 .Four bands close to the Fermi energy which are mainly contributed by p z orbitals of four carbon sublattices are taken into account.Sufficient numbers of long-range hoppings t isjs ′ s are included in order to precisely describe the dispersion of low-energy bands obtained from DFT calculations as shown in Appendix A. The DFT band structures of BBG are calculated by the full-potential linearized augmented plane-wave method 56 within the density functional theory as implemented in the WIEN2K package 57 , where The exchange-correlation interactions are treated by the local-density approximation 58 .A shifted 60 × 60 × 1 special k-point mesh with a modified tetrahedron integration scheme 59 for the sampling of the Brillouin zone is employed.The valence and core states are separated by an energy of 6.0 Ry and the plane-wave cutoff parameter R mt × k max is set to be 7.00, where R mt = 1.33 a.u. is used.The valence wave functions inside the atomic spheres are expanded up to l max = 10 while the charge density is Fourier expanded up to G max = 12.We choose both a charge convergence of 10 −7 e and an energy convergence of 10 −7 Ry as the convergence criteria.A sufficiently large vacuum distance of 17.1 Å is used to eliminate the interactions between periodic images of the layers in the direction perpendicular to the atomic plane.H ef f is the effective Hamiltonian capturing the effect of interlayer VdW interactions, and H E is the aforementioned electric field terms.To obtain the effective Hamiltonian H ef f , we use a VdW potential of interatomic Lennard-Jones type where i(j) and s(s ′ ) correspond to cell and sublattice indices, respectively.V 0 is the strength of interlayer VdW interactions, determining the depth of the potential well.R ss ′ represents the bottom of the sublatticedependent interlayer VdW potential well, including R AA , R BB , R AB , and R BA .Since the interlayer VdW interactions play a dominant role in anchoring the layers at a fixed distance 46 , R ss ′ can be approximated by force equilibrium condition of s sublattice in ith unit cell 13 êss ′ ij = 0, (7)   where êss ′ ij = (r is − r js ′ )/|r is − r js ′ | is the unit vector, and ′ represents the summation over the layer without s sublattice.Thus, for a given V 0 , the potential energy of s sublattices in ith unit cell reads Interestingly, there is a total potential energy difference between A-type and B-type sublattices as namely, the potential energy well of A 1(2) is higher than that of B 1(2) .As a result, electrons redistribute in response to eliminate this difference compared with the case without the interlayer VdW interactions.Thus, interlayer VdW interactions effectively generate a staggered potential between inequivalent sublattices in the layer, namely, H ef f satisfies Thus, the staggered potential is determined as long as the strength of interlayer VdW interactions V 0 is given.

III. RESULTS
Here, we demonstrate firstly the underlying physics for the phenomenon that the gap closes and then reopens with the electric field using the simple model (1) as introduced above.
Starting with the insulating state at zero field (E = 0), we find that the intrinsic gap at zero field is due to the presence of a particular intralayer charge-ordered state where sublattice B is charge-rich while sublattice A is charge-poor.To illustrate this, we calculate the gaps E g and intralayer charge disproportionation (CD) ) as functions of the staggered potential ∆ for different interlayer couplings t ⊥ at zero field, as shown in Fig. 2. As can be seen, the model always predicts a charge-ordered insulating state when ∆ > ∆ c despite the differences in t ⊥ .This insulating ∆ ( e V ) FIG. 2. The gaps Eg and intralayer charge disproportionation δnBA as functions of staggered potential ∆ for three nearest inter-layer coupling with t ⊥ = 0.1 eV, t ⊥ = 0.22 eV, and t ⊥ = 0.3 eV.t0 = 2.7 eV is used here.Gray symbols and left axis describe the band gap while red symbols and right axis depict the charge disproportionation δnBA.
state can be understood using the eigenvalues of Dirac point [Eq.(3)] since they are relevant to the low-energy bands at half filling.We find that when ∆ ≤ t ⊥ , the system remains metallic as the Fermi level lies between the degenerate eigenvalues ε 1 and ε 2 .In contrast, when ∆ > t ⊥ , ε − and ε 1 (ε 2 ) are inverted compared with the disordered case [ see Fig. 1(b) and 3(c) ], resulting in a band-inverted charge-ordered insulating state with a gap of E g = ∆ − t ⊥ .The latter may be relevant to the zero-field insulating state observed in BBG.For example, t ⊥ ≈ 0.22eV 24 in BBG, a weakly intralayer CD with critical δn BA ≈ 0.07 can make it an insulator (Fig. 2).Since a reduced threefold symmetry characteristic is observed by high-resolution scanning tunneling microscopy 43 , implying an intralayer CD, we argue that the insulating state of BBG at zero field is probably due to the presence of this band-inverted intralayer charge-ordered state.
Proceeding to analyze the effect of an electric field E on the band-inverted charge-ordered insulating state, we find that the gap decreases for E < E c , increases for E > E c , with the gap closing at critical value E = E c , which is reminiscent of the intriguing phenomenon observed experimentally that resistance varies non-monotonically with an electric field 26 .A downward electric field can drive the electrons from the bottom to the top layers, resulting in interlayer CD. Figure 3(a) shows δn 21 and δn BA (E) − δn BA (0) as functions of the electric field, where δn 21 = n A2 +n B2 -n A1 -n B1 is the order parameter of the interlayer CD while δn BA (E) and δn BA (0) ≈ 5.84512 × 10 −3 are the order parameters of the field-dependent and zero-field intralayer CDs, respectively.Noticeably, as the band gap is relatively small, we have chosen t ⊥ and ∆ comparable to the bandgap to clearly demonstrate the detailed evolution of all four relevant low-energy bands with the small electric field, which does not alter the underlying physics.As can be seen, distinct behaviors of δn 21 and δn BA (E) − δn BA (0) in E < E c and E > E c imply a phase transition from phase I to II at E = E c .As a small electric field mainly affects the low-energy bands of the system, the phase transition can also be understood using the eigenvalues of the Dirac point.Figure 3(b) shows the eigenvalues of the Dirac point varying with the electric field.We find ε 2 and ε + raise up, whereas ε 1 and ε − go down with the increase of the electric field.This is because ε 2 , ε + , ε 1 , and ε − are mainly contributed by p z orbitals of B 1 , A 1 , B 2 , and A 2 sites, respectively, where the on-site potentials of B 1 and A 1 sites increase, whereas those of B 2 and A 2 sites decrease when an electric field is applied.Given that ε − is higher than ε 2 at E = 0, ε − and ε 2 will move towards and then cross each other as shown in Fig. ⊥ )/(ed∆).Thus, the gap of the bandinverted intralayer charge-ordered insulating state varies non-monotonously with the electric field.In brief, the findings above can be summarized as a schematic illustration in Fig. 4, demonstrating the following physics.(i) An intralayer charge-ordered state can open an intrinsic gap in the BBG lattice at zero field when band inversion between the two touched bands and the lower band of the two untouched bands occurs, compared with the disordered case [see Fig. 4(b) and 4(a)].(ii) When an electric field is applied to this charge-ordered insulating state, the upper band of the two touched bands and the lower band of the two untouched bands move towards and then cross each other as shown from Fig. 4(b) to 4(c), resulting in a nonmonotonic behavior of the gap that closes and then reopens with the electric field.As we proposed in Section II that the interlayer VdW interactions can effectively generate a staggered potential between intralayer inequivalent sublattices.Besides, the ripples which naturally exist in BBG may favor a charge-ordered state.Thus, this band-inverted intralayer charge-ordered insulating state may be the key to the non-monotonic resistance phenomenon observed in BBG, and it is imperative to study the effect of previously ignored interlayer VdW interactions and ripples.
To validate our proposal that the intriguing phenomenon observed in BBG is due to the presence of the aforementioned band-inverted charge-ordered state, the model presented in Section II, specifically Eq.( 4), which includes the effect of interlayer VdW interactions and ripples, is employed to calculate the gap of BBG.For a given strength of interlayer VdW interactions V 0 , the effective Hamiltonian H ef f , especially the staggered potential ∆, should be self-consistently determined.Thus, we derive firstly the potential energy difference between A and B sublattices using the method presented in Section II.We find δU AB ≈ 1.08N V 0 and δU AB ≈ 2.08N V 0 for flat structure and Peierls-type rippled structure with ∆d = 0.01nm, respectively, where N is the number of unit cells.Using Eq.( 10), the effective staggered potential is calculated as ∆ = δU AB /(N δn BA ).Thus, δn BA and ∆ are self-consistently determined by combining this equation with Eq.( 4) as long as V 0 is given, indicating that H ef f can be determined self-consistently.Therefore, the gaps and CDs of bilayer graphene for a given strength of interlayer VdW interactions are calculated.Figure 5 shows the gaps of BBG as functions of the electric field for flat and rippled structures with the strength of interlayer VdW interactions of V 0 = 30.18meV and V 0 = 20.22 meV, respectively, where the gap decreases for E < E c , increases for E > E c , with the closure of the gap at E c = 20 mV/nm.As the resistance R ∝ exp[E g /(k B T )], our results are in excellent agreement with the experimental phenomenon that the resistance decreases and then increases with the electric field, where minimal resistance is at critical electric field E c ≈ 20 mV/nm 26 and the zero-field gap is E g ≈ 2 meV 27,28,60,61 .It is necessary to mention that a weakly spontaneous chargeordered state occurs for both cases with intralayer CD δn BA = 0.10 ∼ 0.12, where δn BA changes very little with the electric field while δn 21 increases from 0 to 3 × 10 −4 .In addition, a smaller V 0 can lead to the same behavior of gap for the rippled structure compared with the flat case, suggesting that the ripple concerned and interlayer VdW interactions are cooperative.Thus, interlayer VdW interactions cooperating with ripples can effectively generate a staggered potential between inequivalent sublattices, which induces an intralayer charge-ordered insulating state, resulting in the experimental phenomenon observed in BBG.The strength of the interlayer VdW interactions is of the order of 10 meV.Please note that R ss ′ s are not tunable parameters which are determined by the force equilibrium condition.The critical values of V 0 that causes the metal-to-insulator phase transition for flat and ripple structures are 29.85 meV and 20.00 meV, respectively.
IV. DISCUSSION Here, a simple model has been used to demonstrate the underlying physics for the intriguing phenomenon observed in BBG that the resistance decreases and then increases with the electric field at low temperatures.We ascribe this phenomenon to the presence of a band-inverted charge-ordered insulating state in BBG.Our proposal is further confirmed by combining DFT calculations with model calculations, where we take into account the effect of the interlayer VdW interactions and ripples.Our results are reliable because our calculations include not only the effects of the nonlocal Coulomb interactions and remote hoppings at the DFT level but also the previously ignored ingredients, namely interlayer VdW interactions and ripple in the layer.Noticeably, although the interlayer Coulomb interactions are stronger than the interlayer VdW interactions, the potential difference between inequivalent sublattices generated by the interlayer Coulomb interactions is much smaller than that generated by the interlayer VdW interactions as presented in Fig. 6(a).Thus, the interlayer VdW interactions play a major role in determining the intralayer charge-ordered state of BBG.However, this key ingredient is often ignored by previous studies.
The validity of modeling BBG by a bilayer honeycomb lattice with staggered potentials is as follows: Although DFT can provide a reasonable energy difference between AA-stacking bilayer graphene and BBG due to cancellation of the uncertainty from interlayer interactions simultaneously 62 , it fails to capture the interlayer interactions for each structure individually.Therefore, it is necessary to model the uncertainty arising from interlayer interactions, which may induce a staggered potential between A 1(2) and B 1(2) sublattices due to their inequivalent interlayer environment.Indeed, our results suggest the presence of a staggered potential due to the interlayer interactions, e.g., VdW interactions.The calculated small band gap, which is consistent with the experimental observations in bilayer graphene devices 27,28,60,61 , strongly indicates that the system is in the critical region of insulator-to-metal transition.
Although varieties of corrections such as vdW-DF 50 , TS-vdW 51 , vdW-DF-C 52 , DFT-D 53 and so on have been proposed, the electronic properties of the BBG and graphite obtained from these corrections are diverse from each other 63,64 .Thus, in order to gain insight into the effect of VdW interactions and eliminate the uncertainty from different choices of corrections, it is necessary to treat the interlayer VdW interactions as free parameters, namely modeling the effect of interlayer VdW interactions.In our calculations, the model Hamiltonian for interlayer VdW interactions is obtained naturally from the interatomic Lennard-Jones potential, where the strength of interlayer VdW interactions V 0 serves as the sole free parameter except for the electric field in Eq.( 4).Noticeably, the isotropic nature of the Lennard-Jones VdW potential can not capture the anisotropic properties of BBG.Thus, we also employ the Kolmogorov-Crespi potential 62 to take into account the interlayer interactions with anisotropy as introduced in Appendix B. Similar behavior that the band gap varies non-monotonically with the electric field has also been observed as shown in Fig. 6(b).
Although several correlated symmetry-broken gapped states with parabolic dispersion relation have been proposed at zero field, namely LAF, CAF, QAH, and QSH, the low-energy bands of BBG observed experimentally at a small applied magnetic field cannot be explained within the framework of parabolic bands, which predicts roughly equidistant Landau levels at low temperatures 65 .Besides, there is a pronounced asymmetry in the cyclotron mass between hole-and electron-doping 24 .These findings raise doubts about the candidates which exhibit a parabolic dispersion relation with particle-hole symmetry near the Fermi level.Moreover, as the temperature increases from zero, two resistance transitions occur.One transition is observed at ∼ 12K 66 , while the other occurs at 200 ∼ 250K 67 , which is suggested to be caused by the interlayer ripple scattering effect.As the charge-ordered state we proposed still exists even when the gap is closed, it may suggest that the former transition corresponds to the evolution from the charge-ordered insulating state to the charge-ordered metallic state, while the latter transition corresponds to the change from the charge-ordered metallic state to the disordered state.
Ripples are inherent features of BBG, arising from the natural undulations of suspended graphene sheets.It has been proposed that suspended graphene sheets are not perfectly flat showing ripples with an amplitude of about 1 nm [43][44][45] with dislocations 68 .Here, for simplicity, we take the Peierls-type ripple into account, which is energetically favored by elastic effects 69 .However, it should be noted that the ripples in BBG exhibit a complex nature.Thus, it is interesting to study how the experimentally observed ripples cooperate with interlayer VdW interactions to affect the properties of the intralayer charge-ordered state.
Although the charge-ordered state we studied has been previously investigated 70 , the properties of this chargeordered state under an external field have not been explored before.Noticeably, a low-energy theory based on a 2 × 2 Hamiltonian matrix with consideration of the charge-ordered characteristic is used to study the properties of BBG 36 .However, it fails to deal with the properties of the charge-ordered state concerned here, where the low-energy bands are inverted.Besides, although lowenergy theory based on a 4 × 4 Hamiltonian matrix is also proposed, it does not take into account the effect of a charge-ordered state 23,71 .
Here, we study the phenomenon observed in BBG that the gap varies non-monotonously with the electric field.Our results strongly suggest that the ground state of BBG is a charge-ordered insulating state.Therefore, further experiments are needed to confirm the ground state of BBG.There are two experimental approaches to identify this state: one is angle-resolved photoelectron spectroscopy and the other is scanning tunneling spectroscopy.An experiment based on angular-resolved photoemission spectroscopy should be done at low temperatures, without external perturbations, to detect the low-energy bands of BBG.If the low-energy bands are inverted, the ground state of BBG is a charge-ordered state.Alternatively, the scanning tunneling spectroscopy would be sensitive to the charge ordering at atomic scale, allowing one to measure spatial variations of the local density of states to determine the ground state of BBG.

V. CONCLUSION
In conclusion, an intriguing phenomenon that the resistance varies non-monotonously with an electric field applied perpendicular to the plane has been observed at low temperatures in BBG.Here, we suggest that this phenomenon is probably due to the presence of a spontaneous charge-ordered insulating state in BBG.The underlying physics is illustrated by a simple model on BBG lattice with staggered potential between inequivalent sublattices.To validate our proposal, we combine DFT calculations with model calculations to include the effect of the interlayer VdW interactions and ripples.We find that the interlayer VdW interactions cooper-ating with ripples can effectively generate a staggered potential in BBG.Remarkably, we have successfully reproduced the gap amplitude and the critical electric field when the strength of the interlayer VdW interactions is on the order of 10 meV.Our results provide a new perspective on the non-monotonic resistance phenomenon in BBG and suggest that the ground state of BBG is likely a charge-ordered state.We argue that angular-resolved photoemission spectroscopy studies at zero field or scanning tunneling spectroscopy can be used to identify the ground state at low temperatures.Due to the fact that only including t 0 and t ⊥ is not sufficient to well describe the low-energy DFT bands of BBG, we establish a tight binding model H D with longrange hoppings to accurately describe the entire dispersions of low-energy DFT bands.For both flat and ripple structures, the fitted bands provided by the tightbinding Hamiltonian of all p z orbitals (blue) and DFT bands (red) match well as shown in Fig. 7.

ACKNOWLEDGEMENT
Appendix B: The total potential energy difference between A-type and B-type sublattices derived from the Kolmogorov-Crespi potential It has been pointed out that the isotropic nature of the Lennard-Jones VdW potential can not capture the anisotropic properties of the graphitic systems.Then, Kolmogorov and Crespi take into account the in-plane and out-of-plane anisotropy, proposing the socalled Kolmogorov-Crespi potential to describe the interlayer interactions in graphite systems 62 .The interatomic

FIG. 1 .
FIG. 1.(a) The structure of BBG, where the nearest intralayer and interlayer hoppings, namely, t0 and t ⊥ are presented.(b) The low-energy bands of the disordered case, where H∆ = 0 and HE = 0.
3(c) to 3(f) [or in Fig. 4(b) and 4(c)].Consequently, the gap decreases for E < E c and increases for E > E c .Notably, the gap closes at the critical electric field E c = (∆ 2 − t 2

FIG. 4 .
FIG. 4. The charge distribution and bands in the vicinity of the Dirac point for different cases with (a) absence of electric field and staggered potential, (b) staggered potential is larger than the critical value but without an electric field, (c) both electric field and staggered potential are larger than the critical value.The size of the ball represents the charge population in corresponding site.

FIG. 5 .
FIG.5.The gaps of flat and rippled structures as functions of the electric field.d = 0.34 nm and a = 0.142 nm are used.∆d = 0.01 nm in rippled structure which is presented in inset.

FIG. 6 .
FIG. 6.(a) The potential energy differences between inequivalent sublattices generated individually by the interlayer Coulomb interactions (δU N C BA ) and that generated individually by the interlayer VdW interactions (δUBA) for the flat structure as functions of interlayer distance d, where V0 is the strength of the interlayer VdW interactions and V ⊥ is the strength of the nearest neighbor interlayer Coulomb interaction.V (r) = V ⊥ d r is used for the interlayer Coulomb interactions.N is the number of unit cells.(b) The calculated band gap of model (4) as functions of the electric field, where the Kolmogorov-Crespi potential is used to derive H ef f and α is the correction factor for the repulsive term of the Kolmogorov-Crespi potential in BBG as introduced in Appendix B.
This work is supported by National Natural Science Foundation of China (Grants No. 12274324, No. 12004283) and Shanghai Science and Technology Commission (Grant No. 21JC405700).

FIG. 7 .
FIG. 7. The bands obtained from DFT calculations and the corresponding wannier fittings for flat structure (a) and for ripple structure with ∆d = 0.01nm (b).
, H K , H τ , H ∆ , and H E denote the nearest intralayer