Electrostatic superlattices on scaled graphene lattices

Electrostatic superlattices have been known to significantly modify the electronic structure of low-dimensional materials. Studies of graphene superlattices were triggered by the discovery of moiré patterns in van der Waals stacks of graphene and hexagonal boron nitride (hBN) layers a few years ago. Very recently, gate-controllable superlattices using spatially modulated gate oxides have been achieved, allowing for Dirac band structure engineering of graphene. Despite these rapid experimental progresses, technical advances in quantum transport simulations for large-scale graphene superlattices have been relatively limited. Here, we show that transport experiments for both graphene/hBN moiré superlattices and gate-controllable superlattices can be well reproduced by transport simulations based on a scalable tight-binding model. Our finding paves the way to tuning-parameter-free quantum transport simulations for graphene superlattices, providing reliable guides for understanding and predicting novel electric properties of complex graphene superlattice devices. https://doi.org/10.1038/s42005-020-0335-1 OPEN

H exagonal boron nitride (hBN), one of the most popular dielectric materials due to its atomic flatness 1 , has been playing a crucial role in studying two-dimensional materials using so-called van der Waals heterostructures 2 . The graphene/hBN moiré pattern 3 arising from the large-scale lattice interference led to the discovery of graphene superlattices. First experiments revealing new transport phenomena, such as the emergence of the Hofstadter butterfly, were reported in 2013 4-6 . In the following years, other exciting transport experiments have been reported [7][8][9][10][11][12][13] , as well as a dynamic band structure tuning 14,15 . More recently, another approach for inducing a superlattice potential in graphene has been demonstrated by using patterned dielectrics 16 , allowing for mini-band structure engineering. On the theory side, most works related to graphene superlattices focus either on calculations for the superlatticeinduced mini-band structures [17][18][19][20][21][22][23] , or on predicting transport properties by solving the Dirac equation with simplified superlattice model potential 24,25 . Quantum transport simulations considering realistic experimental conditions have been relatively rare in the literature 26,27 , not to mention a theory work combining transport simulations and mini-band structure calculations, together with transport experiments.
This article aims at providing a straightforward method to perform reliable quantum transport simulations for graphene superlattices. As shown in the following, our transport simulations based on the real-space Green's function method for twoterminal structures with the superlattice potential arising either from the graphene/hBN moiré pattern or from periodically modulated gating are consistent with transport experiments as well as mini-band structures based on the continuum model. Our method is applicable equally well to multi-terminal structures for simulating, for example, four-probe measurements using the Landauer-Büttiker approach 28 .

