Probing miniband structure and Hofstadter butterfly in gated graphene superlattices via magnetotransport

We first simulate the four-point longitudinal resistance $R_{xx}$. We consider a four terminal device shown in the right inset of Figure 1(c), where the system length $L=1152$ nm, width $W=385$ nm, and the top lead width $w=245$ nm. With the four leads labeled in the right inset of Figure 1(c), we calculate $R_{xx}=R_{14,23}$ (for details see Methods) and show its dependence on the top and back gate voltages in the main panel of Figure 1(c). As can be seen from the map, the strength of the backgate mostly influences the superlattice modulation.

Figure 1(d) presents the linecuts of $R_{xx}$ marked in Figure 1(c) with the respective colors. Whereas at $V_{\mathrm{bg}}\approx 0$ only a single Dirac peak is visible, for increasing $|V_{\mathrm{bg}}|$ second and higher order satellite Dirac peaks start to appear, as the periodic modulation gets stronger. At $V_{\mathrm{bg}}=\pm 40$ V [linecuts in Figure 1(d)], a few higher-order Dirac points are resolved. On the other hand, changing the top gate voltage mostly tunes the carrier density in the device. The left inset of Figure 1(c) shows a close-up of the boxed region with $30\ \mathrm{V}\leq V_{\mathrm{bg}}\leq 45$ V and $-3.5\ \mathrm{V}\leq V_{\mathrm{tg}}\leq-1.6$ V, where two sharp lines are visible at both sides of the main Dirac peak, corresponding to the secondary Dirac points, and several fainter lines, corresponding to higher order Dirac points. Similar results based on the same capacitance model function have been reported in Chen2020 , where two-terminal transport simulations were performed. The character of the bands can be verified in magnetotransport, as shown in the following subsections.

For a general band dispersion $\epsilon(\mathbf{k})$, the semiclassical equations of motion for an electron are given by Ashcroft1976

	$\displaystyle\mathbf{\dot{r}}$	$\displaystyle=\frac{1}{\hbar}\nabla_{\mathbf{k}}\epsilon(\mathbf{k}),$		(1)
	$\displaystyle\hbar\mathbf{\dot{k}}$	$\displaystyle=-e(\mathbf{E}+\mathbf{\dot{r}}\times\mathbf{B}),$		(2)

where $-e$ is the electron charge, $\mathbf{E}$ is the external electric field, and $\mathbf{B}$ the magnetic field. In the presence of constant out-of-plane magnetic field $\mathbf{B}=(0,0,B)$ only, one can obtain the relation between the shape of the Fermi contour in the momentum space $\Delta\mathbf{k}(t)$ and the carrier trajectory in real space $\Delta\mathbf{r}(t)$ Ashcroft1976

$\displaystyle\Delta\mathbf{r}(t)=\frac{\hbar}{eB}\Delta\mathbf{k}(t)\times\hat{z},$

(3)

meaning that the cyclotron orbit is obtained by rotating the orbit in the momentum space by $90^{\circ}$ clockwise, as illustrated in Figure 2(b)–Figure 2(c). Carriers encircle closed orbits of electron type or hole type in the counterclockwise or clockwise direction, respectively, as determined via the group velocity, $\mathbf{v}=\nabla_{\mathbf{k}}\epsilon(\mathbf{k})/\hbar$, and the equation $\hbar\mathbf{\dot{k}}=-eB\mathbf{v}\times\hat{z}$. In pristine graphene at low energy, cyclotron orbits exhibit a circular shape [Figure 2(b)], but the bands formed in systems with SL modulation can be highly distorted from the original, conical shape, and thus, noncircular Fermi contours can be observed [e.g. Figure 2(c)].

In a typical device designed for transverse magnetic focusing measurement, an emitter and collector [in Figure 2(a) marked as 1 and 2, respectively] are located on the same edge at a center to center distance $\ell=d+w$, and other contacts act as absorbers. The current injected from an emitter flows along cyclotron orbits with a radius $R_{\mathrm{c}}$ depending on the magnetic field strength, and at the boundary propagates along skipping orbits. The current can end up in the collector, when the diameter or its multiples match $\ell$ [Figure 2(b)], or otherwise in absorber contacts. In the nonlocal resistance measurement [Figure 2(a)], this results in maximum or minimum, respectively, of the resistance $R_{\mathrm{f}}=R_{14,23}$ (see Methods). For electron-like (hole-like) orbits, the resistance maximum condition can be obtained for positive (negative) $B$.

