Electrostatic Superlattices on Scaled Graphene Lattices

TL;DR

This paper develops a scalable tight-binding model to simulate quantum transport in graphene superlattices, successfully matching experimental results for gate-controlled superlattices and exploring Dirac cone emergence.

Contribution

The study introduces a large-scale tight-binding approach for graphene superlattices, accurately modeling transport phenomena and revealing new Dirac cones at strong superlattice modulations.

Findings

Good agreement with experiments at low magnetic fields for moiré superlattices

Identification of multiple extra Dirac cones at high superlattice modulation

Simulation accuracy for gate-controlled superlattices without higher-order terms

Abstract

A scalable tight-binding model is applied for large-scale quantum transport calculations in clean graphene subject to electrostatic superlattice potentials, including two types of graphene superlattices: moir\'e patterns due to the stacking of graphene and hexagonal boron nitride (hBN) lattices, and gate-controllable superlattices using a spatially modulated gate capacitance. In the case of graphene/hBN moir\'e superlattices, consistency between our transport simulation and experiment is satisfactory at zero and low magnetic field, but breaks down at high magnetic field due to the adopted simple model Hamiltonian that does not comprise higher-order terms of effective vector potential and Dirac mass terms. In the case of gate-controllable superlattices, no higher-order terms are involved, and the simulations are expected to be numerically exact. Revisiting a recent experiment on graphene subject to a gated square superlattice with periodicity of 35 nm, our simulations show excellent agreement, revealing the emergence of multiple extra Dirac cones at stronger superlattice modulation.

Equations

$$H=-t\sum _{⟨i,j⟩}{c}_i^{†}{c}_j+\sum _{j}U({\mathbf{r}}_j)c_j^{†}{c}_j,$$

$$U(x,y)=U_s(x,y)F_s(x,y)+U_0(x,y),$$

$$G_r(E)=\frac{\text{openone}}{E-[H+{\Sigma }_1(E)+{\Sigma }_2(E)]},$$

$$T(E)=\mathrm{Tr}[{\Gamma }_1(E)G_r(E){\Gamma }_2(E)G_r^{†}(E)],$$

$$\mathcal{H}(\mathbf{k})=\hslash v_F\left(\begin{array}{cc}0& -i{k}_x-k_y\\ i{k}_x-k_y& 0\end{array}\right)+U_s(\mathbf{r})\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),$$

$$[E-{\epsilon }_s(\mathbf{k})]c_s(\mathbf{k})=\sum _{s^{\prime },\mathbf{G}}\frac{1+s{s}^{\prime }{e}^{-i{\theta }_{\mathbf{k},\mathbf{k}-\mathbf{G}}}}{2}\mathcal{U}(\mathbf{G})c_{s^{\prime }}(\mathbf{k}-\mathbf{G}),$$

$$D(\epsilon )=\frac{2}{(2\pi {)}^2}{\int }_{1_{\mathrm{st}}\mathrm{BZ}}\delta (E-\epsilon )d{k}_{x}d{k}_y,$$

$$U_s(\mathbf{r})=V\sum _{j=1,2,3}\cos \left({\mathbf{G}}_j\cdot \mathbf{r}\right),$$

$$\begin{array}{rl}\lambda & =\frac{1+ϵ}{\sqrt{{ϵ}^2+2(1+ϵ)(1-\cos \varphi )}}{a}_G\\ \theta & =\arctan \frac{-\sin \varphi }{1+ϵ-\cos \varphi }\end{array},$$

$$V_{\mathrm{bg}}=\frac{e}{\pi C}{\left(\frac{E}{\hslash v_F}\right)}^2\mathrm{sgn}(E),$$

Full text

I Introduction