Results
Superlattice on a scaled graphene lattice. To perform quantum transport simulations for graphene working in real-space, the scalable tight-binding model 29 has been proved to be a very convenient numerical tool [30][31][32][33][34] : the physics of a real graphene system can be captured by a graphene lattice scaled by a factor of s such that the lattice spacing and nearest-neighbor hopping parameter are given by a = sa 0 and t = t 0 /s, respectively, where a 0 ≈ 0.142 nm and t 0 ≈ 3 eV are the tight-binding parameters for a genuine graphene lattice. The scaling is justified as long as a remains much shorter than all important physical length scales in the considered graphene system.
In dealing with graphene superlattices (such as the one sketched in Fig. 1a), the newly introduced physical length scale not previously considered 29 is the periodicity λ of the superlattice. The advantage of the scaling can be easily appreciated by comparing Fig. 1b, c: The former considers a genuine graphene lattice involving lots of carbon atoms, while the latter depicts a scaled graphene lattice (here s = 2 for illustrative purposes) involving a strongly reduced number of lattice sites that downscales with s 2 . As long as a ≪ λ is satisfied, a reasonably large area covering enough superlattice periods can be implemented in realspace quantum transport simulations (see Methods) to reveal transport properties arising from the superlattice effects.
In the following, we demonstrate the reliability of transport simulations for graphene superlattices based on scaled graphene lattices, considering both graphene/hBN moiré superlattice and gate-controlled superlattice.
Graphene/hBN moiré superlattice model potential. Formation of the moiré pattern due to the stacking of hBN and graphene lattices has been understood in one of the earliest experiments 3 . Following their model, the moiré pattern results in a triangular periodic scalar potential described by where V = 0.06 eV is the amplitude of the model potential and G 1 ðλ;θÞ is the reciprocal primitive vector of the moiré pattern corresponding to the primitive vector L 1 (λ, θ) = λ(cos θ, sin θ) in real space. The orientation angleθ and wavelengthλ are related to those in real-space throughθ ¼ θ þ π=2 andλ ¼ 4π= ffiffi ffi 3 p λ. The other two reciprocal vectors are given by G 2 ðλ;θÞ ¼ G 1 ðλ;θ þ π=3Þ and G 3 ðλ;θÞ ¼ G 1 ðλ;θ þ 2π=3Þ. Following Moon and Koshino 23 with the zigzag lattice direction arranged along the x axis, the moiré wavelength λ and the orientation angle θ of the pattern are given by and θ ¼ arctan respectively, where a g ¼ ffiffi ffi 3 p a 0 % 0:246 nm is the graphene lattice constant, ϵ = (a hBN − a g )/a g ≈ 1.81% is the lattice constant mismatch with a hBN ≈ 0.2504 nm the hBN lattice constant, and ϕ is the twist angle of the hBN lattice relative to the graphene lattice. An illustrative example with ϕ = 5°is sketched in Fig. 2a, where an overlay of U s (x, y) given by Eq. (1) is shown to match perfectly the lattice structure of the resulting graphene/hBN moiré pattern. For self-containment, λ and θ as functions of the twist angle ϕ are plotted in Fig. 2b, c, respectively, where the hollow squares mark Transport experiment. To confirm the validity of our simulation scheme illustrated above, the most convincing way is to fabricate a real device, perform a transport measurement, and finally simulate the experiment. Our single-gated two-terminal device was fabricated with a graphene/hBN stack (see the inset of Fig. 2d and Methods) and measured at low temperature (≈4.1 K), using standard low-frequency (≈13 Hz) lock-in technique. Figure 2d shows the two-terminal differential conductance of our sample as a function of the back gate voltage V bg . Conductance dips at around −60 and +58 V are basic characteristics of the moiré superlattice potential [4][5][6] , showing that our graphene and hBN lattices are nearly aligned. By analyzing from the Brown-Zak oscillation 12,13 , the moiré wavelength λ ≈ 10.9 nm of our graphene/hBN sample was deduced. Note that although our graphene sample was encapsulated by two hBN layers, the single dip structure (in each of the electron and hole branches) of our conductance measurement (Fig. 2d) suggests that only one hBN layer exhibits a measurable moiré superlattice effect arising from a small twist angle, contrary to a recent report on doubly aligned hBN/graphene/hBN heterostructures 35 . In the following simulations, therefore, we model our experiment with graphene/hBN, instead of hBN/graphene/hBN. Transport simulation. To simulate such a conductance measurement, we have calculated the transmission T(E) as a function of Fermi energy E at zero temperature, and hence the conductance G(E) = (2e 2 /h)T(E), based on an s = 4 tight-binding model Hamiltonian, for a two-terminal device similar to considers leads with on-site and Fermi energies identical to those at the lattice sites attached to the leads; see Fig. 1a. In this way, the Fermi energy in the leads "floats" with that in the scattering region such that the interface between the lead and the scattering region is as transparent as possible. Compared to Fig. 2d, the electron-hole asymmetry is less pronounced due to the simple model of Eq. (1) from Yankowitz et al. 3 . Although accounting for an electron doping from the metal contacts simply by fixing the Fermi energy in the leads with a positive value can make the conductance curve (black dashed curve in Fig. 2e with Fermi energy 0.32 eV in the leads) even more similar to the experiment, the nature of the electron-hole asymmetry observed in Fig. 2d comes from more subtle interactions between graphene and hBN lattices which require more advanced model Hamiltonians [21][22][23] beyond the simple model of Eq. (1). We will not further address the electron-hole asymmetry but continue our discussions with calculations based on leads with "floating" Fermi energies.
Density of states and mini-band structures. Without transforming to the gate voltage axis, the original conductance data of Fig. 2e as a function of energy is reported in Fig. 2f with a wider energy range up to ±0.4 eV. Compared to the density of states (Fig. 2g) and the band structure ( Fig. 2h) which are calculated based on the same moiré superlattice model potential but within the continuum model (see Methods), consistent features in the energy spectrum can be seen. In view of Fig. 2d-h, our calculations significantly capture basic properties of the graphene/hBN moiré superlattice, at least at zero magnetic field.
Dependence on moiré orientation angle. In all the transport simulations for the moiré superlattice shown in the main text of this article, we have considered realistic moiré orientation angle θ prescribed by Eq. (2b), relative to the zigzag direction. We have further tested that the simulation results are insensitive to θ. See Supplementary Fig. 2 and Supplementary Note 3. Magnetotransport experiment and simulation. We continue our comparison of the experimentally measured and theoretically calculated conductance G(V bg ) with finite magnetic field B perpendicular to the graphene plane, which can be modeled by associating the Peierls phase 28,36   This suggests that the neglected higher-order terms of a more complete model Hamiltonian 21-23 become important when the magnetic field is strong. Interestingly, we note that in the theory map of Fig. 3b, some unusual plateaus in the energy range around the electron-branch secondary Dirac point can be observed. We magnify this region in Fig. 3e with a horizontal line cut shown in Fig. 3f, where the gate voltage range showing quantized conductance plateaus is highlighted by a yellow background.
This V bg range, transformed back to the energy through Eq. (S1) of Supplementary Note 2, corresponds to an energy window where part of the electron branch of the secondary Dirac cones at M points of the superlattice mini-Brillouin zone are completely isolated. The respective energy window is highlighted also by yellow in Fig. 2h. Since there are effectively three such Dirac cones (six cones on six M points within each mini-Brillouin zone but each cone shared by two neighboring mini-Brillouin zones), the degeneracy factor is expected to be 3 × 2 × 2 = 12 with ×2 accounting for spin and another ×2 for valley. Indeed, in the quantum Hall regime, the calculated conductance is quantized to 6, 18, 30, 42 e 2 /h as shown in Fig. 3f. Outside this energy (and hence back gate voltage) range, the higher-order Dirac cones are always mixed with background bands, so that no quantized conductance is observed. However, such a special energy window leading to the 12-fold-degeneracy of the Landau levels at the secondary DP is never observed in transport experiments with graphene/hBN moiré superlattices 4,5,10-12 , including ours shown in Fig. 3a, indicating once again that the simplified model 3 of Eq. (1) containing only the electrostatic scalar potential term is not sufficient to capture transport properties of graphene/hBN moiré superlattices at higher magnetic fields.
At this stage, we comment that the periodic scalar potential enters the tight-binding model Hamiltonian through the on-site energy term (see Methods), so that it is readily compatible with the scaling method 29 . Including higher-order terms 21-23 of the graphene/hBN moiré superlattice should be possible but is beyond the scope of the present study. As we will see below, when the graphene superlattice potential arises solely from the electrostatic gating, our method becomes even more precise because in such systems the scalar potential is the only term comprising the superlattice.
Gate-controllable superlattices. To observe any superlattice effects in graphene, the mean free path must exceed enough periods of the superlattice potential. This means that either the sample quality must be extraordinary, or the superlattice periodicity must be short enough. When the periodicity is too short, however, the resulting extra Dirac cones appear at too high energy, exceeding the experimentally reachable range. This is why the discovery of the graphene/hBN moiré pattern 3 led to studies on the graphene superlattice of its first kind-the periodicity corresponding to small twist angles turns out to be naturally in a suitable range for experiments; see Fig. 2b. A more flexible approach to design artificial graphene superlattice structures for band structure engineering was pursued with the realization of electrostatic gating schemes 16 . To create an externally controllable periodic potential, the most intuitive way is to pattern an array of periodic fine metal gates on top of the graphene sample 37,38 . However, due to technical difficulties such as instabilities of nanometer-scale local gates, the low adhesion between metal gates and the inert hBN, etc., such superlattice graphene devices often suffer the problem of very low sample yield 39 . The basic idea of the new technical breakthrough is to keep the hBN/graphene/hBN sandwich intact, while periodically modulating the gate capacitance. This can be achieved either by using few-layer graphene as a local gate which is subsequently etched with a periodic pattern 39,40 , or by etching the dielectric layer with a periodic pattern using a standard uniform back gate underneath the modulated substrate 16 . In the following, we revisit the experiment of the latter by Forsythe et al. 16 and reproduce the observed transport features by our transport simulations.
Electrostatic simulation. Following the geometry of the device subject to a square superlattice potential with periodicity λ = 35 nm presented in the work by Forsythe et al. 16 , we have performed our own electrostatic simulation to obtain the back gate capacitance showing periodic spatial modulation. We consider an hBN/  16 . The bottom gate capacitance therefore shows a spatial modulation with a square lattice symmetry, as shown in the lower left inset of Fig. 4a.
With the electrostatically simulated position-dependent back gate capacitance per unit area C bg (x, y), contributing carrier density n bg (x, y) = [C bg (x, y)/e]V bg , together with the uniform n tg , the resulting superlattice potential is given by U s ðx; yÞ ¼ Àsgn½nðx; yÞ hv F ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi π nðx; yÞ j j p with n = n bg + n tg , in order to set the global transport Fermi level at zero 41 . Slightly different from the case of the graphene/hBN moiré superlattice where the model potential U s (x, y) given by Eq. (1) is independent of the gating, we consider U(x, y) = U s (x, y) for the on-site energy term (4), and implement it in the tight-binding Hamiltonian (3) with s = 6 (such that a = sa 0 ≈ 0.85 nm ≪ λ = 35 nm) to perform quantum transport simulations over a two-terminal structure with L = 420 nm and W = 385 nm; see the upper right inset of Fig. 4a.
Transport simulation. To compare with the resistance measurements reported by Forsythe et al. 16 , we plot the inverse transmission 1/T as a function of n tg and V bg in the main panel of Fig. 4a, where most areas show high transmission (white regions correspond to low 1/T). Along the diagonal dark thick line showing high 1/T values due to the main Dirac point, multiple satellite peaks can be seen when increasing |V bg | and hence the magnitude of the square superlattice potential, signifying the emerging multiple extra Dirac points due to the gate-controlled square superlattice potential. Exemplary line cuts are plotted in Fig. 4b-d to show clearly the single-and multiple-peak structures, in excellent agreement with the experiment (see Figure 2b of Forsythe et al. 16 ). Note that the peak positions of the secondary Dirac points at the electron and hole side seen in Fig. 4b, d are not symmetric due to the asymmetric superlattice potential shape obtained from the electrostatic simulation. When implementing a model periodic function symmetric in its potential profile, the density spacing from the main Dirac point to the secondary Dirac points at electron and hole sides becomes identical (not shown).
Gate-dependent mini-band structures. We have also checked the consistency between the calculated mini-band structures and the simulated inverse transmission. Overall, we obtain band structures similar to those reported in the work by Forsythe et al. 16 , but since each (n tg , V bg ) point corresponds to a different U s (x, y) profile and hence a different mini-band structure, an overview consistency-check like in Fig. 2f-h is technically not possible. Instead, the consistency can be checked by comparing the 1/T peaks and their corresponding mini-band structure around E = 0. We have chosen three particular (n tg , V bg ) configurations corresponding to the three black arrows in Fig. 4b marking three of the 1/T peaks, at which the E = 0 Fermi level is expected to hit either the main or the extra Dirac points.
These mini-band structures, along with the actual U s (x, y) profiles implemented individually in the continuum model (see Methods) are shown in Fig. 4e-g. Going from low to high n tg (left, middle, and right black arrow in Fig. 4b), the highest filled energy rises relative to the main Dirac point, corresponding to the sinking of the whole band structure due to our choice of fixing the Fermi level at E = 0 (Fig. 4e-g). As expected, the highest peak in Fig. 4 Gate-controlled square superlattice on graphene. a Inverse transmission 1/T as a function of top-gate-contributed carrier density n tg and back gate voltage V bg for the simulated virtual device similar to one of those reported by Forsythe et al. 16 . Upper inset: Device geometry with the modulated back gate capacitance profile C bg showing a spatial modulation of periodicity λ = 35 nm. The boxed region is magnified in the lower inset with the color bar calibrating the value of C bg in units of 10 10 cm −2 V −1 . The three horizontal color lines marked on the main panel of (a) correspond to the line cuts shown in (b-d). e-g Mini-band structures based on the continuum model, showing six bands closest to the main Dirac point, labeled from high to low energy by C 3 (yellow), C 2 (blue), C 1 (red), V 1 (green), V 2 (orange), and V 3 (purple), implementing U s (x, y) displayed over a range of 2λ × 2λ above each main panel. Panels (e-g) correspond to (n tg , V bg ) configurations labeled by the left, middle, and right black arrows marked in (b), respectively. The black dashed boxes mark the E = 0 Fermi level hosting low-bias transport, where the intersected sub-band Fermi contours are shown below the corresponding mini-band structure: V 1 (green) and V 2 (orange) in (e), C 1 (red) in (f), and C 1 (red) and C 2 (blue) in (g). Symmetry points on the square mini-Brillouin zone (Γ for the central point, M for the corner points, and X for the midpoints on the edges) are labeled in a way that avoids disturbing the Fermi contours shown.  Fig. 4e, g, the two satellite 1/T peaks seen in Fig. 4b are mainly contributed by the secondary Dirac points at X, labeling the midpoints on the edges of the square mini-Brillouin zone. Note that the mini-band structures shown in Fig. 4e-g, though corresponding to an increasing uniform n tg , do not exhibit simply an energy shift without changing the band shape. Compare, for example, the shapes of the lowest shown subbands (purple).

