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

Abstract

The presence of periodic modulation in graphene leads to a reconstruction of the band structure and formation of minibands. In an external uniform magnetic field, a fractal energy spectrum called Hofstadter butterfly is formed. Particularly interesting in this regard are superlattices with tunable modulation strength, such as electrostatically induced ones in graphene. We perform quantum transport modeling in gate-induced square two-dimensional superlattice in graphene and investigate the relation to the details of the band structure. At low magnetic field the dynamics of carriers reflects the semi-classical orbits which depend on the mini band structure. We theoretically model transverse magnetic focusing, a ballistic transport technique by means of which we investigate the minibands, their extent and carrier type. We find a good agreement between the focusing spectra and the mini band structures obtained from the continuum model, proving usefulness of this technique. At high magnetic field the calculated four-probe resistance fit the Hofstadter butterfly spectrum obtained for our superlattice. Our quantum transport modeling provides an insight into the mini band structures, and can be applied to other superlattice geometries.

Introduction

Graphene, a 2D material characterized by a linear low-energy dispersion relation, hosts charge carriers named Dirac fermions due to the resemblance of relativistic (massless) particles described by the Dirac equation. Modifying the underlying graphene lattice by a smooth periodic potential can affect the band structure through folding of the pristine graphene Dirac cone into mini bands¹, formation of the secondary Dirac points, and anisotropic renormalization of velocity^2,3,4,5,6. Periodic modulation has been obtained in graphene through chemical functionalization⁷, placing graphene on self-assembled nanostructures⁸, and by stacking graphene together with aligned or slightly misoriented hexagonal boron nitride (hBN), resulting in periodic moiré modulation which generates hexagonal superlattices (SLs)^{9,10,11,12,13,14,15}. Moiré SLs were also created in low-angle twisted graphene bilayers, followed by van der Waals structures made up of few layers of graphene^{16,17,18,19,20,21} and other 2D materials^22,23. SLs are suitable for the observation of the self-similar energy spectrum called Hofstadter butterfly^13,14 and Brown-Zak oscillations^24,25 that occur when the magnetic flux through the superlattice unit cell is of the order of the magnetic flux quantum, ϕ₀ = h/e, and in pristine 2D crystals require unattainable magnetic fields. Also worth mentioning are the many-body phenomena present in moiré SLs^26,27.

Despite their potential, artificial lattices tailored by chemical methods or moiré SLs suffer from inability to tune the strength of the periodic potential. In moiré SLs the period can be tuned to some extent via the rotation angle between the stacked layers, but they are inherent of a hexagonal symmetry. Precise control over the SL geometry, period, and strength is vital for the band structure engineering. The above limitations can be circumvented in electrostatic gate-induced SLs that allow an arbitrary design via the gates geometry, with the gate voltage being a knob for the potential strength. The experimental attempts to create gated SLs in graphene included 1D arrays of metal gates^28,29,30, followed by patterned dielectric substrates^31,32, and few-layer graphene patterned bottom gates for 1D^33,34 and 2D SLs^35,36 with down to sub-20 nm periods³⁴ manifesting the flexibility of this approach. With the recent advance in the fabrication techniques, gated SLs with the period of a few tens of nanometers are achievable with good device quality and long electron mean free path³². This enabled observation of commensurability oscillations^32,33,36, Hofstadter butterfly^34,35 and Brown-Zak oscillations³⁶ at magnetic fields of the order of a few tesla. This is more affordable compared to about 25 T required for graphene/hBN SL, where, due to small lattice constant mismatch of 1.8%, the SL period can reach up to 14 nm for aligned lattices.

While transport in quantizing magnetic field in gated SLs has been thoroughly studied, the intermediate magnetic field regime remained mostly unexplored. In the semiclassical treatment, at low magnetic field fermions undergo cyclotron motion that can be probed in transport measurements via transverse magnetic focusing (TMF). This technique has been used to experimentally probe the band structure in pristine mono-, bi-, and tri-layer graphene³⁷, graphene/WSe₂ heterostructures³⁸, and moiré superlattices^39,40, and theoretically considered in graphene pn junctions^41,42 and 1D gated SLs⁶.