Graphene, a single layer of carbon atoms arranged in a honeycomb lattice, was first successfully isolated from a single crystal in 2004 Novoselov et al. (2004); Geim and Novoselov (2007), which subsequently triggered further investigations on the intriguing properties of relativistic Dirac fermions Castro Neto et al. (2009). However, to further uncover more novel electronic properties of the first truly two-dimensional material, the limited mobility of graphene on standard SiO2 substrates turned out to be the main factor restricting mean-free path and phase-coherence length. The discovery of hexagonal boron-nitride (hBN) as an ideal atomically flat substrate for graphene Dean et al. (2010) boosted the development of high-quality graphene devices. Fabrication of devices that involves the encapsulation of graphene between two thin hBN multilayers has become a standard protocol since then Yankowitz et al. (2019). At the same time, the combination of these two different 2D materials in a so-called van der Waals heterostructure Geim and Grigorieva (2013) led to the subsequent discovery of the graphene/hBN moiré pattern Yankowitz et al. (2012) arising from the large-scale lattice interference due to the slight lattice constant mismatch.

At small twist angles, the resulting moiré pattern provides a natural source of superlattice potential on graphene with periodicity in the order of $\sim 10\mathord{\thinspace\rm nm}$, leading to the formation of new superlattice minibands in the electronic band structure of graphene at energies reachable by standard electrostatic gating. First experiments revealing new transport phenomena (such as the emergence of the Hofstadter butterfly) were reported in 2013 Ponomarenko et al. (2013); Dean et al. (2013); Hunt et al. (2013). In the following years, other exciting transport experiments have been reported Gorbachev et al. (2014); Wang et al. (2015a, b); Lee et al. (2016); Handschin et al. (2017); Krishna Kumar et al. (2017); Chen et al. (2017), as well as a dynamic band structure tuning Yankowitz et al. (2018); Ribeiro-Palau et al. (2018). More recently, another approach for inducing a superlattice potential in graphene has been demonstrated by using patterned dielectrics Forsythe et al. (2018). Such a gate-tunable superlattice structure allows for the design of arbitrary superlattice geometries with defined periodicity but suffers from technical restrictions.

On the theory side, most works related to graphene superlattices focus either on calculations for the superlattice-induced mini-band structures Park et al. (2008); Brey and Fertig (2009); Wu et al. (2012); Ortix et al. (2012); Kindermann et al. (2012); Wallbank et al. (2013); Moon and Koshino (2014), or on predicting transport properties by solving the Dirac equation with oversimplified superlattice model potential Barbier et al. (2010); Burset et al. (2011). On the other hand, quantum transport simulations considering realistic experimental conditions have been relatively rare in the literature Diez et al. (2014); Hu et al. (2018), not to mention a theory work combining quantum transport simulations and mini-band structure calculations, together with transport experiments. This work aims at providing a straightforward method to perform reliable quantum transport simulations, covering both graphene/hBN moiré superlattices and gate-induced superlattices. As shown in the following, our transport simulations based on the real-space Green’s function method for two-terminal structures as sketched in Figure 1(a) 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 Datta (1995).

This paper is organized as follows. In section II, we briefly introduce the theoretical methods used in section III, including graphene/hBN moiré superlattice (subsection III.1) and gate-induced superlattice (subsection III.2). A summary of this work is provided in section IV.

II Theoretical Approaches

II.1 Real-space tight-binding model for quantum transport

To perform quantum transport simulations for graphene working in real-space, the scalable tight-binding model Liu et al. (2015) has been proved to be a very convenient numerical trick (see, for example, Rickhaus et al. (2015); Mreńca-Kolasińska and Szafran (2016); Kolasiński et al. (2017); Petrović et al. (2017); Ma et al. (2018)): the physics of a graphene system can be captured by a graphene lattice scaled by a factor of $s_{f}$ such that the lattice spacing is given by $a=s_{f}a_{0}$ with $a_{0}\approx 0.142\mathord{\thinspace\rm nm}$ the carbon-carbon distance and the nearest neighbor hopping is $t=t_{0}/s_{f}$ with $t_{0}\approx 3\mathord{\thinspace\rm eV}$ the hopping parameter for a genuine graphene lattice, as long as the scaled lattice spacing $a$ remains much shorter than all important physical length scales in the graphene system of interest.

