Mixing of moiré-surface and bulk states in graphite – Nature

Device fabrication

To make hBN/graphite/hBN heterostructures, graphite flakes were encapsulated by hBN through dry transfer as described elsewhere^38,39. In brief, graphite and hBN flakes were mechanically exfoliated onto oxidized silicon substrates. The target hBN flake was picked up by a polymer film made of polydimethylsiloxane and polymethylmethacrylate and then used to pick up a graphite flake of known thickness. The obtained stack is then released onto another hBN flake on an SiO₂/Si wafer, completing the heterostructure. To make aligned hBN/graphite structures, the straight edges of hBN and graphite flakes, which are usually along their crystallographic axes, were aligned in parallel.

Top-gate electrode and metal contacts to graphite (3 nm Cr/80 nm Au) were patterned using electron beam (e-beam) lithography and reactive ion etching, followed by an e-beam evaporation process. These devices were then shaped into Hall bar geometry using a thermally evaporated aluminium film as etch mask, which was later removed by 0.1 M NaOH solution. Alternatively, for Corbino devices, we utilized e-beam overexposure of polymethylmethacrylate resist to form a crosslinked bridge, which separates the inner contact, the top gate and the outer contact.

Graphite capacitor devices used to study surface states of non-aligned heterostructures were fabricated similarly, with the hBN flake intentionally misaligned to the underlying graphite film on a quartz substrate. A quartz substrate was chosen to minimize the parasitic capacitance, known to be a feature of SiO₂/Si substrates. Graphite flakes of around 50 nm thickness were used, guaranteeing the 3D graphite electronic spectrum. Relatively thick hBN flakes (>40 nm) were also chosen to eliminate the inhomogeneity of electrostatic potential introduced by a relatively rough metal electrode.

Transport and capacitance measurements

The longitudinal and Hall voltages of Hall bar devices were recorded with lock-in amplifiers (SR830 or MFLI) on applying a small low-frequency a.c. current of 10 nA (except where a higher current is specified). For Corbino devices, a small ac bias (40–100 μV) was applied to the inner contact, and the current was recorded from the outer one using SR830 in current input mode (lock-in amplifier input resistance of 1 kΩ and any in-line filters were subtracted from the measured resistance to account for any voltage drop across these components). The conductivity of Corbino devices was calculated using σ_xx = 1/(2π)ln(r_o/r_i)G, where G is the measured conductance and r_o is the outer radius and r_i is the inner radius of the graphite channel. Magnetic fields up to 18 T were generated by superconducting magnets, while data above 18 T were obtained in a 20 MW resistive magnet at the LNCMI-Grenoble.

To confirm the alignment or misalignment of the top and bottom graphite/hBN interfaces and to extract the moiré wavelength λ, we used the following two methods: (1) measuring the Landau fan diagrams of surface states at each interface by sweeping either the top or bottom gate voltages, respectively, and (2) measuring high-temperature Brown–Zak oscillations. Because the two surface states are electronically decoupled, they can feel only the potential in the vicinity of the corresponding interface because of electrostatic screening. For doubly aligned device D3, both surface states show Brown–Zak oscillations with conductivity peaks at nearly the same B fields, indicating the same moiré period for the top and bottom interfaces. We fitted multiple oscillations corresponding to integer flux fractions ϕ/ϕ₀ from 1/2 to 1/8 to σ_xx and derivatives \(\frac{{\rm{d}}{\sigma }_{xx}}{{\rm{d}}B}\) and \(\frac{{{\rm{d}}}^{2}{\sigma }_{xx}}{{\rm{d}}{B}^{2}}\) (Extended Data Fig. 3c,d), which yields a value of B₀ for each sequence of oscillations, where B₀ is the magnetic field at which ϕ₀ = B₀A₀ with magnetic flux quantum ϕ₀ and superlattice unit cell area A₀. Using this value of B₀, the moiré wavelength λ is calculated as λ = \(\sqrt{2{A}_{0}/\sqrt{3}}\). Furthermore, B₀ can also be extracted from fractal QHE states fitting at low temperatures. The difference in λ between the two interfaces calculated using these two methods was less than or equal to 0.1 nm. Given the similarity in the measured λ from Brown–Zak oscillations and the appearance of only one set of fractal states in the dual-gate maps in Fig. 3 (where distinct moiré periods would be expected to generate multiple sets of fractal states), we confirm the alignment of both interfaces with a precision of ±0.1 nm.