Discussion
We have shown that quantum transport simulations based on the scalable tight-binding model 29 correctly capture transport properties of electrostatic graphene superlattices. In the case of graphene/hBN moiré superlattice, the consistency of our simulation and experiment at zero and low magnetic field is rather satisfactory but breaks down at too strong magnetic field due to the adopted simple moiré model potential 3 that neglects higher order terms. In the other case of gated superlattices, without such higher order terms the simulations are expected to be precise for all magnetic field range. Compared to the recent transport experiment on a gate-controlled square superlattice device reported by Forsythe et al. 16 , our simulations show an excellent agreement in revealing the emergence of multiple extra Dirac cones at zero magnetic field. Transport simulations at finite magnetic field for the gated superlattices are expected to reveal also consistent behaviors compared to the experiment, but are left as future work.
Our finding significantly lowers the computation burdens and hence paves the way to tuning-parameter-free quantum transport simulations for graphene superlattices, providing reliable guides for understanding and predicting novel electric properties of complex graphene superlattice devices. We note some recent studies working on developing numerical techniques that allow large-scale efficient transport simulations [42][43][44] , but properly scaling the graphene lattice seems to be of least technical complexity and is readily applicable to anyone who is familiar with quantum transport using, for example, real-space Green's function methods 28 or the open-source python package KWANT 45 .