In dealing with graphene superlattices, the newly introduced physical length scale not mentioned in Ref. Liu et al., 2015 is the periodicity $\lambda$ of the superlattice. The advantage of the scaling can be easily appreciated by comparing Figure 1(b) and Figure 1(c): The former considers a genuine graphene lattice involving lots of carbon atoms, while the latter considers a scaled graphene lattice (here $s_{f}=2$ for illustrative purposes) involving a much reduced number of lattice sites. As long as $a\ll\lambda$ is satisfied, a reasonably large area covering enough superlattice periods can be implemented in real-space quantum transport simulations to reveal transport properties arising from the superlattice effects.

The model Hamiltonian including the superlattice potential $U_{s}(x,y)$ using the scaled graphene lattice can be written as

[TABLE]

where the operator $c_{j}$ ($c_{j}^{\dagger}$) annihilates (creates) an electron at site $\mathbf{r}_{j}=(x_{j},y_{j})$. The first term in Eq. (1) represents the clean part of the Hamiltonian which contains nearest neighbor hoppings summing over site indices $i$ and $j$ with $\langle i,j\rangle$ standing for $|\mathbf{r}_{i}-\mathbf{r}_{j}|=a$, and the second term is the on-site energy

[TABLE]

containing the superlattice potential $U_{s}(x,y)$ smoothed by a model function $F_{s}(x,y)=f_{s}(x,-L/2+d)f_{s}(-x,L/2-d)f_{s}(y,-W/2+d)f_{s}(-y,W/2-d)$ with $f_{s}(z,z_{0})=\{1+\tanh[(z-z_{0})/\ell_{s}]\}/2$, where $\ell_{s}$ is a smoothing parameter typically taken as $\ell_{s}=\lambda/4$. The purpose of smearing off the superlattice potential function $U_{s}$ to zero at a distance $d$ (typically taken as $\lambda$) away from the edges and the leads [see Figure 1(a) 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. (2).

With the model Hamiltonian Eq. (1) constructed, together with self-energies $\Sigma_{1}$ and $\Sigma_{2}$ describing the attached two leads (following, for example, Ref. Wimmer, 2008), the retarded Green’s function at energy $E$ is given by

[TABLE]

leading to the transmission function

[TABLE]

where $\Gamma_{j}=i(\Sigma_{j}-\Sigma_{j}^{\dagger})$ with $j=1,2$ is the broadening function. In the low-temperature 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, Ref. Datta, 1995. Note that in most simulations, the full matrix of Eq. (3) 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 (4) can be implemented as the way it reads.

II.2 Continuum model for mini-bands and density of states

To calculate the mini-band structure of graphene in the presence of a superlattice potential $U_{s}(\mathbf{r})$, we consider an infinitely large two-dimensional pristine graphene described by $\mathcal{H}_{0}$ in $k$-space. Following Ref. Park et al., 2008, we start with the continuum model Hamiltonian near the $K$ valley:

[TABLE]

where the first term is $\mathcal{H}_{0}$ and superlattice potential in the second term is treated as a perturbation. In Eq. (5), the product of the reduced Planck constant $\hbar$ and Fermi velocity $v_{F}$ is related to the tight-binding parameters through $\hbar v_{F}=(3/2)ta$, and the two-dimensional wave vector $(k_{x},k_{y})=\mathbf{k}$ is small relative to the $K$ point. Using the eigenstates of $H_{0}(\mathbf{k})$ as a new basis, we solve the eigenvalue problem of Eq. (5) to obtain a set of linear equations:

[TABLE]

where $E$ is the energy eigenvalue of the graphene superlattice, $\varepsilon_{s}(\mathbf{k})=s\hbar v_{F}k$ is the eigenenergy of $\mathcal{H}_{0}(\mathbf{k})$ associated with the $s$ branch ($s=1$ for electron above the Dirac point and $s=-1$ for hole below the Dirac point), $\mathcal{U}_{\mathbf{G}}$ is the Fourier component of $U_{s}(\mathbf{r})$ with the reciprocal lattice vector $\mathbf{G}=m_{1}\mathbf{G}_{1}+m_{2}\mathbf{G}_{2}$ of the superlattice potential, $\theta_{\mathbf{k},\mathbf{k}-\mathbf{G}}$ is the angle from $\mathbf{k}-\mathbf{G}$ to $\mathbf{k}$, and $c_{s}(\mathbf{k})$ are the expansion coefficients of the pristine graphene eigenstates.

The infinite-dimensional matrix spanned by the states with wave vectors $\sum_{m_{1},m_{2}}\mathbf{k}+m_{1}\mathbf{G}_{1}+m_{2}\mathbf{G}_{2}$ in Eq. (6) 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}|\leq 3$ and $|m_{2}|\leq 3$ is found to be sufficient to attain the convergence of the band structure.