Differential capacitance C was measured as a function of bias voltage V_b between the metal gate and graphite using an on-chip cryogenic bridge⁴⁰, which reaches a sensitivity of about 10 aF at 1-mV excitation. Excitations of 102.53-kHz frequency and opposite phases were applied to the sample and a reference capacitor. Output signals from these two capacitors were mixed at the gate of a high-electron-mobility transistor, which served as an amplifier. The excitation voltage of the reference capacitor was modulated so that the output signal from the high-electron-mobility transistor becomes zero, and the capacitance of the sample is obtained from the ratio of excitation voltages at the balance point. A typical excitation voltage applied to the samples ranged from 1 to 10 mV, depending on the thickness of the hBN dielectric layer.

Surface states in non-aligned graphite films in a zero-B field

To compute the surface states, we adapted an effective-mass model of a finite-thickness graphite film using the Slonczewski–Weiss–McClure (SWMC) parameterization^26,31,32 combined with the self-consistent potential profile of graphite sandwiched between two gates with carrier densities n_t and n_b. In the Hartree approximation, the potentials on the layers, U_j > 1, are related to layer electronic densities n_l as

where ε = 2.6 accounts for the vertical polarizability of graphene⁴¹, c = 3.35 Å is the interlayer separation and 2N is the number of graphene layers. We temporarily fix the value of U₁, which has the role of a surface chemical potential and then self-consistently calculate Hartree potentials and densities on all the layers. The electronic density in layer \(l\) of graphite, calculated in the Hartree approximation, is

where f is a Fermi–Dirac distribution, and n enumerates the eigenfunctions for a given in-plane momentum k. The constant n₀ is chosen to match \({\sum }_{l=1}^{2N}{n}_{l}=0\) to provide electrical neutrality. After finding the densities on all the layers, we relate U₁ to n_t using \({n}_{{\rm{t}}}=-{n}_{{\rm{b}}}-\,{\sum }_{n=1}^{2N}{n}_{l}\).

To examine the thermodynamic density of states (DOS) at the graphite–hBN interface, we use capacitance spectroscopy, which is a tool that has been applied to 2D systems^40,42. However, its application to study surface states of metals or semimetals is rare. The measured capacitance (C) can be considered as geometric parallel-plate capacitance C_G = εε₀A/d and quantum capacitance C_Q in series, 1/C = 1/C_G + 1/C_Q, where A is the device area, ε₀ is the vacuum permittivity and d and ε are the thickness and relative permittivity of the hBN dielectric layer⁴³. The quantum capacitance reflects the DOS = dn/dU₁ on the surface of graphite: C_Q = Ae²dn/dU₁, where n is the carrier density and e is the electron charge. In our measurements, C follows a V-shaped dependence on n, where \(n({V}_{{\rm{b}}})=\frac{1}{Ae}{\int }_{0}^{{V}_{{\rm{b}}}}C\left(V\right){\rm{d}}V\), with a notable fine structure (Extended Data Fig. 1a).

The capacitance (per unit area) can be calculated self-consistently from equations (3) and (4) as

By comparing the calculated capacitance with experimental data, we obtained a set of SWMC parameters (γ₀ = 3.16 eV, γ₁ = 0.39 eV, γ₂ = −17 meV, γ₃ = −0.315 eV, γ₄ = 44 meV, γ₅ = 38 meV and Δ_AB = 50 meV). Results of this procedure are shown in Extended Data Fig. 1b, showing excellent agreement between theory and experiment. At low doping (|U₁| < |γ₂|, that is, |n| < 6 × 10¹¹ cm⁻²), there are no surface states, and quantum capacitance plotted in Extended Data Fig. 1b is determined by electron and hole screening. Because holes have a slightly larger DOS (shallower dispersion) than electrons, we see larger C_Q at hole doping. When doping reaches U₁ ≈ ±γ₂, type 1 and type 2 surface states appear and contribute to the quantum capacitance. The radius of the surface Fermi line for type 1 states grows with |n|, leading to an increase in the density of surface states and growth of C_Q. Examples of the graphite film dispersion spectra with self-consistently determined layer potentials are shown in Extended Data Fig. 1c–g for n ranging from −6 ×  10¹²  cm⁻² to 6 × 10¹² cm⁻², where the red colour coding represents a high probability for wavefunctions at the first graphene bilayer and the green colour coding represents a high probability for wavefunctions at the second graphene bilayer.