In this work, we perform a theoretical study of the TMF in 2D gate-induced square SLs, and analyze the relation between the observed TMF spectra and the miniband structure. This is complemented by the investigation of magnetotransport in the quantum Hall regime, where we observe signatures of Hofstadter butterfly, matching the numerical Hofstadter spectrum calculated for our gated SL. Our study is performed using quantum transport calculations for multiterminal structures, considering realistic experimental conditions. To model a realistic geometry, we consider SL induced by a patterned dielectric substrate^31,32 with a uniform global backgate underneath the dielectric layer, and the graphene sheet sandwiched between two hBN layers lying on top (Fig. 1a). The hBN/graphene/hBN sandwich is covered by a global top gate. The voltage applied to the back gate V_bg controls the strength of the periodic modulation, while the top gate voltage V_tg is used to tune the carrier density across the SL. We follow the design of Ref. ³¹, with a square lattice etched in SiO₂ substrate, with a lattice period λ = 35 nm (Fig. 1b). For the ease of the calculation, we use a model function C_bg(x, y) (see Fig. 1b) that approximates the electrostatically simulated capacitance, obtained previously by some of us⁴³. Previous studies^35,44 and modeling of previous experiment on magnetic focusing in monolayer graphene³⁷ show good agreement between experiment and simulation for graphene superlattice devices (see Methods), and the present purely theoretical work can be regarded as a guide for further experimental magneto-transport studies.

Fig. 1: Zero magnetic field characterization of the superlattice.

a Sketch of the gated superlattice system. The backgate (dark gray) lays below the patterned substrate (light gray), which induces modulated potential in graphene sandwiched between two hBN layers (blue). The global top gate is shown in yellow. b The backgate capacitance profile with the periodicity λ = 35 nm in the units of 10¹²cm⁻²V⁻¹. c Longitudinal resistance R_xx as a function of the backgate and top gate voltages in a four-terminal device shown in the right inset. Left inset: zoom of R_xx in the boxed region. d R_xx line cuts as marked with the respective line colors in (c).

Results

No external magnetic field

We first simulate the four-point longitudinal resistance R_xx. We consider a four terminal device shown in the right inset of Fig. 1c, 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 Fig. 1c, 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 Fig. 1c. As can be seen from the map, the strength of the backgate mostly influences the superlattice modulation.

Figure 1d presents the linecuts of R_xx marked in Fig. 1c with the respective colors. Whereas at V_bg ≈ 0 only a single Dirac peak is visible, for increasing ∣V_bg∣ second and higher order satellite Dirac peaks start to appear, as the periodic modulation gets stronger. At V_bg = ± 40 V (linecuts in Fig. 1d), 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 Fig. 1c shows a close-up of the boxed region with 30 V ≤ V_bg ≤ 45 V and − 3.5 V ≤ V_tg ≤ − 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⁴³, where two-terminal transport simulations were performed. The character of the bands can be verified in magnetotransport, as shown in the following subsections.

Low magnetic field

For a general band dispersion ϵ(k), the semiclassical equations of motion for an electron are given by⁴⁵

$$\dot{{{{\bf{r}}}}}=\frac{1}{\hslash }{\nabla }_{{{{\bf{k}}}}}\epsilon ({{{\bf{k}}}}),$$

(1)

$$\hslash \dot{{{{\bf{k}}}}}=-e({{{\bf{E}}}}+\dot{{{{\bf{r}}}}}\times {{{\bf{B}}}}),$$

(2)

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

$${{\Delta }}{{{\bf{r}}}}(t)=\frac{\hslash }{eB}{{\Delta }}{{{\bf{k}}}}(t)\times \hat{z},$$

(3)

meaning that the cyclotron orbit is obtained by rotating the orbit in the momentum space by 90° clockwise, as illustrated in Fig. 2b, c. Carriers encircle closed orbits of electron type or hole type in the counterclockwise or clockwise direction, respectively, as determined via the group velocity, v = ∇_kϵ(k)/ℏ, and the equation $\hslash \dot{{{{\bf{k}}}}}=-eB{{{\bf{v}}}}\times \hat{z}$. In pristine graphene at low energy, cyclotron orbits exhibit a circular shape (Fig. 2b), 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. Fig. 2c).

Fig. 2: The TMF schematics.

a Sketch of the geometry considered for magnetic focusing. Fermi contours in the reciprocal space for a circular (b) and noncircular (c) orbit, and the corresponding focused trajectories in real space.

In a typical device designed for transverse magnetic focusing measurement, an emitter and collector [in Fig. 2a marked as 1 and 2, respectively] are located on the same edge at a center to center distance ℓ = d + w, and other contacts act as absorbers. The current injected from an emitter flows along cyclotron orbits with a radius R_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 ℓ (Fig. 2b), or otherwise in absorber contacts. In the nonlocal resistance measurement (Fig. 2a), this results in maximum or minimum, respectively, of the resistance R_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³⁷. 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_bg = 40 V, such that a few higher-order Dirac peaks are present next to the main Dirac point as seen in Fig. 1c.