The density of states $D$ as a function of energy can be calculated by

[TABLE]

where the integration is taken over the first Brillouin zone. Since $D$ is proportional to the number of energy eigenstates, it can be used to compare with the transport calculations.

III Electrostatic superlattices in graphene

III.1 Graphene/hBN moiré superlattices

Formation of the moiré pattern due to the stacking of hBN and graphene lattices has been understood in one of the earliest experiments Yankowitz et al. (2012). Following their model, the moiré pattern results in a triangular periodic scalar potential described by

[TABLE]

where $V=0.06\mathord{\thinspace\rm eV}$ is the amplitude of the model potential and $\mathbf{G}_{1}(\tilde{\lambda},\tilde{\theta})$ is the reciprocal primitive vector of the moiré pattern corresponding to the primitive vector $\mathbf{L}_{1}(\lambda,\theta)=\lambda(\cos\theta,\sin\theta)$ in real space. The orientation angle $\tilde{\theta}$ and wavelength $\tilde{\lambda}$ are related with those in real-space through $\tilde{\theta}=\theta+\pi/2$ and $\tilde{\lambda}=4\pi/\sqrt{3}\lambda$. The other two reciprocal vectors are given by $\mathbf{G}_{2}(\tilde{\lambda},\tilde{\theta})=\mathbf{G}_{1}(\tilde{\lambda},\tilde{\theta}+\pi/3)$ and $\mathbf{G}_{3}(\tilde{\lambda},\tilde{\theta})=\mathbf{G}_{1}(\tilde{\lambda},\tilde{\theta}+2\pi/3)$. Following Moon and Koshino (2014) with the zigzag lattice direction arranged along the $x$ axis, the moiré wavelength $\lambda$ and the orientation angle $\theta$ of the pattern are given by

[TABLE]

where $a_{G}=\sqrt{3}a_{0}\approx 0.246\mathord{\thinspace\rm nm}$ is the graphene lattice constant, $\epsilon=(a_{\mathrm{hBN}}-a_{G})/a_{G}\approx 1.81\%$ is the lattice constant mismatch with $a_{\mathrm{hBN}}\approx 0.2504\mathord{\thinspace\rm nm}$ the hBN lattice constant, and $\phi$ is the twist angle of the hBN lattice relative to the graphene lattice. An illustrative example with $\phi=5^{\circ}$ is sketched in Fig. 2(a), where an overlay of $U_{s}(x,y)$ given by Eq. (8) is shown to match perfectly the lattice structure of the resulting graphene/hBN moiré pattern. For completeness, $\lambda$ and $\theta$ as functions of the twist angle $\phi$ are plotted in Figure 2(b) and (c), respectively, where the hollow squares mark the $\phi=5^{\circ}$ example of Figure 2(a) and the hollow circles mark the $\phi=0.9^{\circ}$ case corresponding to our transport experiments and simulations to be elaborated below.