Methods
Real-space Green's function method. The model Hamiltonian including the superlattice potential U s (x, y) using the scaled graphene lattice can be written as where the operator c j (c y j ) annihilates (creates) an electron at site r j = (x j , y j ). The first term in Eq. (3) represents the clean part of the Hamiltonian which contains nearest neighbor hoppings summing over site indices i and j with 〈i, j〉 standing for |r i − r j | = a, and the second term is the on-site energy containing the superlattice potential U s (x, y) smoothed by a model function F s (x, where ' s is a smoothing parameter typically taken as ' s ¼ λ=4. The purpose of smearing off the superlattice potential function U s to zero at a distance d (typically taken as λ) away from the edges and the leads (see Fig. 1a and its inset) is to avoid any spurious effects due to the combination of the superlattice potential and the physical edges of the graphene lattice, as well as to avoid oversized unit cells for the lead self-energies. Any contributions to the on-site energy term other than the superlattice potential are collected in the U 0 term in Eq. (4). With the model Hamiltonian Eq. (3) constructed, together with self-energies Σ 1 and Σ 2 describing the attached two leads (following, for example, Wimmer 46 ), the retarded Green's function at energy E is given by leading to the transmission function where Γ j ¼ iðΣ j À Σ y j Þ with j = 1, 2 is the broadening function. In the lowtemperature low-bias limit, the conductance across the modeled scattering region is given by the Landauer formula G = (2e 2 /h)T, where the factor of 2 accounts for the spin degeneracy. For a pedagogical introduction to the above outlined real-space Green's function, see, for example, Datta 28 . Note that in most simulations, the full matrix of Eq. (5) is not needed, suggesting that a partial inversion should be implemented in the numerics to avoid wasting computer memories and CPU time. On the other hand, the matrix version of the Fisher-Lee relation (6) can be implemented as the way it reads.
Device fabrication. Our two-terminal device contains a hBN/graphene/hBN stack on a Si/SiO 2 substrate, where the crystallographic axis of the graphene flake is aligned with respect to one of the hBN flakes. Electric contact to the graphene is made from the edge of the mesa 47 with self-aligned Ti/Al electrodes. We use the Si wafer as an overall back gate with a two-layer dielectric consisting of SiO 2 with thickness d SiO 2 = 300 nm and the bottom hBN flake with thickness d hBN = 20 nm. A typical exemplary junction similar to the measured device is shown by the atomic-force microscope (AFM) image in the inset of Fig. 2d and marked by the white dashed box.
Continuum model for mini-bands and density of states. To calculate the miniband structure of graphene in the presence of a superlattice potential U s (r), we consider an infinitely large two-dimensional pristine graphene described by H 0 in k-space. Following Park et al. 17 , we start with the continuum model Hamiltonian near the K valley: where the first term is H 0 and superlattice potential in the second term is treated as a perturbation. In Eq. (7), the product of the reduced Planck constant ħ and Fermi velocity v F is related to the tight-binding parameters through ħv F = (3/2)ta, and the two-dimensional wave vector (k x , k y ) = k is small relative to the K point. Using the eigenstates of H 0 ðkÞ as a new basis, we solve the eigenvalue problem of Eq. (7) to obtain a set of linear equations: ½E À ε s ðkÞc s ðkÞ ¼ X s0;G 1 þ ss0e Àiθ k;kÀG 2 UðGÞc s0 ðk À GÞ; where E is the energy eigenvalue of the graphene superlattice, ε s (k) = sħv F k is the eigenenergy of H 0 ðkÞ associated with the s branch (s = 1 for electron above the Dirac point and s = −1 for hole below the Dirac point), UðGÞ is the Fourier component of U s (r) with the reciprocal lattice vector G = m 1 G 1 + m 2 G 2 of the superlattice potential, θ k,k−G is the angle from k − G to k, and c s (k) are the expansion coefficients of the pristine graphene eigenstates. The infinite-dimensional matrix spanned by the states with wave vectors P m 1 ;m 2 k þ m 1 G 1 þ m 2 G 2 in Eq. (8) allows for solving for E and hence calculating the band structure. Since we focus on the low-energy region, a matrix involving states with |m 1 | ≤ 3 and |m 2 | ≤ 3 is found to be sufficient to attain the convergence of the band structure. With the eigenenergies obtained, the density of states D as a function of energy can be calculated by where the integration is taken over the first Brillouin zone.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.