To provide a qualitative understanding of the surface states in graphite, we analytically solve the spectrum of graphite³², considering the boundary conditions (Ψ = 0 for surface carbon atoms) and plot the eigenstates for homogeneous bulk graphite (Extended Data Fig. 1h,i), which consist of a propagating mode (black curves, real k_z) and an evanescent mode (orange curves, complex k_z). There are no complex k_z solutions at zero doping, as only real k_z solutions satisfying zero boundary conditions can be normalized. Electrostatic doping of graphite surface creates an inhomogeneous z direction potential near the surface, which does not preserve k_z momentum, allowing real k_z solutions near the surface, which then turn into evanescent modes decaying into the bulk. This provides a heuristic picture of the origin of non-trivial surface-state solutions (Extended Data Fig. 1h,i).

There are three types of propagating mode: majority electron and hole bands with bandwidth 2γ₂ and a minority carrier band near ck_z = π/2. These propagating bands cross the bulk Fermi level at a small in-plane momentum k_x,k_y (Extended Data Fig. 1h, 1D metal regime in the z direction) but are spread away from the Fermi level at large k_x,k_y (Extended Data Fig. 1i, 1D semiconductor regime in the z direction). When a potential near the surface is introduced by doping, the propagating modes in the 1D semiconductor region start to cross the Fermi level. With potential abating away from the surface, these modes evolve into evanescent modes in the gap (Extended Data Fig. 1i, green arrows for electron doping and blue arrows for hole doping). The dispersion of these evanescent modes, which we denote as type 1, crosses the Fermi level, forming a surface Fermi line with a radius larger than the Fermi surface of propagating carriers (Extended Data Fig. 1c–g, yellow contours). These states are similar to surface states in doped semiconductors, with the difference that they exist for only in-plane momenta larger than the projection of the bulk Fermi surface of graphite (no surface states observed for zero doping in Extended Data Fig. 1g). In the 1D metal regime, another type of evanescent mode, denoted as type 2, appears for |E| > |γ₂| and never crosses the Fermi level (Extended Data Fig. 1h).

Surface states in graphite films in the presence of a moiré superlattice

The spectrum for graphite aligned with hBN was calculated by treating the periodic moiré potential as a perturbation applied to only the top graphene layer. We followed the standard procedure¹⁹, using the mirror-symmetric superlattice coupling Hamiltonian \(\delta H={\sum }_{m=0}^{5}{{\rm{e}}}^{{{\rm{i}}{\bf{g}}}_{m}\cdot {\bf{r}}}\,\left[{U}_{0}+{(-1)}^{m}({\rm{i}}\,{U}_{3}{\sigma }_{3}+{U}_{1}\frac{{{\bf{g}}}_{m}\times \hat{z}}{|g|}\sigma ){\tau }_{3}\right]\) applied to the two top-layer components of the graphite film wavefunction, where Pauli matrices σ operate on top-layer sublattices and τ operates on valleys, \({{\bf{g}}}_{m}={{\mathcal{R}}}_{\pi (m-1)/3}\{0,\,4\pi \delta /(3a)\}\) are six reciprocal lattice vectors of superlattice (where \({\mathcal{R}}\) is a rotation matrix), \(\delta =0.018\) is a lattice mismatch, a = 1.42 Å is carbon–carbon distance, and we use the parameters U₀ = 8.5 meV, U₁ = −17 meV and U₃ = −14.7 meV (refs. ^44,45). The results do not significantly depend on the values of superlattice couplings, and it was sufficient to restrict the momentum space to the first star of the superlattice reciprocal lattice vectors to achieve convergence.