According to Ref. Park et al., 2008, the expected secondary Dirac points in the presence of a triangular superlattice potential occur at the $M$ points of the hexagonal mini-Brillouin zone in $k$ space, with a distance $\tilde{\lambda}/2$ to the $\Gamma$ point (where the normal Dirac point resides). Taking $k=\sqrt{\pi|n|}$ as an estimate for the corresponding carrier density $n$, this suggests that the example of Fig. 1(a) with $\phi=5^{\circ}$ requires a density above $5\times 10^{13}\mathord{\thinspace\rm cm^{-2}}$, which is beyond a reasonable density range from typical electrostatic gating. Indeed, observable graphene/hBN superlattice effects are typically found only in nearly aligned graphene/hBN stacks with a very small twist angle. On the other hand, the $M$ points for the case of $\phi=0.9^{\circ}$ are expected at a density no more than $4\times 10^{12}\mathord{\thinspace\rm cm^{-2}}$, lying in the typical density range using standard electrostatic gating.

III.1.1 Transport at zero magnetic field

To test the validity of the quantum transport simulation illustrated in subsection II.1 using the above moiré superlattice model potential Eq. (8), we compare our simulations with the experimental results obtained from a two-terminal device based on a hBN/graphene/hBN stack on a Si/SiO2 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 Wang et al. (2013) with self-aligned Ti/Al electrodes. We use the Si wafer as an overall back gate with a two-layer dielectric consisting of SiO2 with thickness $d_{\mathrm{SiO_{2}}}=300\mathord{\thinspace\rm nm}$ and the bottom hBN flake with thickness $d_{\mathrm{hBN}}=20\mathord{\thinspace\rm nm}$. A typical exemplary junction similar to the measured device is shown by the atomic-force mircroscope (AFM) image in the inset of Fig. 2(d) and marked by the white dashed box. Figure 2(d) shows the two-terminal differential conductance of our sample as a function of the back gate voltage $V_{\mathrm{bg}}$, measured at low temperature ($\approx 4.1\mathord{\thinspace\rm K}$), using standard low-frequency ($\approx 13\mathord{\thinspace\rm Hz}$) lock-in technique.

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 a $s_{f}=4$ tight-binding model Hamiltonian, for a two-terminal device similar to Figure 1(a) with $L=W=500\mathord{\thinspace\rm nm}$, implementing the moiré model potential (8) with a twist angle $\phi$. To compare with the experiment on the same voltage axis, we adopt the parallel-plate capacitor formula for the carrier density, $n=(C/e)V_{\mathrm{bg}}$, and relate $n$ with the Fermi energy through $E=\mathop{\rm sgn}\nolimits(n)\hbar v_{F}\sqrt{\pi|n|}$. These give us

[TABLE]

where $C$ is the back gate capacitance per unit area with $C/e=(\epsilon_{0}/e)(d_{\mathrm{SiO_{2}}}/3.9+d_{\mathrm{hBN}}/4.2)^{-1}\approx 6.77\times 10^{10}\mathord{\thinspace\rm cm^{-2}V^{-1}}$. Using this capacitance value, the exhibiting conductance dips at $V_{\mathrm{bg}}=+54.5\mathord{\thinspace\rm V}$ and $-59.9\mathord{\thinspace\rm V}$ observed in Figure 2(d) correspond to densities about $+3.7\times 10^{12}\mathord{\thinspace\rm cm^{-2}}$ and $-4.0\times 10^{12}\mathord{\thinspace\rm cm^{-2}}$, respectively, so that the twist angle is estimated to lie in the range of $0.85^{\circ}<\phi<0.95^{\circ}$. Indeed, when choosing $\phi=0.9^{\circ}$ for the moiré model potential Eq. (8), the simulated conductance $G(V_{\mathrm{bg}})$ transformed from $G(E)$ and reported in Figure 2(e) is found to show excellent agreement with the experiment Figure 2(d) in the positions of the conductance dips.

