Introduction

Strain in graphene can be served as an effective tuning parameter in tailoring its exotic properties reflected by different degrees of freedom such as spin and valley1,2,3,4,5,6. In particular, it has been predicted that a non-uniform strain distribution imposed on graphene arising from either external stress1,2 or interfacial coupling3,4 can induce a strong gauge potential that in turn can act as a pseudo-magnetic field (PMF) on the Dirac electrons. Indeed, it has been demonstrated experimentally that locally strained graphene in the form of nanobubbles can induce PMFs, as manifested by the emergent pseudo-Landau levels observed using the scanning tunneling microscopy/spectroscopy (STM/S)5,6. More recently, existence of large-area, spatially distributed PMFs has also been reported for single-layer graphene on a black phosphorus substrate, which provides strain textures due to the mismatched symmetries and lattice constants of the heterostructure7. Separately, intensive theoretical and experimental studies have demonstrated that low-angle twisted bilayer graphene (TBG) can exhibit a variety of exotic phenomena, including superconductivity and quantum phase transitions8,9,10,11,12,13,14,15,16,17,18,19. In addition, it has been reported that a small uniaxial heterostrain in low-angle TBG can suppress Dirac cones and lead to the emergence of a zero energy resonance20. In practice, low-angle TBG itself possesses large-scale, natural, and periodic strains due to the incommensurate moiré superstructures3,4,21. Such strain has been theoretically studied to significantly affect the electronic properties of the large-scale moiré pattern22,23. However, yet to date, experimental explorations of such natural and intrinsic strain effects arising from lattice relaxations in the strongly coupled homostructural system, especially the pseudo-magnetic behaviors, are still lacking22,24,25.

In this letter, by combining STM/S experiments and large-scale tight-binding calculations, we study the electronic properties of TBG with extremely small twist angles (\(\theta\, <\, {1}^{\circ }\)). We directly detect pseudo-Landau levels generated by the intrinsic PMFs, which originate from the interplay between interlayer interactions and in-plane strain fields. The PMFs form a ring-structured vortex lattice associated with the moiré pattern, and are periodically distributed over the whole sample. The vortex lattice, in which all angular incommensurability is focused at specific singularity points, is chiral for the charge carriers at different valleys in k-space, and can be further tuned, together with the PMFs, by the twist angle and uniaxial heterostrain in the bilayer system. Indeed, the existence of periodic PMFs in such hexagon-structured system provides a possible platform to further explore quantum anomalous Hall effects if a tiny external magnetic field is applied to break the time-reversal symmetry26.

Results

Pseudo-magnetic fields with twist angle \(\theta =0.4{8}^{\circ }\)

Our samples were prepared by chemical vapor deposition on Cu–Ni alloyed substrates. Cu–Ni alloyed substrates have been widely used to grow bilayer graphene27. High quality bilayer graphene samples with multiple domains were obtained after appropriate annealing processes (Supplementary Note 1). Figure 1a shows an STM topography image with hexagonal symmetry moiré patterns. The atomic structure of the moiré pattern can be determined explicitly by the topography image20,28,29. Here, we follow the same method to identify the twist angle in Fig. 1a, which gives \(\theta =0.4{8}^{\circ }\) (Supplementary Note 5).

Fig. 1: Electronic properties of TBG with twist angle \({\boldsymbol{\theta }}=0.4{8}^{\circ }\).
figure 1

a STM topography image (\(100\ {\rm{nm}}\times 100\ {\rm{nm}}\)) with hexagonal moiré pattern (top panel) and the atomic model of AA, AB, and BA stacking regions (bottom panel). The three moiré wavelengths are: \(L1\approx L2\approx L3\approx 29.6\ {\rm{nm}}\). Sample bias \(V=100\ {\rm{mV}}\), tunneling current \({I}_{t}=1.0\ {\rm{nA}}\). b Logarithmic dI/dV spectra measured at AA and AB regions. The two solid (dashed) lines were taken at different AA (AB) regions to show reproducibility (curves are vertically shifted for clarity). c Calculated local density of states in the AA and AB regions for deformed (upper panel) and rigidly twisted (bottom panel) cases. d, e Calculated local potential \(V\) and pseudo-magnetic fields \({\bf{B}}\), respectively.

Full size image