For transverse magnetic focusing, we choose the system geometry shown in Fig. 2a, 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 3a shows the R_f map as a function of B and V_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 Fig. 3b, where the main and secondary Dirac peaks are seen. The sign change of the focusing peaks in Fig. 3a 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_tg, occurring as an effect of the band reconstruction due to the superlattice potential. To understand the result in detail, in Fig. 3c–f we plot the miniband structures calculated at B = 0, as described in Ref. ⁴³, and the Fermi contours at E = 0, at selected values of V_tg marked by arrows in Fig. 3a. Additionally, in Fig. 3g–k we show zoomed regions of the R_f map marked with the colored rectangles in Fig. 3a.

Fig. 3: Magnetic focusing in the superlattice.

a Nonlocal resistance R_f as a function of magnetic field and top gate voltage at V_bg = 40 V. b R_xx calculated for the same system at B = 0. c–f Top row: the miniband structures calculated at V_tg marked by the arrows in (a). In (f) C₁, C₂, and V₁-V₄ label selected minibands. Bottom row: the corresponding Fermi contours at E = 0. g–k The zoomed map in (a) at the rectangles marked with the corresponding colors. The dashed lines mark the focusing condition ℓ/2n = ℏk/eB, n = 1, 2, … .

In Fig. 3c, at V_tg = − 2.75 V, the Fermi level is located at the C₁ subband, and the Fermi contour has a rounded shape (Figs. 2b, 3c). In the semiclassical description, electrons injected from lead 1 (Fig. 2a) with an initial velocity v = (0, v) and 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_c = 2ℏk/eB along the x direction (Fig. 2b). For ℓ = 2R_c, it leads to the first maximum of R_f. For smaller R_c, the beam is reflected at the edge, and can flow to the collector when ℓ = 4R_c, 6R_c, … , giving rise to higher order R_f peaks. In general, we can evaluate the field at which the jth maximum occurs as B_j = 2ℏjk/eℓ, j = 1, 2, … . We find k(V_tg) numerically and plot B_j with dashed lines in Fig. 3(g). The C₁ subband cone is within − 2.9V ≲ V_tg ≲ − 2.5 V. We find a very good agreement with the R_f signal for up to j ≈ 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.1V ≲ V_tg ≲ − 2.9 V, for the V₁ subband, R_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_tg ≈ − 2.5 for the electron subband, and V_tg ≈ − 3.1 for the hole subband), the focusing signal vanishes, and smaller fans reappear. Based on the miniband structures (Fig. 3d), we interpret them as due to focusing of the secondary Dirac cones fermions. For example, at V_tg = − 3.12 V (Fig. 3d) at E = 0 there are tiny Dirac cones around the M and X points of the Brillouin zone.

For V_tg ≲ − 3.2 V and V_tg ≳ − 2.4 V, the miniband structures around the Fermi level get more complex, with many overlapping subbands. Nevertheless, we find ranges of V_tg where a single isolated higher-order Dirac cone is present, giving clear focusing signal [see Fig. 3e for the V₂ subband, the corresponding R_f zoom in Fig. 3h, and the zoom in Fig. 3i for the C₂ 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₄ subband around the M point, see Fig. 3f at V_tg = − 3.85 V, and zoomed R_f in Fig. 3j. A similar feature is resolved in Fig. 3k 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 4a–c shows the line cuts of R_f at V_tg and B marked by the respective colors in Fig. 3g, h, j. In Fig. 4a for the C₁ subband, the first two focusing peaks are marked by ∘ and ⊲. The respective current density maps are presented in Fig. 4d, e, revealing typical current densities found for TMF calculations in pristine graphene^46,47. In Fig. 4b the line cut for the V₂ subband is shown. For the focusing peaks marked with ◇ and ▿, the current density maps are shown in Fig. 4f, g. The trajectories acquire a shape close to a rectangle, consistent with the Fermi contour (Fig. 3e). For the line cut in the hole-like V₄ subband (purple) (Fig. 4c), the current densities of the first two peaks are shown in Fig. 4h, i. The orbits acquire a rhombus shape, matching well the Fermi contour in the corner of the Brillouin zone (Figs. 2c, 3f). 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 k for which v_x ∝ ∂E/∂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 Fig. 4h). Note that in Fig. 4f–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 Fig. 4a–c originates from scattering and resonant states due to SL which are irregular and complex due to the wave-like nature of the carriers.