Note that the red curve reported in Figure 2(e) considers leads with on-site and Fermi energies identical to those at the attaching lattice sites; see Figure 1(a). Compared to Figure 2(d), the electron-hole asymmetry is less pronounced due to the simple model of Eq. (8) from Ref. Yankowitz et al., 2012. Although accounting for an electron doping from the metal leads simply by fixing the Fermi energy in the leads with a positive value can make the conductance curve [black dashed curve in Figure 2(e) with Fermi energy $0.32\mathord{\thinspace\rm eV}$ in the leads] even more similar to the experiment, the nature of the electron-hole asymmetry observed in Figure 2(d) comes from higher-order terms such as the effective vector potential and Dirac mass terms Kindermann et al. (2012); Wallbank et al. (2013); Moon and Koshino (2014) needed for a model Hamiltonian that can better describe the graphene/hBN moiré superlattice. We will continue our discussions with calculations based on leads with “floating” Fermi energies.

Without transforming to the gate voltage axis, the original conductance data of Figure 2(e) as a function of energy is reported in Figure 3(a) with a wider energy range up to $\pm 0.4\mathord{\thinspace\rm eV}$. Compared to the density of states [Figure 3(b)] and the band structure [Figure 3(c)] which are calculated based on the same moiré superlattice model potential but within the continuum model (subsection II.2), consistent features in the energy spectrum can be seen. In view of Figure 2(d)–(e) and Figure 3, our calculations significantly capture some of the basic properties of the graphene/hBN moiré superlattice, at least at zero magnetic field.

III.1.2 Transport at finite magnetic field

We continue our comparison of the experimentally measured and theoretically calculated conductance $G(V_{\mathrm{bg}})$ with, however, finite magnetic field $B$ perpendicular to the graphene plane, which can be modeled by associating the Peierls phase Peierls (1933); Datta (1995) to the hopping $t\rightarrow te^{i\Phi}$ in Eq. (1), where $\Phi=(e/\hbar)\int_{\mathbf{r}_{j}}^{\mathbf{r}_{i}}\mathbf{A}\cdot d\mathbf{r}$, choosing the Landau gauge $\mathbf{A}=(-yB,0,0)$ for the vector potential $\mathbf{A}$; see the axes shown in Figure 1. Conductance maps of $G(V_{\mathrm{bg}},B)$ are reported in Figure 4(a)/(b) from the experiment/theory showing magnetic field up to $5\mathord{\thinspace\rm T}$ and the gate voltage range same as Figure 2(d)/(e). Note that the red curves of Figure 2(d)/(e) correspond exactly to the horizontal line cuts at $B=0$ of Figure 4(a)/(b). Within the gate voltage range of about $-45\mathord{\thinspace\rm V}\lesssim V_{\mathrm{bg}}\lesssim+45\mathord{\thinspace\rm V}$, typical relativistic Landau fans can be seen in both experiment and theory. To have a closer look, we magnify the regions marked by the black dashed box on Figure 4(a)/(b) in Figure 4(c)/(d) with a different color map to highlight the quantized conductance plateaus. Numbers $-6,-10,\cdots,-38$ on Figure 4(d) label the filling factor $\nu$ on the corresponding plateau with the expected conductance $|\nu|e^{2}/h$. Good agreement between experiment [Figure 4(a) and (c)] and theory [Figure 4(b) and (d)] within the main Dirac cone can be seen.

At gate voltages $|V_{\mathrm{bg}}|\gtrsim 45\mathord{\thinspace\rm V}$, transport properties are dominated by the extra Dirac cones arising from the modulating moiré superlattice. Discrepancies between the experiment and theory are self-evident. This suggests that the neglected higher-order terms of a more complete model Hamiltonian (see, for example, Refs. Kindermann et al., 2012; Wallbank et al., 2013; Moon and Koshino, 2014) become important when the magnetic field is strong. Interestingly, we note that in the theory map of Figure 4(b), some unusual plateaus in the energy range around the electron-branch secondary Dirac point can be observed. We magnify this region in Figure 4(e) with a horizontal line cut shown in Figure 4(f), where the gate voltage range showing quantized conductance plateaus is highlighted by a yellow background.