In pristine graphene, a typical TMF spectrum as a function of magnetic field and voltage contains two fans of focusing peaks, one for the electron band and the other for the hole band Taychatanapat2013 . In a superlattice, the emerging replicas of the Dirac cone cause a substantial modification of the TMF signal. In the following discussion, we choose $V_{\mathrm{bg}}=40$ V, such that a few higher-order Dirac peaks are present next to the main Dirac point as seen in Figure 1(c).

For transverse magnetic focusing, we choose the system geometry shown in Figure 2(a), with the distance between the bottom leads $d=1200$ nm, their widths $w=100$ nm, the side leads width $W=1636$ nm, and the top lead width $L=1792$ nm.

Figure 3(a) shows the $R_{\mathrm{f}}$ map as a function of $B$ and $V_{\mathrm{tg}}$. One can see several series of fans, with the focusing peaks appearing at $B>0$ or $B<0$. The map is put together with the $R_{xx}$ calculated at $B=0$, shown in Figure 3(b), where the main and secondary Dirac peaks are seen. The sign change of the focusing peaks in Figure 3(a) occurs either at the Dirac points, or van Hove singularities, and for the former, it coincides with the $R_{xx}$ peaks. This confirms the sign change of the carriers when tuning $V_{\mathrm{tg}}$, occurring as an effect of the band reconstruction due to the superlattice potential. To understand the result in detail, in Figure 3(c)–Figure 3(f) we plot the miniband structures calculated at $B=0$, as described in Ref. Chen2020 , and the Fermi contours at $E=0$, at selected values of $V_{\mathrm{tg}}$ marked by arrows in Figure 3(a). Additionally, in Figure 3(g)–Figure 3(k) we show zoomed regions of the $R_{\mathrm{f}}$ map marked with the colored rectangles in Figure 3(a).

In Figure 3(c), at $V_{\mathrm{tg}}=-2.75$ V, the Fermi level is located at the $C_{1}$ subband, and the Fermi contour has a rounded shape [Figure 2(b) and Figure 3(c)]. In the semiclassical description, electrons injected from lead 1 [Figure 2(a)] with an initial velocity $\mathbf{v}=(0,v)$ and $\mathbf{k}=(0,k)$, in a moderate magnetic field travel along a rounded trajectory, and after half a period, $t=T/2$ encircle half of the closed orbit. From Eq. (3), this corresponds to traveling a distance equal to the diameter, $2R_{\mathrm{c}}=2\hbar k/eB$ along the $x$ direction [Figure 2(b)]. For $\ell=2R_{\mathrm{c}}$, it leads to the first maximum of $R_{\mathrm{f}}$. For smaller $R_{\mathrm{c}}$, the beam is reflected at the edge, and can flow to the collector when $\ell=4R_{\mathrm{c}},6R_{\mathrm{c}},\dots$, giving rise to higher order $R_{\mathrm{f}}$ peaks. In general, we can evaluate the field at which the $j$th maximum occurs as $B_{j}=2\hbar jk/e\ell$, $j=1,2,\dots$. We find $k(V_{\mathrm{tg}})$ numerically and plot $B_{j}$ with dashed lines in Figure 3(g). The $C_{1}$ subband cone is within $-2.9\ \mathrm{V}\lesssim V_{\mathrm{tg}}\lesssim-2.5$ V. We find a very good agreement with the $R_{\mathrm{f}}$ signal for up to $j\approx 7$. Higher $j$ are not resolved as the system enters the quantum Hall regime, and semiclassical description of the skipping orbits at the edge no longer applies. At $-3.1\ \mathrm{V}\lesssim V_{\mathrm{tg}}\lesssim-2.9$ V, for the $V_{1}$ subband, $R_{\mathrm{f}}$ is noisy due to scattering of low-energy carriers by the periodic potential.