Next, we characterize the electronic properties of the TBG samples by the dI/dV tunneling spectra. Generally, dI/dV conductance is proportional to the local density of states in the vicinity of the STM tip position. To rule out the possible tip effects on the measured dI/dV conductance, we have performed dI/dV measurements on TBG with large twist angles. The obtained spectra show typical angle-dependent van Hove singularities which are highly consistent with previous reports30 (Supplementary Note 2). Figure 1b shows dI/dV spectra taken in both AA and AB stacking regions. The AA region spectra show a series of pronounced peaks (indicated by numbers) with nearly equal energy spacing (Supplementary Note 3). Surprisingly, these peaks are typically different from those of TBG in the absence of external magnetic fields10,20,30,31, but rather similar to the Landau levels of a two dimensional electron gas in graphene under strong external magnetic fields32,33. In stark contrast, these resonances are hardly found in the AB region. For extremely low-angle TBG, the charge density in the AB region resembles that of ideal Bernal-stacking bilayer graphene4.

To interpret our experimental observations, we have performed theoretical calculations of the superstructure shown in Fig. 1a. The as shown moiré pattern unit cell contains 57964 atoms that is too large even for state-of-the-art first-principles methods. Therefore, we adopted a widely used full tight-binding model of bilayer graphene with an angle-dependent Hamiltonian9,20. Figure 1c shows the calculated local density of states by utilizing the Lanczos recursion method in real space34. Considering the lattice deformation of the bilayers (upper panel), the result reproduces the energy peaks that appear in the measured spectra (Fig. 1b) whereas the appearance of high energy resonances are dramatically different in the rigidly twisted one (bottom panel). The implementation of the lattice distortion in a rigidly twisted bilayer graphene will be discussed later. As a side issue, our calculations (Fig. 1c) also reproduce the observed resonance peak at the Fermi energy. Recent investigations of such system have provided experimental evidence of such strong electronic correlations near the magic angle TBG by means of STM/S16,17,18,19.

For bilayer graphene under external magnetic fields, the massive Dirac electrons are quantized with energies \({E}_{N} =\pm\! \hslash {\omega }_{{\rm{c}}}\sqrt{N(N-1)}\), \({\omega }_{c}=eB/{m}^{* }\), N = 0, 1, 2,… , where \(\hslash\) is Planck’s constant and \({m}^{* }\) is the effective mass of electrons24,33. Here, the resonant peak positions are linear with \(\sqrt{N(N-1)}\), indicating the feature of pseudo-Landau levels in bilayer graphene. For the TBG with twist angle of \(0.4{8}^{\circ }\), the fitted PMF value is about 9 T, as shown in Fig. 2a. In addition, the PMF in low-angle TBG does not break the time-reversal symmetry and we further calculated the pseudo-Landau levels under external magnetic fields (Fig. 2b). An interesting observation is that, under the strong external magnetic fields, the pseudo-Landau levels are splitted into two peaks with an energy gap of about 15 meV due to the breaking of valley degeneracy. Similar phenomena have been reported in earlier studies of strained monolayer graphene35 except that in the present study, such valley splitting is not as large as theoretically expected. The suppression can be understood, because the theoretically expected valley splitting is valid only for well-defined K and K′ valleys, where the PMF is uniform. While in our systems the observed PMF has drastically varying ring-like structures in real space. As a result, the two valleys are not well defined; therefore, the theoretically expected valley splitting value is no longer applicable. Instead, the valley splitting observed in this study is extrapolated from the regions where the pseudo-magnetic field only varies mildly in real space.

Fig. 2: Pseudo-Landau levels in the deformed twisted bilayer graphene with \({\boldsymbol{\theta }}=0.4{8}^{\circ }\).
figure 2

a Linear fit of the equation \({E}_{N}\propto \sqrt{N(N-1)}\) and the obtained pseudo magnetic fields is about 9 T. b Calculated LDOS at AA region under the external magnetic fields, in which we can confirm the splittings of the pseudo-Landau level due to the break of the valley degeneracy.

Full size image

Theoretical calculations

To illustrate the delicate physics involved in reaching the agreement between the theory and experiment shown in Fig. 1, we note that, in graphene-based van der Waals structures, compelling evidences have revealed that the lattice relaxation effects are crucial in determining their electronic properties3,4,15,29. Especially, in low-angle TBG, the superstructure undergoes lattice distortions and structural transformations, which can be described as a lattice deformation of the moiré superlattice due to the interplay between the interlayer interaction and the in-plane strain fields. That is, on one hand, the in-plane forces move atoms to maximize the area of AB/BA stacking domains, which has the minimum binding energy. On the other hand, the strain induced by the atomic displacements hinders such in-plane atomic rearrangement3,4. By fitting the atomic structures obtained from molecular dynamics simulations of TBG in refs. 3,4, the in-plane displacement \(\Delta d\) and out-of-plane displacement \(\Delta z\) of individual atoms can be approximately expressed as:

$$\Delta d(r)= \,\Delta D\cdot \sin\left(\frac{\pi r}{2{l}_{D}}\right)[1-\Theta (r-{l}_{D})]\\ +\, \Delta D\cdot {{\rm{e}}}^{[-{(r-{l}_{D})}^{2}/{\sigma }_{D}]}\Theta (r-{l}_{D}),\\ \Delta z(r)= \,\pm \Delta Z\cdot {{\rm{e}}}^{(-{r}^{2}/{\sigma }_{Z})}[1-\Theta (r-{l}_{Z})],$$
(1)

where \(r\) is the distance between the individual atoms and the nearest AA point in the x–y plane, \(\Theta\) is the Heaviside step function, \(\Delta D\) is the maximum in-plane displacement and \(\Delta Z\) is the maximum out-of-plane displacement that are dependent on the rotation angle4. lD, σD, lZ, σZ are constants.

The deformed samples are constructed using the following procedure: First, the rigidly twisted bilayer graphene is constructed according to the commensurability (Supplementary Note 5) and the position of each atom is well known. Next, we modify the in-plane and out-of-plane positions of atoms according to the position-dependent expression in Eq. (1) (Supplementary Note 6). For instance, for the in-plane deformation, an individual atom moves \(\Delta d\) towards the nearest AA point. In the out-of-plane deformation, each individual atom has a displacement in the \(z\) direction with \(+\) sign for the top layer and − sign for the bottom layer. The maps of \(| \Delta d|\) and \(| \Delta z|\) for TBG with \(\theta =0.4{8}^{\circ }\) are plotted in Fig. 3c, d, respectively. The movement of the atoms have the same tendency as that obtained from molecular dynamics simulations in reference4. In this part, we assume that the deformed sample keeps the period of the rigid TBG. In the rigid sample, we use the AA point as the rotation center. The moiré supperlattice has point group \({D}_{3}\) generated by a three-fold rotation about the axis formed by the AA point (\({C}_{3}\)) and a two-fold rotation about the axis perpendicular to the AA point (\(C^{\prime}_{2}\)). The deformed sample maintains the \({D}_{3}\) point group symmetry36,37. For another sample with uniaxial heterostrain, the supperlattice for both rigidly twisted and deformed cases lose the \({D}_{3}\) symmetry. In Eq. (1), the parameters \(\Delta D\), \({l}_{D}\), \({\sigma }_{D}\), \(\Delta Z\), \({l}_{Z}\), and \({\sigma }_{Z}\) determine the profile of the deformed structure and are fitting parameters to generate proper structures with properties matching the experimental observations. In fact, in our samples, the substrate suppresses the lattice distortion, and the movements of the atoms are weaker than those in a free-standing twisted bilayer graphene38. Moreover, the in-plane deformation has a dominant effect on the electronic properties, which can be seen from the LDOS in Fig. 3e. In the deformed sample in Fig. 1c, the fitting parameters are the sets B (in-plane) and E (out-of-plane) labeled in Fig. 3.

Fig. 3: Calculated results of the distorted TBG structure with \({\boldsymbol{\theta }}=0.4{8}^{\circ }\).
figure 3

a Schematic model of the moiré pattern. b Absolute magnitude of different in-plane atomic displacements (top panel) and out-of-plane displacements for the deformed system along the path MNP defined in a. For the in-plane displacements, the parameters (\(\Delta D\), \({l}_{D}\), \({\sigma }_{D}\)) are: A (0.04 nm, 5.91 nm, 1.14), B (0.07 nm, 5.91 nm, 1.14), C (0.07 nm, 3.70 nm, 1.14), D (0.07 nm, 3.70 nm, 11.44). For the out-of-plane displacements, the parameters (\(\Delta Z\), \({l}_{Z}\), \({\sigma }_{Z}\)) are: E (0.0115 nm, 4.93 nm, 0.7), F (0.023 nm, 4.93 nm, 0.7), G (0.0058 nm, 4.93 nm, 0.7). c, d Maps of the absolute magnitude of the in-plane atomic displacement \(| \Delta d|\) and out-of-plane displacement \(\Delta z\) upon deformation of rigidly TBG calculated from Eq. (1), respectively. e, f Local density of states in the AA region of deformed systems with different in-plane and out-of-plane displacements, respectively.