This $V_{\mathrm{bg}}$ range, transformed back to the energy through Eq. (10), 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 Figure 3(c). Since there are effectively 3 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\times 2\times 2=12$ with $\times 2$ accounting for spin and another $\times 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 Figure 4(f). 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 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 Ponomarenko et al. (2013); Dean et al. (2013); Lee et al. (2016); Handschin et al. (2017); Krishna Kumar et al. (2017), including ours shown in Figure 2(d), indicating once again that the simplified model of Eq. (8) containing only the electrostatic scalar potential term is not sufficient to capture transport properties of graphene/hBN moiré superlattices at high magnetic fields, i.e. in the quantum Hall regime. As we will see below, when the graphene superlattice potential comes solely from the electrostatic gating, our method becomes exact because in that system no such higher order terms are involved.

III.2 Gate-controlled superlattices

To observe any superlattice effects in graphene, the mean free path must be able to cover enough periods of the superlattice potential. This means, 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 Yankowitz et al. (2012) led to first of its kind studies on graphene superlattices – the periodicity corresponding to small twist angles turns out to be naturally in a suitable range for experiments; see Figure 2(b). On the other hand, the superimposed superlattice potential due to the graphene/hBN moiré pattern is defined as a hexagonal lattice emerging from the two host lattices.

A more flexible approach to design artificial graphene superlattice structures for band structure engineering was pursued with the realization of electrostatic gating schemes Forsythe et al. (2018). 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 Dubey et al. (2013); Drienovsky et al. (2014). 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 Drienovsky et al. (2017). 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 Drienovsky et al. (2017, 2018), or by etching the dielectric layer with a periodic pattern using a standard uniform back gate underneath the modulated substrate Forsythe et al. (2018). The latter will be our focus in the rest of this section.

Following the geometry of the device subject to a square superlattice potential with periodicity $\lambda=35\mathord{\thinspace\rm nm}$ presented in Ref. Forsythe et al., 2018, we have performed our own electrostatic simulation to obtain the back gate capacitance showing periodic spatial modulation. We consider an hBN/graphene/hBN sandwich (showing no measurable moiré superlattice effects) gated by a global top gate contributing a uniform carrier density $n_{\mathrm{tg}}$, and a bottom gate at voltage $V_{\mathrm{bg}}$ with a pre-patterned SiO2 substrate in between. See Fig. 1 of Ref. Forsythe et al., 2018. The bottom gate capacitance therefore shows a spatial modulation with a square lattice symmetry, as shown in the lower left inset of Figure 5(a).

With the electrostatically simulated position-dependent back gate capacitance per unit area $C_{\mathrm{bg}}(x,y)$, contributing carrier density $n_{\mathrm{bg}}(x,y)=[C_{\mathrm{bg}}(x,y)/e]V_{\mathrm{bg}}$, together with the uniform $n_{\mathrm{tg}}$, the resulting superlattice potential is given by $U_{s}(x,y)=-\mathop{\rm sgn}\nolimits[n(x,y)]\hbar v_{F}\sqrt{\pi|n(x,y)|}$ with $n=n_{\mathrm{bg}}+n_{\mathrm{tg}}$, in order to set the global transport Fermi level at zero Liu and Richter (2012). Slightly different from the case of the graphene/hBN moiré superlattice (subsection III.1) where the model potential $U_{s}(x,y)$ given by Eq. (8) is independent of the gating, we consider $U(x,y)=U_{s}(x,y)$ for the on-site energy term (2), and implement it in the tight-binding Hamiltonian (1) with $s_{f}=6$ to perform quantum transport simulations over a two-terminal structure with $L=420\mathord{\thinspace\rm nm}$ and $W=385\mathord{\thinspace\rm nm}$; see the upper right inset of Figure 5(a).