When the Fermi level is tuned to the van Hove singularity (at $V_{\mathrm{tg}}\approx-2.5$ for the electron subband, and $V_{\mathrm{tg}}\approx-3.1$ for the hole subband), the focusing signal vanishes, and smaller fans reappear. Based on the miniband structures [Figure 3(d)], we interpret them as due to focusing of the secondary Dirac cones fermions. For example, at $V_{\mathrm{tg}}=-3.12$ V [Figure 3(d)] at $E=0$ there are tiny Dirac cones around the $M$ and $X$ points of the Brillouin zone.

For $V_{\mathrm{tg}}\lesssim-3.2$ V and $V_{\mathrm{tg}}\gtrsim-2.4$ V, the miniband structures around the Fermi level get more complex, with many overlapping subbands. Nevertheless, we find ranges of $V_{\mathrm{tg}}$ where a single isolated higher-order Dirac cone is present, giving clear focusing signal [see Figure 3(e) for the $V_{2}$ subband, the corresponding $R_{\mathrm{f}}$ zoom in Figure 3(h), and the zoom in Figure 3(i) for the $C_{2}$ subband]. When there are more overlapping subbands, the signal gets very faint. Nevertheless, one can spot fans that fit well to the hole-like orbit within the $V_{4}$ subband around the $M$ point, see Figure 3(f) at $V_{\mathrm{tg}}=-3.85$ V, and zoomed $R_{\mathrm{f}}$ in Figure 3(j). A similar feature is resolved in Figure 3(k) for the electron-like orbits.

To further illustrate the relation of the real-space orbits to the subbands, we calculate the current density maps at selected focusing peaks. Figure 4(a)–Figure 4(c) show the line cuts of $R_{\mathrm{f}}$ at $V_{\mathrm{tg}}$ and $B$ marked by the respective colors in Figure 3(g), Figure 3(h) and Figure 3(j). In Figure 4(a) for the $C_{1}$ subband, the first two focusing peaks are marked by $\circ$ and $\triangleleft$. The respective current density maps are presented in Figure 4(d) and Figure 4(e), revealing typical current densities found for TMF calculations in pristine graphene Stegmann2015 ; Petrovic2017 . In Figure 4(b) the line cut for the $V_{2}$ subband is shown. For the focusing peaks marked with $\diamond$ and $\triangledown$, the current density maps are shown in Figure 4(f) and Figure 4(g). The trajectories acquire a shape close to a rectangle, consistent with the Fermi contour [Figure 3(e)]. For the line cut in the hole-like $V_{4}$ subband (purple) [Figure 4(c)], the current densities of the first two peaks are shown in Figure 4(h) and Figure 4(i)]. The orbits acquire a rhombus shape, matching well the Fermi contour in the corner of the Brillouin zone [Figure 2(c) and Figure 3(f)]. The red dashed lines show semiclassical trajectories calculated using the Fermi contours obtained from our band structures and Equation 3 with the contour starting at $\mathbf{k}$ for which $v_{x}\propto\partial E/\partial k_{x}=0$. These semiclassical orbits show similarity with the current density, but the current density is obtained from quantum calculation so they are expected to be similar but not strictly identical, in particular, counter-intuitive patterns may occur at certain resonant conditions (e.g. a vertical blob on top of the rhombus-like pattern in Figure 4(h)). Note that in Figure 4(f)–Figure 4(g) the current density is non-zero in the area between the vertical segments, as it contains contributions from multiple initial $k_{x}$ for which $v_{x}$ is low. The noisy background visible in Figure 4(a)–Figure 4(c) originates from scattering and resonant states due to SL which are irregular and complex due to the wave-like nature of the carriers.

Although the focusing spectrum is not symmetric with respect to the main Dirac point, it is qualitatively similar in minibands above and below the main Dirac point, except for the noisy signal for low-energy valence subbands. Let us note that the modulation induced by electrostatic gates is a complex function of $V_{\mathrm{bg}}$ and $V_{\mathrm{tg}}$, and the band structure is significantly modified merely by changing $V_{\mathrm{tg}}$ [Figure 3(c)–Figure 3(f)]. This is in contrast to the SLs induced via low-angle twisted hBN substrate or twisted graphene layers, where the top gate only sweeps the carrier density without affecting the periodic modulation, and the band structure remains unchanged. This leads to the shape of the fans in $R_{\mathrm{f}}$ closer to straight lines, unlike Refs. Lee2016 ; Berdyugin2020 that observed ones close to parabolic.