Full size image

The structure deformation leads to an electrostatic potential \(V\), proportional to the local compression/dilatation, and a pseudo-vector potential \({\bf{A}}\), associated with the shear deformation at one site. They can be written as2,3:

$$V= \,2{g}_{1}\frac{({d}_{1}+{d}_{2}+{d}_{3})/3-d}{d},\\ {A}_{x}= \, c\frac{\beta {\gamma }_{0}}{d}({u}_{xx}-{u}_{yy}),\\ {A}_{y}= \, -c\frac{2\beta {\gamma }_{0}}{d}{u}_{xy},$$
(2)

where \({g}_{1}\approx 4\ {\rm{eV}}\) is the deformation potential for graphene, \({d}_{1}\), \({d}_{2}\), and \({d}_{3}\) are the first neighbor inter-atomic distances in the deformed lattice, d = 1.42 Å is the equilibrium carbon–carbon distance, \(c\approx 1\) is a numerical factor depending on the detailed model of chemical bonding, \(\beta \approx 2\) is the electron Grüneisen parameter, \({\gamma }_{0}=3.2\ {\rm{eV}}\) is the strength of first-neighbor interaction in the plane, \({u}_{ij}\) is the deformation tensor. The PMF is given by \({\bf{B}}=\nabla \times {\bf{A}}\), which can be estimated as3:

$$| {\bf{B}}| \approx \frac{\hslash c}{e}\left(\frac{2\beta }{3}\right)\frac{\overline{u}}{{a}_{m}d},$$
(3)

where \({a}_{m}\) is the period of the moiré pattern. And \(\overline{u}\) is the magnitude of the shear deformation, which can be calculated from the first neighbor interatomic distances in the deformed sample as:

$$\overline{u}=\frac{\sqrt{3{({d}_{3}-{d}_{2})}^{2}+{({d}_{2}+{d}_{3}-2{d}_{1})}^{2}}}{2d}$$
(4)

The local deformation potential \(V\) and the magnitude of the pseudo-magnetic fields \({\bf{B}}\) are plotted in Fig. 1d, e, respectively. Such structural deformation results in a non-uniform \({\bf{B}}\) with a maximum value of about 6 T. One of the effects of applying a strong external magnetic field in a two-dimensional electronic structure is the Landau quantization of the eigenstates and, subsequently, the quantum Hall effect in the transport properties. In our samples, the PMFs are not uniform, but form a ring-structured vortex lattice around the centers of the AA regions, similar to an Abrikosov vortex in superconductors3. Here, these pseudo-Landau levels are strongly correlated to the space-dependent PMFs, similar to the Landau levels that appear in the presence of a real magnetic field. It is important to emphasize that the PMFs in low-angle TBG intrinsically arise from the lattice deformation. We expect that, as the displacements of the atoms are more significant in a structure with smaller twist angle3, the induced PMFs will also be stronger.

Twist angle and heterostrain effetcs

The pseudo-Landau level behavior with twisted angle of \(0.6{5}^{\circ }\) has also been studied by STM/S, in which the observed pseudo-magnetic fields are fitted to be 8 T (Supplementary Note 7). The calculated PMFs are in qualitative agreement with the fitted one from the \({E}_{N}\) formula. The fitted PMFs decrease more smoothly in the space comparing to the calculated values, and this difference can be explained in the following: (i) The expression of the Landau energies \({E}_{N}\) is valid for the case of bilayer graphene under a uniform external magnetic field. On the contrary, in our sample, the deformation-induced PMFs are non-uniform. (ii) As can be seen from the Landau energies \({E}_{N}\), the fitted PMFs can be modulated by the effective mass \({m}^{* }\), of which the value is taken from other experimental results24,33, and the effective mass can also be space-dependent due to the local stacking structure. All in all, in this paper, both the theoretical and experimental values of maximum PMFs increase when the twist angle decreases, which is consistent with the fact that the strength of the lattice distortion becomes stronger for a TBG with smaller rotation angles.