To compare with the resistance measurements reported in Ref. Forsythe et al., 2018, we plot the inverse transmission $1/T$ as a function of $n_{\mathrm{tg}}$ and $V_{\mathrm{bg}}$ in the main panel of Figure 5(a), 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_{\mathrm{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 Figure 5(b)–(d) to show clearly the single- and multiple-peak structures, in excellent agreement with the experiment Forsythe et al. (2018).

We have also checked the consistency between the calculated mini-band structures and the simulated inverse transmission. Overall, we obtain band structures similar to that reported in Ref. Forsythe et al., 2018, but since each $(n_{\mathrm{tg}},V_{\mathrm{bg}})$ point corresponds to a different $U_{s}(x,y)$ profile and hence a different mini-band structure, an overview consistency-check like in Figure 3 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_{\mathrm{tg}},V_{\mathrm{bg}})$ configurations corresponding to the three black arrows in Figure 5(b) 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 (subsection II.2) are shown in Figure 6. Going from low to high $n_{\mathrm{tg}}$ [left, middle, and right black arrow in Figure 5(b)], 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$ [Figure 6(a), (b), and (c)]. As expected, the highest peak in Figure 5(b) marked by the middle black arrow corresponds to Figure 6(b), where the main Dirac point is nearly hit; see the lower sub-panel therein and the relevant caption. From the $E=0$ Fermi contours of Figure 6(a) and (c), the two satellite $1/T$ peaks seen in Figure 5(b) 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 Figure 6, though corresponding to an increasing uniform $n_{\mathrm{tg}}$, do not exhibit simply an energy shift without changing the band shape. Compare, for example, the shapes of the lowest subbands. In addition, in this simulated case, no energy windows accommodating completely isolated extra Dirac points can be found. We note, however, that by properly designing the gate capacitance geometry, it is possible to find isolated extra Dirac points, even at $\Gamma$. When isolated extra Dirac points are found, electronic transport is supported solely by the isolated extra Dirac cone, and more novel transport properties of band-engineered graphene superlattices can be explored. This is beyond the scope of the present work and is left as a future direction to further elaborate.

IV Concluding Remarks

In summary, we have shown that quantum transport simulations based on the scalable tight-binding model Liu et al. (2015) correctly capture transport properties of electrostatic graphene superlattices. In the case of graphene/hBN moiré superlattice (subsection III.1), the consistency of our simulation and experiment at zero and low magnetic field is rather satisfactory but breaks down at strong magnetic field due to the neglected higher-order terms in a more complete model Hamiltonian Kindermann et al. (2012); Wallbank et al. (2013); Moon and Koshino (2014). In the other case of gated superlattices (subsection III.2), without such higher order terms the simulations are expected to be exact. Indeed, compared to the transport experiment with a gate-controlled square superlattice reported in Ref. Forsythe et al., 2018, 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 a future work.

Our work shows that real-space transport simulations based on the scalable tight-binding model Liu et al. (2015) can be extended to treat electrostatic superlattices, whether of the graphene/hBN type or the modulated gate capacitance type, providing consistent results compared to experiments. The method can be immediately applied to take into account, for example, complex local gating or multi-probe transport, in order to make further analysis for transport experiments or even reliable predictions. We note some recent studies working on developing numerical techniques that allow large-scale efficient transport simulations Beconcini et al. (2016); Calogero et al. (2018); Papior et al. (2019), but scaling the graphene lattices with an appropriately chosen scaling factor depending on the superlattice periodicity 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 method Datta (1995) or the popular open-source python package KWANT Groth et al. (2014).

Acknowledgements.

We thank R. Huber, J. Eroms, and C. J. Kent for valuable and stimulating discussions. Financial supports from Taiwan Minister of Science and Technology (grant No. 107-2112-M-006-004-MY3), Deutsche Forschungsgemeinschaft (projects A07 and Ri681/13-1 under SFB 1277), and Helmholtz Society (program STN) are gratefully acknowledged.