At low fields (B < 1 T), the onset of 2.5D QHE is strongly altered by the kaleidoscopic band structure of the surface states (Fig. 1e). We compare the low field transport for aligned (D1) and non-aligned (D4) devices of similar graphite thickness (approximately 8  nm) (Extended Data Fig. 2a–d). In a non-aligned graphite device, we observe that a Landau fan develops for finite densities |n_b| > 10¹² cm⁻², and all QHE states can be traced back to n_b ≈ 0 as B approaches 0. By contrast, for aligned graphite similar QHE features are also overlaid by oscillations emanating from LTs at |n| ≈ 2.0 and 3.7 × 10¹² cm⁻² resulting in the diamond-like features in σ_xx occurring at flux fractions ϕ/ϕ₀ = p/q. Comparison of low field conductivity as a function of tuning aligned and non-aligned interfaces in the same device also shows pronounced differences, as shown in Extended Data Fig. 2e,f, where the most visible features occur only at |n_b| > 2 × 10¹² cm⁻², independent of n_t doping.

To highlight Brown–Zak oscillations across a large range of magnetic fields, we also calculated Δσ_xx by subtracting a smooth background from the σ_xx data. In comparison to graphene–hBN systems in which the background conductivity can be fitted with polynomials²¹, we find that even-higher-order (>10) polynomials are insufficient as many oscillatory artefacts are present. Instead, we use a two-carrier Drude model of σ_xx(B) and σ_xy(B) and fit both simultaneously to yield carrier densities and mobilities n_h = 2.2 × 10¹² cm⁻², µ_h = 24,000 cm² V⁻¹ s⁻¹, n_e = 2.8 × 10¹² cm⁻² and µ_e = −19,000 cm² V⁻¹ s⁻¹ for zero gate bias at T = 60 K. This two-carrier model fit, \({\sigma }_{xx}^{{\rm{fit}}}(B)\), is then used to calculate \({\Delta \sigma }_{xx}\left({n}_{{\rm{b}}},B\right)={\sigma }_{xx}\left({n}_{{\rm{b}}},B\right)-{\sigma }_{xx}^{{\rm{fit}}}\left(B\right)\). Oscillations in Δσ_xx occurring at \({B}_{1/q}\) visible for q ≤ 11 (Fig. 2b and Extended Data Fig. 3a) were cross-examined against raw σ_xx data to confirm they were not introduced by the subtraction process.

Bulk states in graphite films in the presence of a surface moiré superlattice

To model the transport gaps in our aligned devices at high B and low T, we treat moiré superlattice potential as a weak perturbation; each 2.5D QHE Landau level (Fig. 4b) is split into q subbands, at a given ϕ/ϕ₀ = p/q. Levels in Fig. 4b were calculated from the tight-binding description of Landau bands in graphite from ref. ²⁶ using the same set of SWMC parameters as stated above, with an adjustment to the splitting between Landau bands 0 and 1 attributed to the effects of self-energy in high-B fields^46,47. Hofstadter’s butterfly for a honeycomb lattice (Fig. 4a) was calculated from the finite-difference equation in ref. ³⁶, in which the energy scale has been normalized, and is given in arbitrary units, ε = ±1. A limit of q ≤ 100 was used for the computation to give a balance between plot density and speed. However, this results in the apparent absence of states near ϕ/ϕ₀ = 1, 1/2 and 1/3 (Fig. 4a,c), and it should be noted that this is a feature of the computation, not gaps in the spectra.

Analysis of the thermal activation of gaps for device D3 at B = 13.5 T (Fig. 3e) indicates the largest fractal gaps are ΔE_fractal ≈ 0.1 meV. We assign ΔE_fractal to the largest gap in Hofstadter’s butterfly at the flux value ϕ/ϕ₀ = 0.57 (corresponding to B = 13.5 T; Fig. 4a, dashed line), which spans 0.32 < ε < 0.56. This yields a scaling factor S = 0.42 meV. The full spectrum is then calculated using equation (2) and shown in Fig. 4c, in which we use the periodicity of Hofstadter’s butterfly (such that ε(ϕ/ϕ₀ + ρ) = ε(ϕ/ϕ₀) for any integer ρ) to plot states at ϕ/ϕ₀ > 1.