Fig. 4: Non-circular motion of Dirac fermions.

Line cuts at V_tg and B ranges as marked in Fig. 3g, h, j by the respective colors at (a) V_tg = − 2.75 V, (b) V_tg = − 3.3 V, and (c) V_tg = − 3.85 V. d–i are the current density maps in arbitrary units at B marked in (a–c) by the respective symbols. The red dashed lines show the expected semi-classical trajectory.

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_bg and V_tg, and the band structure is significantly modified merely by changing V_tg (Fig. 3c–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_f closer to straight lines, unlike Refs. ^39,40 that observed ones close to parabolic.

High magnetic field

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, ϕ = AB, is of the order of the flux quantum ϕ₀ = 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 Fig. 5a with the actually considered geometric dimensions shown. We compute all the 5 × 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^48,49, according to the lead labels shown in Fig. 5a. For all the following discussions, we convert our V_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.

Fig. 5: High-field transport.

a The 5-terminal system used for the calculation. b, c are ${G}_{{{{\rm{H}}}}}={R}_{xy}^{-1}$ as a function of average carrier density $\bar{n}$ and magnetic field B with V_bg fixed at 0 and 40 V, respectively, where horizontal lines at B = 2 T and B = 0.2 T mark the line cuts shown in (d) for G_H and (e) for R_xy. f Longitudinal resistance R_xx as a function of $\bar{n}$ and B. g Calculated Hofstadter spectrum corresponding to (c) and (f).

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

Without taking the inverse of R_xy, Fig. 5e shows the Hall resistance at B = 0.2T with V_bg = 0 [red, corresponding to the $\bar{n}$ range marked on Fig. 5b] and V_bg = 40 V [cyan, corresponding to the $\bar{n}$ range marked on Fig. 5c]. 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{{\rm{n}}}$ and B as Fig. 5b and c, Fig. 5f shows the longitudinal resistance R_xx with the back gate voltage fixed at V_bg = 40 V. Since the period of our gate-controlled superlattice is λ = 35 nm, and hence the square SL unit cell area A = λ² = 1225 nm², the condition ϕ/ϕ₀ = AB/(h/e) = 1 is reached at B ≈ 3.376 T ≡ B₁, which is in sharp contrast with the graphene/hBN moiré superlattice that requires B ≈ 24 T in order to reach ϕ/ϕ₀ = 1, because of its periodicity limited to λ ≲ 14 nm and hence the area $A=\sqrt{3}{\lambda }^{2}/2\lesssim 170\,{{{{\rm{nm}}}}}^{2}$. As is visible on Fig. 5f, the seemingly complicated R_xx map does exhibit certain features at B₁ and B₁/2, and vaguely at B₁/3 (marked by black triangles), corresponding to ϕ/ϕ₀ = 1, 1/2, 1/3, respectively, consistent with our calculation of the magnetic energy subbands shown in Fig. 5g, where the vertical axis of ϕ/ϕ₀ corresponds to exactly the same B range as Fig. 5f, as well as the average density range. Good consistency between the R_xx map of Fig. 5f and the Hofstadter butterfly shown in Fig. 5g can be seen. For methods adopted to calculate Fig. 5g, see “Methods”.

Discussion

In summary, we theoretically investigated transport in gated superlattices based on monolayer graphene. Our zero and low magnetic field transport calculations remain in a good agreement with the continuum model band structure calculated in presence of periodic modulation. We explored the potential of TMF for probing the intricate band structure of graphene with periodic modulation. It offers possibilities to study a plethora of phenomena in superlattices, and opens the door for studies of strongly correlated systems in twisted bilayer graphene⁵² or in bilayer graphene superlattices⁵³. By exploring the reconstructed band structure via magnetotransport calculations it is possible to engineer devices relying on directed electron flow due to the distortion of Fermi contour, as well as for other applications based on mini band electron optics. We also obtained high-magnetic-field response consistent with the Hofstadter spectrum calculated for a gated superlattice as a function of the gate voltage. Our modeling can be generalized to other superlattice geometries, and is promising for the investigations of future band structure engineered devices working in a broad range of magnetic fields.

Methods

Gated superlattice model

To model a realistic geometry, we consider SL induced by a patterned dielectric substrate^31,32 with a uniform global backgate underneath the dielectric layer, and the graphene sheet sandwiched between two hBN layers lying on top [Fig. 1a of the main text]. The hBN/graphene/hBN sandwich is covered by a global top gate. The voltage applied to the back gate V_bg controls the strength of the periodic modulation, while the top gate voltage V_tg is used to tune the carrier density across the SL. We follow the design of Ref. ³¹, with a square lattice etched in SiO₂ substrate, with a lattice period λ = 35 nm. We use a model function C_bg(x, y)

$${C}_{{{{\rm{d}}}}}(x,y)={C}_{{{{\rm{out}}}}}-({C}_{{{{\rm{out}}}}}-{C}_{{{{\rm{in}}}}})\tanh \left[\exp \left(-\frac{{x}^{2}+{y}^{2}}{{d}_{{{{\rm{smooth}}}}}^{2}}\right)\right],$$

(4)

$${C}_{{{{\rm{bg}}}}}(x,y)={C}_{{{{\rm{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_smooth = 7.5 nm is the smoothness of the modulation, C_out/e = 0.77 × 10¹¹ cm⁻²V⁻¹ and C_in/e = 0.22 × 10¹¹ cm⁻²V⁻¹, and a%b means the remainder after dividing a by b. The top gate capacitance is assumed to be C_tg/e = 1 × 10¹² cm⁻²V⁻¹. Using the parallel capacitor model, this roughly corresponds to the top hBN thickness d_t ≈ 16.6 nm, from C_tg/e = ε₀ε_hBN/ed_t, where ε₀ is the vacuum permittivity, ε_hBN = 3 is the dielectric constant of hBN, and − e is the electron charge. Figure 1b 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}_{{{{\rm{bg}}}}}(x,y)}{e}{V}_{{{{\rm{bg}}}}}+\frac{{C}_{{{{\rm{tg}}}}}}{e}{V}_{{{{\rm{tg}}}}}.$$

(6)

Assuming that the carrier energy in graphene is given by E = ± ℏv_Fk, where ℏ is the reduced Planck constant, v_F ≈ 10⁶ m/s is the Fermi velocity of graphene, and using $\hslash {v}_{{{{\rm{F}}}}}\approx 3\sqrt{3}/8\,{{{\rm{eV}}}}\,{{{\rm{nm}}}}$, the onsite potential energy can be obtained from

$$U=-{{{\rm{sgn}}}}(n)\hslash {v}_{{{{\rm{F}}}}}\sqrt{\pi | n| },$$

(7)

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

Transport calculation

For transport calculation, we use the tight-binding Hamiltonian

$$H=-\mathop{\sum}\limits_{\left\langle i,j\right\rangle }{t}_{ij}{c}_{i}^{{\dagger} }{c}_{j}+\mathop{\sum}\limits_{j}U({{{{\bf{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 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 B = (0, 0, B), the hopping integral is modified by ${t}_{ij}\to {t}_{ij}\exp (i\phi )$, where the Peierls phase $\phi =(-e/\hslash )\int\nolimits_{{{{{\bf{r}}}}}_{i}}^{{{{{\bf{r}}}}}_{j}}{{{\bf{A}}}}\cdot d{{{\bf{r}}}}$, with A being the vector potential that satisfies ∇ × A = B, and the integral going from the site at r_i to the site at r_j. For a feasible simulation of realistic devices, we use the scalable tight-binding model⁵⁴, where the hopping parameter becomes t_ij = t₀/s_f, and the lattice spacing a = a₀s_f, s_f is the scaling factor, and we use t₀ = − 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²/h)T_ji, where the transmission probability T_ji is evaluated as a sum over the propagating modes ${T}_{ji}=\mathop{\sum}\limits_{q}{T}_{ji}^{q}$, and

$${T}_{ji}^{q}=\mathop{\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 ${{{\boldsymbol{{{{\mathcal{G}}}}}}}}$^48,49 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}=\mathop{\sum }\nolimits_{j = 1}^{N}{{{{\mathcal{G}}}}}_{ij}{V}_{j}$. For an N-terminal system, the matrix elements are

$${{{{\mathcal{G}}}}}_{ij}=-{G}_{ij},\,\quad i\,\ne\, j,$$

(10)

$${{{{\mathcal{G}}}}}_{ii}=\mathop{\sum }\limits_{j=1,j\ne 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) × (N − 1) matrix $\bar{{{{\boldsymbol{{{{\mathcal{G}}}}}}}}}$ can be inverted to get ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}={\bar{{{{\boldsymbol{{{{\mathcal{G}}}}}}}}}}^{-1}$, where the ${{{\boldsymbol{{{{\mathcal{R}}}}}}}}$ matrix satisfies

$${V}_{i}=\mathop{\sum }\limits_{j=1,j\ne l}^{N}{{{{\mathcal{R}}}}}_{ij}{I}_{j}.$$

(12)

With the elements of matrix ${{{\boldsymbol{{{{\mathcal{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.

Transverse magnetic focusing

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³⁷, considering the same probe spacing (ℓ = 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_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 Fig. 6a, showing a map rather consistent with the experiment. Due to the isotropic low-energy dispersion of graphene, the resulting cyclotron trajectory is a simple circle of radius R_c = ℏk/eB, which is simplified from Eq. (3). By requiring the probe spacing to be equal to an integer times the cyclotron diameter, ℓ = j ⋅ 2R_c, j = 1, 2, ⋯ , together with $k=\sqrt{\pi | n| }$ for graphene, one can solve for carrier density corresponding to the jth peak of the TMF on the B-n map of Fig. 6a:

$${n}_{j}(B)=\frac{1}{\pi }{\left(\frac{eB\ell }{2\hslash j}\right)}^{2}\,.$$

(14)

The dashed lines on Fig. 6a show n₁, n₂, n₃, n₄, matching very well with the patterns of the simulated R_61,54, which requires totally 6 × 5 = 30 transmission functions for such a 6-terminal device, as explained above. Figure 6b, c shows two exemplary maps of transmission functions, which can look generally very different from the resulting four-point resistance.

Fig. 6: Revisiting the TMF experiment on graphene.

a Computed four-point resistance R_61,54 as a function of magnetic field B and carrier density n revisiting the TMF experiment on graphene³⁷, and two exemplary transmission functions required for R_61,54: (b) T₂₃ and (c) T₄₁. The inset in (a) shows the orientation of the considered Hall bar and the configuration of the leads for ground, injector, and voltage probes.

Choosing the gauge

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, A = ( − yB, 0, 0) or 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⁵⁵. Assuming that the proper gauge in the 1st lead is A₁(r), for another lead that is at an angle θ_n with respect to lead 1, the gauge can be chosen as

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

(15)

with

$${f}_{n}({{{\bf{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 ∇ × A = B. As an illustration of the transformation, consider A₁(r) = ( − yB, 0, 0), and θ_n = 90°. Then, f_n(r) = Bxy, ∇ f_n(r) = (By, Bx, 0), and A_n(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 ζ_n(r) which is only nonzero in the translationally invariant part of lead n

$${\zeta }_{n}({{{\bf{r}}}})=\left\{\begin{array}{ll}1,\quad &{{{\bf{r}}}}\,{{\mbox{in lead}}}\,n,\\ 0,\quad &{{{\bf{r}}}}\,{{\mbox{in lead}}}\,m\ne n,\\ \,{{\mbox{smooth interpolation}}}\,\quad &{{{\bf{r}}}}\,{{\mbox{elsewhere}}}\,.\end{array}\right.$$

(17)

Then, in (15) we substitute f_n(r) → ζ_n(r)f_n(r) for lead n. In general, for the entire system we define

$$f({{{\bf{r}}}})=\mathop{\sum }\limits_{n=2}^{N}{\zeta }_{n}({{{\bf{r}}}}){f}_{n}({{{\bf{r}}}}),$$

(18)

this completes our gauge transformation. Adopting the vector potential

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

(19)

we have B = ∇ × 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 ζ_n function.

As an example, for the system used for the TMF modeling (Fig. 2a), in the vertical leads we choose the same gauge A = (0, xB) with $\zeta ({{{\bf{r}}}})={(\exp (-(y-{y}_{1})/{d}_{{{{\rm{step}}}}})+1)}^{-1}+{(\exp (-({y}_{2}-y)/{d}_{{{{\rm{step}}}}})+1)}^{-1}$, y₁ = 1607 nm, y₂ = − 49 nm, and d_step = 2 nm. In the rest of the system, A = ( − yB, 0) is used. The resulting A_x and A_y profiles are shown in Fig. 7.

Fig. 7: Spatial profiles of the vector potential.

a A_x and (b) A_y vector potential components used for the 5-terminal TMF simulation.

Hofstadter butterfly calculation

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⁵⁴ with s_f ~ 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 (λ). Here, in order to ensure superlattice period equal to a multiple of 3a, s_f is not an integer, and the ribbon width is larger than 20λ 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_tg values and collect the Fermi states to construct the gate-dependent Hofstadter butterfly spectrum to compare with R_xx and G_H obtained from the transport calculation. Note that the spectrum in Fig. 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.