The results presented so far have demonstrated that the twist angle can be an effective tuning parameter in tailoring the lattice deformation and subsequently tuning the electronic and magnetic properties in low dimensional systems. Next, we investigate heterostrain effects on the PMFs and their distributions in the low-angle TBG entity. The heterostrain is uniaxial, which originates from the pinning of the graphene layer at its grain boundaries or step edges during the epitaxial growth (Supplementary Note 4 and Note 9)20. Figure 4a shows the topography image with severely distorted moiré pattern. Such distorted moiré structures are due to the existence of uniaxial heterostrain, where the moiré pattern can serve as a magnifying glass20,28,29,39. Namely, when the honeycomb lattice is slightly strained, the resulting moiré pattern is significantly distorted, in which the twist angle and strain can be determined by the STM images. For the sample in Fig. 4a, we obtain the twist angle and heterostrain to be \(\theta =0.9{8}^{\circ }\) and \(\sigma =0.78 \%\), respectively (Supplementary Note 5).

Fig. 4: Electronic properties and PMF behaviors of TBG with uniaxial heterostrain.
figure 4

a STM topography image (\(65\ {\rm{nm}}\times 50\ {\rm{nm}}\)) with distorted moiré pattern. Twist angle \(\theta =0.9{8}^{\circ }\) with heterostrain \(\sigma =0.78 \%\). Sample bias \(V=200\ {\rm{mV}}\), tunneling current \({I}_{t}=1.0\ {\rm{nA}}\). The three moiré wavelengths are: \(L1=20.3\ {\rm{nm}}\), \(L2=12.3\ {\rm{nm}}\), and \(L3=19.2\ {\rm{nm}}\). b Linear fittings of the peak energies and Landau index obtained from both the measured and calculated datas. The estimated PMF value is about 6 T. c dI/dV spectra with interval of 1 nm taken from AA to AB region indicated by black dots in a. The curves have been vertically shifted for clarity. d Local density of states computed from the same locations in a. The obtained features are in good agreement with our experimental results. e, f Calculations of the magnitude of the pseudo-manetic fields \({\bf{B}}\) produced by the atomic distortion of the top layer (unstrained one, e) and the bottom layer (strained one, f), respectively.

Full size image

 Figure 4c, d shows dI/dV spectra and calculated LDOS curves taken along a line from an AA to an AB region. Similar to the findings in the \(\theta =0.4{8}^{\circ }\) sample, pseudo-Landau levels are also observed as a series of pronounced peaks. The energy spacings between the pseudo-Landau levels vary from the AA to AB region, which indicates the non-uniform characteristics of the PMFs. Significantly, the two low-energy pseudo-Landau levels near the Fermi level are similar to the correlation-induced splitting of flat bands reported recently in STM/S studies of the magic-angle TBG16,17,18. Such enhanced electronic correlation effects can be understood within the context of heterostrain in low-angle TBG (Supplementary Note 10)20. Figure 4b shows the corresponding linear fitting of pseudo-Landau levels, which gives a PMF value of 6 T. By linear extrapolating the LDOS at high energy in Fig. 4c, we deduced the center of the states, which changes with the positions. For our sample with a considerable heterostrain of \(\sigma =0.78 \%\), it is still hard to visualize the atomic displacements from the high resolution STM image. Despite this, the topography image of large area moiré patterns shows distorted six-fold symmetry, which is in good agreement with our atomic simulations (Supplementary Fig. 5).

In accordance with recent reports20,25, high energy peaks appear in AB regions, as shown in Fig. 4c, d. These peaks, associated with a partial band gap opening in a high energy moiré band20,40 or localized domain wall modes arising from strongly confined moiré potentials25, suffer minor changes in the presence of heterostrain and atomic deformations (Supplementary Note 8 for the effects of heterostrain and deformation on the electronic properties of TBG). This is expected since these lattice deformations and pseudo-magnetic fields only occur around AA regions, as shown in Fig. 4e, f. Due to the considerable heterostrain applied to the bottom layer here, the shape and amplitude of the arising PMFs in each layer are no longer identical. Obviously, the strained layer has been more extended than the unstrained one, resulting in larger PMFs, which go up to 7 T. This tendency is in good agreement with previous findings3,4. Besides, the PMFs are stronger along the direction of the uniaxial heterostrain within each circular regions. As a consequence of uniaxial strain, the local deformation potential \(V\) on the strained layer is 0.035 eV higher than that of unstrained one (Supplementary Note 9). Thereby, both the rotation angle and uniaxial heterostrain can be utilized to fine-tune the geometry, amplitude, and periodicity of the PMFs. In Fig. 5, we plot the twisted angle dependent PMFs and their spatial distributions around the AA regions. The locations in experiments are equally taken with interval of 1 nm or 2 nm. The experimental fitted results show that the PMFs increase with the decreasing twist angles and the PMF area occurs near the AA regions which reaches its maximum at the AA/AB transitions. The experiments and calculations show good agreements.