Now, we turn our attention to the high-field regime, where we expect the Hofstadter spectrum to emerge when the magnetic flux per SL unit cell of area $A$, $\phi=AB$, is of the order of the flux quantum $\phi_{0}=h/e$. To simulate longitudinal resistance $R_{xx}$ and Hall resistance $R_{xy}$ at the same time, while keeping calculations for the four-point resistances minimized, we consider a 5-terminal Hall bar sketched in Figure 5(a) with the actually considered geometric dimensions shown. We compute all the $5\times 4=20$ transmission functions between all pairs among the five leads, and process the data to obtain longitudinal resistance $R_{xx}=R_{14,23}$ and Hall resistance $R_{xy}=R_{14,52}$ following the Büttiker formalism Buttiker1986 ; Datta1995 , according to the lead labels shown in Figure 5(a). For all the following discussions, we convert our $V_{\mathrm{tg}}$ axis to the numerically obtained average carrier density $\bar{n}$ over the entire lattice, in order for a more transparent presentation of our high-field transport simulations.

Figure 5(b) and Figure 5(c) show the Hall conductance $G_{\mathrm{H}}=R_{xy}^{-1}$ as a function of $\bar{n}$ and magnetic field $B$, with the back gate voltage fixed at $V_{\mathrm{bg}}=0$ and $V_{\mathrm{bg}}=40\mathord{\thinspace\rm V}$, respectively. To better highlight the Landau fans arising from the quantum Hall effect of graphene, we limit the color range to $-40\leq G_{\mathrm{H}}/G_{0}\leq+40$ in Figure 5(b), where $G_{0}=e^{2}/h$ is the conductance quantum. The same limit of the color range is applied to Figure 5(c) in order for a direct comparison. Line cuts at $B=2\mathord{\thinspace\rm T}$ marked by the black line on Figure 5(b) and orange line on Figure 5(c) are shown and compared in Figure 5(d). The lowest few quantum Hall conductance plateaus $G_{\mathrm{H}}=\pm 2G_{0},\pm 6G_{0},\pm 10G_{0},\cdots$ can be clearly seen in the $V_{\mathrm{bg}}=0$ case free of superlattice potential (black line), while the combined strong magnetic field ($B=2\mathord{\thinspace\rm T}$) and strong superlattice modulation ($V_{\mathrm{bg}}=40\mathord{\thinspace\rm V}$) result in a more complex transport feature, in accordance with the predictions Thouless1982 ; Streda1982 for a periodic 2D modulation, which can be better understood by checking $R_{xx}$, instead of $R_{xy}$ (or $G_{\mathrm{H}}$), to be discussed soon below.

Without taking the inverse of $R_{xy}$, Figure 5(e) shows the Hall resistance at $B=0.2\mathord{\thinspace\rm T}$ with $V_{\mathrm{bg}}=0$ [red, corresponding to the $\bar{n}$ range marked on Figure 5(b)] and $V_{\mathrm{bg}}=40\mathord{\thinspace\rm V}$ [cyan, corresponding to the $\bar{n}$ range marked on Figure 5(c)]. The former (without superlattice) shows a single sign change at $\bar{n}=0$, typical for graphene, while the latter (with superlattice) shows multiple sign changes at positive and negative $\bar{n}$, in addition to the main Dirac point at $\bar{n}=0$, consistent with our previous low-field magnetotransport simulations discussed above.