For comparison, the fractal energy spectrum was also computed for device D2, which has a different alignment to hBN and layer parity to that of device D3 (device D2 is 21 layers in thickness and aligned to only one encapsulating hBN). Odd-layer parity lifts the ±KH valley degeneracy in 2.5D QHE in graphite²⁶ and, therefore, the gap size is significantly reduced (Extended Data Fig. 6a–c) and the maximal gap size is about 0.9 meV (compared with 1.8 meV in device D3; Fig. 3e). In Extended Data Fig. 6d, we focus on the evolution of gap size at filling ν = 0 between two level crossings at B = 10 T and B = 16 T, with a maximal observed gap of about 0.48 meV. Extended Data Fig. 6f shows the Landau levels without moiré perturbation for the 21-layer graphite. Both the extent of the ν = 0 gap (8.5 T < B < 17 T) and its maximal size (1.3 meV) in the model are notably larger than those observed in the experiment. On application of Hofstadter’s butterfly to each Landau level (using the same S = 0.42 meV as in Fig. 3e), each level has effectively broadened and thus the ν = 0 gap in the model is reduced to about 0.6 meV (Extended Data Fig. 6g), in closer agreement with our experiments. However, the broadening of energy levels from Hofstadter’s butterfly leads to many overlapping states and hence gap closures, which were not observed in our experiment. This is probably because of inadequate treatment of moiré perturbation to states hosted on even and odd layers. As the moiré reconstruction is limited to a very short penetration depth, the perturbation will be much larger on the outermost layer (odd) than on subsequent layers. We plot a revised model with S = 0.42 meV for odd layers and S = 0.12 meV for even layers (Extended Data Fig. 6h), yielding fewer gap closures, whereas the same ν = 0 gap remains, which is in overall better agreement with our experiment (Extended Data Fig. 6e).

Surface states in non-aligned graphite films in finite B fields

In the B field, surface states manifest in the capacitance spectra as pronounced magnetocapacitance oscillations (Extended Data Fig. 7a). For bulk Landau bands that cross the Fermi level, the associated surface states would coexist and mix with them. However, bulk Landau bands away from the Fermi level can become occupied at the surface when electrostatically doped, giving rise to surface Landau levels. On filling these surface Landau levels, regions of high compressibility appear as peaks in the capacitance spectra. Note that the width of these high-compressibility regions does not correspond to integer degeneracy (>4), because some fraction of gate-voltage-induced charge is sunk into the bulk to support the self-consistent screening potential near the surface (Extended Data Fig. 8b).

Having determined the geometric capacitance from the fitting of zero field data, we can convert C(n) into DOS(U₁) using U₁ = eV_b − e²n/C_G (ref. ⁴⁰). As shown in Extended Data Fig. 7b, peaks in the DOS correspond to metallic-surface Landau levels, which are separated by relatively low DOS regions (cyclotron gaps of the surface states). In contrast to true 2D systems, the DOS in these cyclotron gaps is non-zero, because charges can be injected into the bulk graphite. At 12 T, three minima are further developed on top of most peaks, indicating that the fourfold degeneracy (spin and valley) of the surface Landau levels is lifted.

Experimental results are better visualized and more informative when presented as a C(n, B) map (Extended Data Fig. 7c). The branches of surface states spawn out from the neutrality points at B = 7.5 T, 3 T, 2 T and so on. These B fields correspond to the critical fields above which the bulk Landau bands no longer cross the Fermi level and appear as only surface Landau levels. For instance, according to our SWMC model, at 7.5 T, the bulk Landau band 2⁺ is just above the Fermi level (Extended Data Fig. 7d). Thus, a branch of surface states spawned out around this field is labelled as S²⁺. The same happens with the electron bulk Landau band 3⁺ at 3 T and hole bulk Landau band 2⁻ at 2 T (Extended Data Fig. 7c).