Fig. 5: The fitted pseudo-magnetic field of TBGs with different twisted angles around the region of AA/AB transtion.
figure 5

The obtained PMFs increase with the decreasing twisted angles and the PMF areas are distributed near the AA regions with its maximum value occuring at the AA/AB transitions, which is highly consistent with our calculated results.

Full size image

Discussion

Collectively, the central findings achieved here may have important far-reaching implications. First of all, the existence of PMFs in TBG systems may finally offer testing grounds for observing quantum anomalous Hall effects as originally proposed26. Secondly, such PMFs are opposite for electrons at the K and K’ valleys, leading to clockwise or anticlockwise valley vortex currents around the AA regions, which can be exploited for establishing quantum valley Hall effects1. Thirdly, the present study has demonstrated that the atomic deformations play a critical role in determining the electronic properties of the low-angle TBG, which in turn should also influence the emergent physics of correlated electrons13,14,15,16,17,18,19. Furthermore, similar physical phenomena, including the appearance of PMFs, may also be present in other van der Waals systems (for example, transition metal dichalcogenides bilayers or heterostructures41,42,43), and their effects on the electronic properties remain to be fully explored.

Methods

Sample fabrication and characterization

The bilayer graphene samples were prepared by chemical vapor deposition on Cu–Ni alloy substrates using \({{\rm{CH}}}_{4}\) as the carbon feedstock and \({{\rm{H}}}_{2}\) as the reduction gas. Graphene growth was achieved at 1040 °C for 10–30 min after substrate was annealed for 30 min at 1000 °C under vacuum. The as-grown sample was then annealed in ultra-high vacuum at around 600 °C for several hours. The STM and STS measurements were conducted in a home-built STM with RHK R9 controller and RHK Pan-style scan head at 77 K with base pressure below \(1.0\, \times 1{0}^{-10}\ {\rm{Torr}}\). The scanning tips were obtained by electrochemical etching of tungsten wires. STS measurements were carried out with standard lock-in technique with bias modulation of 5 mV at 2023 Hz.

Tight-binding calculation

The full tight-binding model by considering \({p}_{z}\) orbitals of all carbon atoms in twisted bilayer graphene is used to compute the local density states and other electronic properties of all the presented calculations. The model is developed in20, and has been shown to reproduce accurately the differential conductance \({\mathrm{dI}}/{\mathrm{dV}}\) measured from low-angel twisted bilayer graphene. In this model, the hopping integral \({t}_{ij}\) between \({p}_{z}\) orbitals at site \(i\) and \(j\) is described by:

$${t}_{ij}={n}^{2}{V}_{pp\sigma }({r}_{ij})+(1-{n}^{2}){V}_{pp\pi }({r}_{ij}),$$

where \({r}_{ij}\) is the distance between site \(i\) and \(j\), \(n\) is the direction cosine along the direction perpendicular to the graphene layer \({{\boldsymbol{e}}}_{{\boldsymbol{z}}}\). The Slater and Koster parameters \({V}_{pp\sigma }\) and \({V}_{pp\pi }\) follow

$${V}_{pp\pi }({r}_{ij})= \, -{\gamma }_{0}{e}^{{q}_{\pi }(1-{r}_{ij}/d)}{F}_{c}({r}_{ij}),\\ {V}_{pp\sigma }({r}_{ij})= \, {\gamma }_{1}{e}^{{q}_{\sigma }(1-{r}_{ij}/h)}{F}_{c}({r}_{ij}),$$

where d = 1.42 Å is the nearest carbon-carbon distance, h = 3.349 Å is the interlayer distance, \({\gamma }_{0}\) and \({\gamma }_{1}\) were chosen to be 3.12 and 0.48 eV, respectively. \({q}_{\sigma }\) and \({q}_{\pi }\) satisfy \(\frac{{q}_{\pi }}{d}=\frac{{q}_{\sigma }}{h}=\) 2.218 Å−1, and \({F}_{c}\) is a smooth function \({F}_{c}(r)={(1+{e}^{(r-{r}_{c})/{l}_{c}})}^{-1}\), in which \({l}_{c}\) and \({r}_{c}\) were chosen to be 0.265 and 5.0 Å, respectively. The local density of states were calculated by using the recursion method in real space based on Lanczos algorithm34. To reach a precision in energy of 3.5 meV, we have computed 10,000 steps in the continued fraction built in a super-cell of 15 millions of orbitals with periodic boundary conditions.