Considering the same range of $\bar{n}$ and $B$ as Figure 5(b) and Figure 5(c), Figure 5(f) shows the longitudinal resistance $R_{xx}$ with the back gate voltage fixed at $V_{\mathrm{bg}}=40\mathord{\thinspace\rm V}$. Since the period of our gate-controlled superlattice is $\lambda=35$ nm, and hence the square SL unit cell area $A=\lambda^{2}=1225\mathord{\thinspace\rm nm^{2}}$, the condition $\phi/\phi_{0}=AB/(h/e)=1$ is reached at $B\approx 3.376\mathord{\thinspace\rm T}\equiv B_{1}$, which is in sharp contrast with the graphene/hBN moiré superlattice that requires $B\approx 24\mathord{\thinspace\rm T}$ in order to reach $\phi/\phi_{0}=1$, because of its periodicity limited to $\lambda\lesssim 14\mathord{\thinspace\rm nm}$ and hence the area $A=\sqrt{3}\lambda^{2}/2\lesssim 170\mathord{\thinspace\rm nm^{2}}$. As is visible on Figure 5(f), the seemingly complicated $R_{xx}$ map does exhibit certain features at $B_{1}$ and $B_{1}/2$, and vaguely at $B_{1}/3$ (marked by black triangles), corresponding to $\phi/\phi_{0}=1,1/2,1/3$, respectively, consistent with our calculation of the magnetic energy subbands shown in Figure 5(g), where the vertical axis of $\phi/\phi_{0}$ corresponds to exactly the same $B$ range as Figure 5(f), as well as the average density range. Good consistency between the $R_{xx}$ map of Figure 5(f) and the Hofstadter butterfly shown in Figure 5(g) can be seen. For methods adopted to calculate Figure 5(g), see Methods.

To model a realistic geometry, we consider SL induced by a patterned dielectric substrate Forsythe2018 ; Li2021 with a uniform global backgate underneath the dielectric layer, and the graphene sheet sandwiched between two hBN layers lying on top [Figure 1(a) of the main text]. The hBN/graphene/hBN sandwich is covered by a global top gate. The voltage applied to the back gate $V_{\mathrm{bg}}$ controls the strength of the periodic modulation, while the top gate voltage $V_{\mathrm{tg}}$ is used to tune the carrier density across the SL. We follow the design of Ref. Forsythe2018 , with a square lattice etched in SiO₂ substrate, with a lattice period $\lambda=35$ nm. We use a model function $C_{\mathrm{bg}}(x,y)$

	$\displaystyle C_{\mathrm{d}}(x,y)$	$\displaystyle=C_{\mathrm{out}}-(C_{\mathrm{out}}-C_{\mathrm{in}})\tanh\left[\exp\left(-\frac{x^{2}+y^{2}}{d_{\mathrm{smooth}}^{2}}\right)\right],$		(4)
	$\displaystyle C_{\mathrm{bg}}(x,y)$	$\displaystyle=C_{\mathrm{d}}\left[\left(x-\frac{\lambda}{2}\right)\%\lambda-\frac{\lambda}{2},\left(y-\frac{\lambda}{2}\right)\%\lambda-\frac{\lambda}{2}\right],$		(5)

where $d_{\mathrm{smooth}}=7.5$ nm is the smoothness of the modulation, $C_{\mathrm{out}}/e=0.77\times 10^{11}\ \mathrm{cm}^{-2}\mathrm{V}^{-1}$ and $C_{\mathrm{in}}/e=0.22\times 10^{11}\ \mathrm{cm}^{-2}\mathrm{V}^{-1}$ , and $a\%b$ means the remainder after dividing $a$ by $b$. The top gate capacitance is assumed to be $C_{\mathrm{tg}}/e=1\times 10^{12}\ \mathrm{cm}^{-2}\mathrm{V}^{-1}$. Using the parallel capacitor model, this roughly corresponds to the top hBN thickness $d_{\mathrm{t}}\approx 16.6$ nm, from $C_{\mathrm{tg}}/e=\varepsilon_{0}\varepsilon_{\mathrm{hBN}}/ed_{\mathrm{t}}$, where $\varepsilon_{0}$ is the vacuum permittivity, $\varepsilon_{\mathrm{hBN}}=3$ is the dielectric constant of hBN, and $-e$ is the electron charge. Figure 1(b) of the main text shows the profile of the above model function (5).

For the dual-gated graphene sample free from intrinsic doping, the carrier density is given by

$$n(x,y)=\frac{C_{\mathrm{bg}}(x,y)}{e}V_{\mathrm{bg}}+\frac{C_{\mathrm{tg}}}{e}V_{\mathrm{tg}}.$$

(6)