We observed oscillations down to B ≈ 0.1 T (Extended Data Fig. 8a), which sets a lower bound of approximately 100,000  cm²  V⁻¹  s⁻¹ for surface-charge carrier mobility. The high electronic quality of surface states also enables fractional features in the Landau quantization of charge carriers. A graphite capacitor device was fabricated to investigate fractional QHE features, with a thicker hBN dielectric to reduce the inhomogeneity of electrostatic potential from the metal electrode. At a high magnetic field, B = 20 T, we observe the formation of two minima on top of singly degenerate surface states of S²⁺ (Extended Data Fig. 9a,b). The Δν between the fractional gap is around 0.27, which is lower than the expected Δν = 1/3 for fractional QHE. To further investigate these fractional QHE states, we used thin (6 nm) graphite (device D9) and studied transport under an applied displacement field, D = (n_t − n_b)e/2ε₀. At D = 0.24 V nm⁻¹, B–n regions can be found in which the energy level of surface states locates in the bulk gap (Extended Data Fig. 9c,d). In these regions, the surface states are isolated from the bulk completely, and vanishing σ_xx and quantized σ_xy indicate the development of fractional QHE with a 1/3 degeneracy. The difference between the capacitance and transport measurements can be reconciled by considering the negative compressibility of the fractional states: the chemical potential of the surface states reduces with the injection of n, acquiring additional charges from the bulk^48,49,50.

Conventional interpretation of Brown–Zak oscillations

The classical dynamics of the electron is set by \(\dot{{\bf{p}}}=eB\hat{{\bf{z}}}{\boldsymbol{\times }}\dot{{\bf{r}}}\) and \(\dot{{\bf{r}}}\equiv {\bf{v}}={{\nabla }}_{{\bf{p}}}\varepsilon ({\bf{p}})\), which implies that the real-space trajectories can be obtained from constant energy contours in momentum space by a 90° rotation and rescaling by 1/eB. Near Van Hove singularities, caused by the saddle points in the dispersion of electrons, the change of the sign of the band mass occurs, which is known as the LT. At the LT, closed cyclotron orbits of electrons transform into open trajectories, forming a network, which, because of the C₃ symmetry of graphite film, looks like a Kagome pattern. This leads to delocalized electron orbits resulting in high conductivity even at strong magnetic fields, even though the electron ballistic motion along such a network has a stochastic element: when reaching the saddle points in dispersion, electron paths can switch between electron- and hole-like segments (Extended Data Fig. 4a,b). This process, known as the magnetic breakdown of cyclotron motion, can be captured^51,52,53,54 by transmission amplitudes, \(\mathop{S}\limits^{\leftharpoonup }\) and \(\mathop{S}\limits^{\rightharpoonup }\),

These magnitudes of \(| \mathop{S}\limits^{\rightharpoonup }| \) and \(| \mathop{S}\limits^{\leftharpoonup }| \) are comparable to each other in the magnetic breakdown⁵⁵ interval of energies, proportional to \(eB{\mathfrak{r}}/(2\pi \hbar )\), which is determined by the strength of B field and Gaussian curvature of the dispersion saddle point, \({\mathfrak{r}}=\hbar \sqrt{| \det \frac{{\partial }^{2}\varepsilon ({\rm{p}})}{\partial {p}_{i}\partial {p}_{j}}| }\) and sets the energy window around E_LT, where the LT network of trajectories is relevant for electron transport.

For any pair of points in the network, there are several distinct equal-length paths connecting them. These paths consist of an equivalent set of segments passed in a different order (for example, green and brown paths in Extended Data Fig. 4b), which—because of the periodicity of the Kagome network—ensures independence of the interference phase between partial waves following those paths, on the exact energy of electrons (similar to the physics of weak localization). As a result, broadening of the Fermi step does not lead to self-averaging of constructive and destructive interference contributions generated by electrons at various energies (as happens with the interference-induced mesoscopic fluctuations). The length of each segment of the trajectory scales as 1/B. So, for low B, only the shortest possible trajectories retain ballistic propagation (see Extended Data Fig. 4b for examples of such pairs of trajectories). The area, \(A=\frac{{A}_{{\rm{BZ}}}}{{(eB)}^{2}}=\frac{{(2{\rm{\pi }})}^{2}}{{{A}_{0}(eB)}^{2}}\), between the pairs of such trajectories is related—by rescaling with B field—to the actual Brillouin zone area, \({A}_{{\rm{BZ}}}={(2{\rm{\pi }})}^{2}/{A}_{0}\), of the superlattice, where \({A}_{0}\) is the unit supercell area. Multiplied by the magnetic field, this determines the encircled magnetic flux \(\phi ={A}_{0}B\) and the Aharonov–Bohm phase, \(\varphi =\hbar eAB=\hbar \frac{{(2{\rm{\pi }})}^{2}}{{A}_{0}eB}=2{\rm{\pi }}\frac{{\phi }_{0}}{\phi }\). The interference between partial waves that undergo ballistic propagation along the Kagome trajectories and stochastic switching at the Kagome network sites produce conductivity oscillations,