Assuming that the carrier energy in graphene is given by $E=\pm\hbar v_{\mathrm{F}}k$, where $\hbar$ is the reduced Planck constant, $v_{\mathrm{F}}\approx 10^{6}\mathord{\thinspace\rm m/s}$ is the Fermi velocity of graphene, and using $\hbar v_{\mathrm{F}}\approx 3\sqrt{3}/8\ \mathord{\thinspace\rm eVnm}$, the onsite potential energy can be obtained from

$$U=-\text{sgn}(n)\hbar v_{\mathrm{F}}\sqrt{\uppi|n|},$$

(7)

in order to set the global Fermi energy at $E=0$ where transport occurs.

For transport calculation, we use the tight-binding Hamiltonian

$$H=-\sum\limits_{\left\langle{i,j}\right\rangle}t_{ij}c_{i}^{\dagger}c_{j}+\sum\limits_{j}{U({\mathbf{r}}_{j})}c_{j}^{\dagger}c_{j},$$

(8)

where $c_{i}$ ($c_{i}^{\dagger}$) is an annihilation (creation) operator of an electron on site $i$ located at $\mathbf{r}_{i}=(x_{i},y_{i})$. The first sum contains the nearest-neighbor hoppings with the hopping parameter $t_{ij}$, and the second sum describes the onsite potential energy profile. In the presence of an external magnetic field $\mathbf{B}=(0,0,B)$, the hopping integral is modified by $t_{ij}\rightarrow t_{ij}\exp(i\phi)$, where the Peierls phase $\phi=(-e/\hbar)\int_{\textbf{r}_{i}}^{\textbf{r}_{j}}\textbf{A}\cdot d\textbf{r}$, with A being the vector potential that satisfies $\nabla\times\mathbf{A}=\mathbf{B}$, and the integral going from the site at $\textbf{r}_{i}$ to the site at $\textbf{r}_{j}$. For a feasible simulation of realistic devices, we use the scalable tight-binding model Liu2015 , where the hopping parameter becomes $t_{ij}=t_{0}/s_{\mathrm{f}}$, and the lattice spacing $a=a_{0}s_{\mathrm{f}}$, $s_{\mathrm{f}}$ is the scaling factor, and we use $t_{0}=-3$ eV and $a_{0}=1/4\sqrt{3}$ nm. Transport calculations based on Hamiltonian (8) are done within the wave-function matching for the TMF, and real-space Green’s function method in other cases, at the global Fermi energy $E=0$ and zero temperature. The conductance from lead $i$ to lead $j$ is obtained from the Landauer formula $G_{ji}=2(e^{2}/h)T_{ji}$, where the transmission probability $T_{ji}$ is evaluated as a sum over the propagating modes $T_{ji}=\sum\limits_{q}T_{ji}^{q}$, and

$$T_{ji}^{q}=\sum\limits_{p}|t_{ij}^{pq}|.$$

(9)

Here, $t_{ij}^{pq}$ is the probability amplitude for the transfer from the incoming mode $p$ in lead $i$ to the outgoing mode $q$ in lead $j$.

In the multiterminal devices, we solve the transport problem for each lead as an input, and build the conductance matrix $\mathbfcal{G}$ Buttiker1986 ; Datta1995 which relates the current $I_{i}$ fed to the system in lead $i$ to the voltage at $j$-th lead $V_{j}$ through $I_{i}=\sum_{j=1}^{N}\mathcal{G}_{ij}V_{j}$. For an $N$-terminal system, the matrix elements are

	$\displaystyle\mathcal{G}_{ij}$	$\displaystyle=-G_{ij},\ \quad i\neq j,$		(10)
	$\displaystyle\mathcal{G}_{ii}$	$\displaystyle=\sum\limits_{j=1,j\neq i}^{N}G_{ij}.$		(11)

We set the voltage at $l$-th lead equal to zero, and eliminate the $l$-th row and column of the matrix. The reduced $(N-1)\times(N-1)$ matrix $\bar{\mathbfcal{G}}$ can be inverted to get $\boldsymbol{\mathbfcal{R}}=\bar{\mathbfcal{G}}^{-1}$, where the $\mathbfcal{R}$ matrix satisfies

$$V_{i}=\sum_{j=1,j\neq l}^{N}\mathcal{R}_{ij}I_{j}.$$

(12)

With the elements of matrix $\boldsymbol{\mathbfcal{R}}$, one can evaluate the resistance

$$R_{kl,mn}=\frac{V_{m}-V_{n}}{I_{k}}=\mathcal{R}_{mk}-\mathcal{R}_{nk}$$

(13)

with the current flowing from lead $k$ to lead $l$, zero current in other terminals, and voltage drop measured between leads $m$ and $n$.

As a numerical example of applying the above outlined Landauer-Büttiker formalism for computing the four-point resistance Eq. (13), we revisit the first TMF experiment on graphene Taychatanapat2013 , considering the same probe spacing ($\ell=500$ nm) and width ($100$ nm) but slightly different geometry of the scattering region (total length $700$ nm and width $500$ nm) for a 6-terminal Hall bar, made of a graphene lattice scaled by $s_{\mathrm{f}}=12$. Choosing the same configuration of the leads for injector, ground, and voltage probes as the revisited experiment, the computed $R_{61,54}$ as a function of the external magnetic field $B$ perpendicular to the plane of graphene and the carrier density $n$ is reported in Figure 6(a), showing a map rather consistent with the experiment. Due to the isotropic low-energy dispersion of graphene, the resulting cyclotron trajector is a simple circle of radius $R_{\mathrm{c}}=\hbar k/eB$, which is simplified from Eq. (3). By requiring the probe spacing to be equal to an integer times the cyclotron diameter, $\ell=j\cdot 2R_{\mathrm{c}}$, $j=1,2,\cdots$, together with $k=\sqrt{\uppi|n|}$ for graphene, one can solve for carrier density corresponding to the $j$th peak of the TMF on the $B$-$n$ map of Figure 6(a):

$$n_{j}(B)=\frac{1}{\uppi}\left(\frac{eB\ell}{2\hbar j}\right)^{2}\ .$$

(14)

The dashed lines on Figure 6(a) show $n_{1},n_{2},n_{3},n_{4}$, matching very well with the patterns of the simulated $R_{61,54}$, which requires totally $6\times 5=30$ transmission functions for such a 6-terminal device, as explained above. Figures 6(b) and (c) show two exemplary maps of transmission functions, which can look generally very different from the resulting four-point resistance.

In the presence of an external magnetic field and semi-infinite leads, the vector potential must satisfy the translational invariance of the leads. For the magnetic field along the $\hat{z}$ axis, the most common choice is the Landau gauge, $\mathbf{A}=(-yB,0,0)$ or $\mathbf{A}=(0,xB,0)$ for the lead which is translationally invariant along the $x$ or $y$ direction, respectively. For other lead orientation, in general, the proper gauge is different. Therefore, in a multi-terminal device, the required vector potential is not uniform in the entire space. This is not a problem since adding an arbitrary curl-free component to the vector potential does not change the magnetic field. Here, we use the approach introduced in Baranger1989 .

Assuming that the proper gauge in the 1st lead is $\mathbf{A}_{1}(\mathbf{r})$, for another lead that is at an angle $\theta_{n}$ with respect to lead 1, the gauge can be chosen as

$$\mathbf{A}_{n}(\mathbf{r})=\mathbf{A}_{1}(\mathbf{r})+\boldsymbol{\nabla}f_{n}(\mathbf{r}),$$

(15)

with

$$f_{n}(\mathbf{r})=Bxy\sin^{2}(\theta_{n})+\frac{1}{4}B(x^{2}-y^{2})\sin(2\theta_{n}).$$

(16)

The addition of a gradient of a scalar function does not influence the requirement $\nabla\times\mathbf{A}=\mathbf{B}$. As an illustration of the transformation, consider $\mathbf{A}_{1}(\mathbf{r})=(-yB,0,0)$, and $\theta_{n}=90^{\circ}$. Then, $f_{n}(\mathbf{r})=Bxy$, $\nabla f_{n}(\mathbf{r})=(By,Bx,0)$, and $\mathbf{A}_{n}(\mathbf{r})=(-yB,0,0)+(yB,xB,0)=(0,xB,0)$.

Applying the transformation (15) so that it only affects lead $n$ is possible by defining a smooth step function $\zeta_{n}(\mathbf{r})$ which is only nonzero in the translationally invariant part of lead $n$