which are 1/B periodic. At low magnetic fields, the length of the paths, \({\mathcal{L}}(B)\propto {B}^{-1}\), would be longer than the shortest of the mean free path and coherence lengths, \({\mathcal{l}}\), which is captured by the exponential factor in equation (6). Here ρ is the thermodynamic density of states and the width, \(\delta n\), of the doping interval around LT in which the oscillations are expected to be visible is determined by both T and B. Note that the described oscillations are related to the ability of the electron to propagate across the superlattice rather than its density of states. Thus, they appear in the conductivity measurements but would be absent in quantum capacitance measurements. We also note that for a graphite surface aligned with hBN, LTs start to appear in surface and mixed bulk-surface bands after surface doping reaches 2 × 10¹² cm⁻². With increasing doping, there is a cascade of LTs (Extended Data Fig. 4a,b). As a result, the interval of densities, in which the above-described 1/B-periodic oscillations are visible, is broadened.

The mean free path, \({\mathcal{l}}\), appearing in the denominator of the exponent of our expression for the amplitude of oscillations in equation (6) is usually assumed to depend on only the temperature, and equation (6) produces an exponential decay, \({\rm{\exp }}\left(-\frac{2{\mathcal{L}}}{{\mathcal{l}}}\right)\), of oscillations at low magnetic fields. However, in the case of graphite, there are a lot of Fermi contours of bulk bands nearer to the K point and the increasing magnetic field could result in magnetic breakdown scattering from the surface to the bulk bands. This additional scattering decreases the lifetime of electrons on surface-state Fermi contours, effectively decreasing the mean free path, \({\mathcal{l}}\), with growing B. This mechanism may lead to non-monotonic dependence of the amplitude of oscillations on 1/B (Extended Data Fig. 4e), reflecting the complexity of our 3D twistronics system.

Raman spectroscopy of aligned graphite films

To characterize the effect of surface superlattice potential on hBN-encapsulated graphite, we performed Raman spectroscopy measurements. A graphite flake with an extended monolayer graphene (MLG) region was selected to benchmark the alignment of the entire graphite film (Extended Data Fig. 10a,b). Raman spectra of MLG/hBN superlattices have been well studied⁵⁶, and the alignment can be traced by the width of the 2D peak of MLG. The 2D peak of MLG broadens with better alignment because of the increased strain inhomogeneity caused by the moiré periodic potential of the hBN substrate. Similar broadening of the 2D peak was also observed in the bilayer graphene–hBN superlattice system⁵⁷, indicating that the superlattice potential of the hBN substrate can propagate through graphene bilayers, and is therefore detectable by Raman spectroscopy. However, how far this superlattice potential can penetrate the bulk of graphite remains unclear.

To clarify this, we fabricated two hBN/graphite/hBN heterostructures at the same time, by transferring graphite onto two adjacent but intentionally misoriented hBN flakes. The graphite flake is controlled to be aligned with one of the hBNs, and as a consequence is misaligned with the other (Extended Data Fig. 10c). The flake alignment is characterized by the full width at half maximum (FWHM) of the MLG 2D peak (Extended Data Fig. 10e). Each spectrum was averaged over ten spectra acquired at different positions and normalized by the intensity of the E_2g hBN peak at 1,363 cm⁻¹. The FWHM is 21 cm⁻¹ and 35 cm⁻¹ for non-aligned and aligned regions of the MLG, respectively, which agrees well with the results in ref. ⁵⁶. Broadening of the 2D peak is expected if the superlattice potential at the interface can propagate through the bulk graphite crystal. We found no appreciable difference on the Raman map of 2D FWHM between aligned and non-aligned graphite regions (Extended Data Fig. 10d,e). This implies that the surface superlattice potential of the hBN substrate does not penetrate through graphite, at least for films of thickness at least 2.6 nm.