$$\zeta_{n}(\mathbf{r})=\begin{cases}1,&\mathbf{r}\text{ in lead }n,\\ 0,&\mathbf{r}\text{ in lead }m\neq n,\\ \text{smooth interpolation}&\mathbf{r}\text{ elsewhere}.\end{cases}$$

(17)

Then, in (15) we substitute $f_{n}(\mathbf{r})\rightarrow\zeta_{n}(\mathbf{r})f_{n}(\mathbf{r})$ for lead $n$. In general, for the entire system we define

$$f(\mathbf{r})=\sum\limits_{n=2}^{N}\zeta_{n}(\mathbf{r})f_{n}(\mathbf{r}),$$

(18)

this completes our gauge transformation. Adopting the vector potential

$$\mathbf{A}(\mathbf{r})=\mathbf{A}_{1}(\mathbf{r})+\boldsymbol{\nabla}f(\mathbf{r}),$$

(19)

we have $\mathbf{B}=\nabla\times\mathbf{A}$ everywhere in the system, and the translation invariance in each lead is guaranteed. Importantly, curl of (19) gives exactly the desired $B$, regardless of the smoothness of the $\zeta_{n}$ function.

As an example, for the system used for the TMF modeling [Figure 2(a)], in the vertical leads we choose the same gauge $\textbf{A}=(0,xB)$ with $\zeta(\textbf{r})=(\exp(-(y-y_{1})/d_{\mathrm{step}})+1)^{-1}+(\exp(-(y_{2}-y)/d_{\mathrm{step}})+1)^{-1}$, $y_{1}=1607$ nm, $y_{2}=-49$ nm, and $d_{\mathrm{step}}=2$ nm. In the rest of the system, $\textbf{A}=(-yB,0)$ is used. The resulting $A_{x}$ and $A_{y}$ profiles are shown in Figure 7.

For the calculation of Hofstadter butterfly one has to consider a magnetic unit cell whose length is the least common multiple of the lattice periodicity and the periodicity introduced by the Peierls phase. For graphene, it contains more than hundreds of thousands of carbon atoms when the magnetic field strength is smaller than 1 T. However, the calculation is greatly simplified by considering a graphene ribbon. For an armchair ribbon with translational invariance along the $x$ direction and finite width along the $y$ direction, in the presence of a perpendicular magnetic field, dispersionless Landau levels appear near $k_{x}=0$, and the dispersive edge states show up at larger $k_{x}$. Calculating $E_{k_{x}=0}$ as a function of magnetic field, we get the Hofstadter butterfly of graphene. Because of the finite width of the ribbon, the spectrum can also contain edge states. With the increase of $B$ the Landau levels elongate, and at some $B$ the edge states are pushed to $k_{x}=0$, which results in the appearance of the states in the gaps. To lower the computational burden, we use the scalable tight-binding model Liu2015 with $s_{\mathrm{f}}\sim 8$ to calculate $E_{k_{x}=0}$ as a function of magnetic field for an armchair ribbon with periodic length along the $x$ axis equal to the superlattice period ($\lambda$). Here, in order to ensure superlattice period equal to a multiple of $3a$, $s_{\mathrm{f}}$ is not an integer, and the ribbon width is larger than $20\lambda$ to show the superlattice effect.

In transport, only the states at the Fermi level contribute to the conductance. Therefore, we calculate Hofstadter butterfly spectra for all $V_{\mathrm{tg}}$ values and collect the Fermi states to construct the gate-dependent Hofstadter butterfly spectrum to compare with $R_{xx}$ and $G_{\mathrm{H}}$ obtained from the transport calculation. Note that the spectrum in Figure 5(g) contains energy levels that appear across the gaps, which are an artifact of the method due to the finite width of the ribbon. As mentioned above, they appear since at some value of $B$ the edge states are pushed to $k_{x}=0$.

\bmhead

Acknowledgments We thank National Science and Technology Council of Taiwan (grant numbers: MOST 109-2112-M-006-020-MY3 and NSTC 112-2112-M-006-019-MY3) for financial supports and National Center for High-performance Computing (NCHC) for providing computational and storage resources. This research was supported in part by PL-Grid Infrastructure, and by the program ,,Excellence Initiative – Research University” for the AGH University of Science and Technology.