Mixing of moiré-surface and bulk states in graphite

Van der Waals assembly enables the design of electronic states in two-dimensional (2D) materials, often by superimposing a long-wavelength periodic potential on a crystal lattice using moiré superlattices1–9. This twistronics approach has resulted in numerous previously undescribed physics, including strong correlations and superconductivity in twisted bilayer graphene10–12, resonant excitons, charge ordering and Wigner crystallization in transition-metal chalcogenide moiré structures13–18 and Hofstadter’s butterfly spectra and Brown–Zak quantum oscillations in graphene superlattices19–22. Moreover, twistronics has been used to modify near-surface states at the interface between van der Waals crystals23,24. Here we show that electronic states in three-dimensional (3D) crystals such as graphite can be tuned by a superlattice potential occurring at the interface with another crystal—namely, crystallographically aligned hexagonal boron nitride. This alignment results in several Lifshitz transitions and Brown–Zak oscillations arising from near-surface states, whereas, in high magnetic fields, fractal states of Hofstadter’s butterfly draw deep into the bulk of graphite. Our work shows a way in which 3D spectra can be controlled using the approach of 2D twistronics.

At the surface of a crystal, the periodic lattice is interrupted.Surface states arise, as itinerant Bloch wavefunctions with real-valued crystal momenta, k, become evanescent states localised at the surface with complex kz values [25].For example, surface charge accumulation in gapped semiconductors leads to distinct two-dimensional subbands easily tuned by electrostatic gating.While in metals, the high charge carrier density makes it difficult to observe and control surface states, as the bulk shunts the surface conductivity.Lying in between these two extremes are semimetals like bismuth and graphite, whose tunable surface states are interesting but not well explored yet.In this work, using van der Waals engineering and twistronics, we explore the gallery of moiré structuring of highly-tunable electronic states at the semimetal graphite surface, by aligning two bulk crystals, hexagonal graphite (Grt) and hexagonal boron nitride (hBN), to each other.
To this end, we prepared hBN/Grt/hBN heterostructures by aligning thin graphite films (typically, 5-10 nm thick) on top of hBN substrate and encapsulating the stack with another hBN crystal.Unless otherwise stated, this encapsulating hBN is intentionally misaligned (see Methods for details).As the lattice constants of hBN and graphite are similar, in a heterostack, they form a moiré superlattice with the periodicity controlled by the lattice mismatch, δ = 1.8%, and a misalignment angle, ϑ (Fig. 1A).In addition to providing moiré superlattice, the hBN encapsulation also preserves high electronic quality of graphite films [26][27][28].Figure 1A-C shows schematics and micrographs of hBN/Grt/hBN heterostructures, fabricated into Hall bar and Corbino geometry devices (in total we have studied 10 graphite heterostructure devices).In these devices, top and bottom electrostatic gates were used to independently control carrier densities nt and nb, at the top and bottom interfaces of hBN/Grt/hBN heterostructure.Hexagonal graphite (Bernal stacking) is a compensated semimetal with Fermi surface occupying only a small fraction of Brillouin zone.The size of the Fermi surface is determined mostly by a through-layer hopping parameter, γ2 ≈ -20 meV [29].Due to its semimetallic nature, graphite does not host surface states (evanescent modes) in the absence of dangling bonds or electric fields on its surfaces.However, when electric field above a critical value is applied perpendicular to graphite basal plane, tunable surface states start to appear, see refs.[26,30] and Supplementary Note 1.
Moiré superlattice drastically modifies graphite surface states, resulting in entirely different transport behaviour between aligned and non-aligned devices, Fig. 1D.Non-aligned graphite interface shows nearly featureless carrier density dependence of longitudinal σxx(n) and transversal σxy(n) conductivities in small magnetic fields.However, when gating aligned graphite interface, σxy(n) displays multiple zero crossings accompanied by the corresponding peaks in σxx(n).We attribute these crossings to the recurrence of electrostatically induced surface states dominated by electron-or hole-like charge carriers.To quantify this, we calculated Fermi surface projections using an effective-mass model with Slonczewski-Weiss-McClure (SWMC) parameterisation of graphite [31,32] subjected to moiré superlattice potential, in combination with self-consistent Hartree analysis (see Methods).Our calculations in superlattice Brillouin zone (SBZ) of hBN/Grt/hBN heterostructure show kaleidoscopic surface states with numerous topological Lifshitz transitions (LTs) across a range of carrier densities, shown in Fig. 1E.Pairs of plots (Fig. 1E, a-d) with drastic changes in the Fermi surface topology demonstrate four LTs, corresponding to the four carrier density ranges highlighted in Fig. 1D.LTs observed at |n|≈ 2 and at |n|≈ 3.7×10 12 cm −2 belong to two different branches of surface states, the former ones mostly residing on the first graphene bilayer in Grt, while the latter are more peaked at the second bilayer (Fig. S1C).As magnetic field increases the surface states in the vicinity of LTs give rise to separate branches of Landau levels; the evolution of σxx(n) and σxy(n) in low magnetic fields is further illustrated in Supplementary Note 2.
At higher B fields, the difference between aligned and non-aligned interfaces of hBN/Grt/hBN system becomes even more prominent.Figure 2A plots σxx as a function B for surface states on aligned and nonaligned interfaces.The curves were measured at T = 60 K to supress the effects of Landau quantisation.When the aligned surface is doped away from electron-hole compensation point, σxx exhibits oscillatory behaviour periodic with 1/B.Peaks in σxx can be traced to values of  1/ = 1   0  0 corresponding to integer number q of superlattice unit cells A0 = √3/2 λ 2 , commensurate with a magnetic flux quantum ϕ0 = h/e, here λ is the wavelength of moiré superlattice, h is the Planck's constant, and e is the elementary charge.The commensurability between ϕ0 and magnetic flux through a moiré unit cell ϕ = BA0 can be interpreted as the manifestation of Brown-Zak quantum oscillations at the superlattice interface, which were recently reported for aligned graphene/hBN heterostructures [21,22].Similarly, the formation of magnetic Bloch states leads to higher conductivity due to straight rather than cyclotron trajectories of surface charge carriers [21,22,33,34], as evidenced by the conductivity peaks at B1/q, Fig. 2A. Figure 2B shows these nbindependent conductivity peaks located at all distinguishable 1/q commensurate fields.Note that not only unit fractions but also second-order fractal states (e.g., B2/5 in Fig. S3A) are evident in Brown-Zak oscillations.Since Brown-Zak oscillations stem from translational invariance of magnetic Bloch states at rational fractions of magnetic flux ϕ/ϕ0 = p/q, they are expected to be insensitive to temperature as long as electrons retain phase coherence within magnetic supercell qA0. Figure 2C shows that at intermediate temperatures (T = 20K), states with ϕ0/ϕ = 24 are clearly visible (and ϕ0/ϕ > 35 states are distinguishable in Fig. S3), providing a lower bound on phase coherence length > 85 nm.Brown-Zak oscillations can also be interpreted as the Aharonov-Bohm interference of electron trajectories in a periodic 2D network formed by classical trajectories of electron drifting around the Fermi contours in a magnetic field that are joined by magnetic breakdown tunnelling in the vicinity of Van Hove singularities [35].This interpretation allows a convenient transition to the regime of low B-fields where we see multiple Lifshitz transitions of the Fermi surface topology, as shown in Fig. 1E, and explains the disappearance of Brown-Zak oscillations for |nb|< 2 × 10 12 cm -2 .
In comparison, gating the non-aligned interface in our hBN/Grt/hBN heterostructures does not result in Lifshitz transitions or Brown-Zak oscillations (cf.Fig. 1D, Fig. 2A, and Supplementary Note 3).This is not surprising, as our previous study has shown that the surface states on the two surfaces of a graphite film are well screened from each other, with the screening depth of the order of 2-3 layers [26].Raman measurements also do not show any qualitative difference in strain distribution or effects of alignment for films thicker than 7-8 graphene layers at the aligned and non-aligned graphite surfaces (see Supplementary Note 4).Our results are further supported by a recent paper on atomic relaxation in multilayer moiré heterostructures, which predicts very short (1 layer) penetration depth of moiré reconstruction for superlattices with λ < 20 nm [23].
However, to our surprise, when we cooled our devices to base temperature (30 mK) and measured Landau fan map, we observed a development of Hofstadter butterfly -the fractal quantum Hall effect (QHE), not just in the surface states, but across the entire bulk of the graphite film, Fig. 3.A low-T high-B Landau fan map of conductivity σxx versus n = nt + nb for device D2, presented in Fig. 3A, reveals multiple QHE features.Figure 3C summarises the results in a Wannier diagram.Analogous map and Wannier diagram for device D3 are plotted in Fig. 3B, D. Although QHE is in principle forbidden in a 3D electronic system like graphite, it has been recently reported for thin (up to 100 nm) graphite films [26].Two main factors contribute to the observed QHE: dimensional reduction of electronic system from 3D semimetal to 1D Landau bands in strong B-fields, and the consequent formation of standing waves in the 1D Landau bands due to finite thickness of graphite film.Standing waves result in the quantisation of 1D Landau bands and the development of minigaps among them, manifesting as a special 2.5D QHE.At large B-fields (above those marked by white dashed lines in Fig. 3A, B), only the two lowest Landau bands (0 and 1) remain crossing the Fermi level and contribute to the transport.These two bands are split by magnetic field δ10 ≈ 0.4 meV/T.In addition to being gapped by the standing waves, they are also spin resolved by Zeeman effect, 2μBB (μB is the Bohr magneton), while lifting their +KH and −KH valley degeneracy depends on the graphite layer parity.The Hofstadter butterfly [36] -a fractal set of energy eigenvalues for magnetic flux ϕ/ϕ0 = p/q commensurate with moiré periodic potential -is plotted in Fig. 4A for a honeycomb lattice [37], matching the geometry of our moiré perturbation.In charge transport measurements, this fractal set manifests in the Diophantine equation for the Landau fan and its Wannier diagram representation: Here integer t is the Landau filling index ν = nh/eB and integer s is the superlattice Bloch band filling index; n0 = 1/A0 is the charge carrier density per superlattice unit cell area.For s = 0, Eq. 1 plots conventional Landau fan diagram with t ≡ ν (grey lines in Fig. 3C, D), while for s,t ≠ 0, it traces Hofstadter states (purple lines in Fig. 3C, D) emanating from B-fields satisfying ϕ/ϕ0 = p/q.Figure 3E shows a clear hierarchy of QHE gaps: those measured at integer filling factors (s = 0) are an order of magnitude larger than the gaps between Hofstadter butterfly states (s,t ≠ 0).This suggests that the effect of moiré superlattice on the QHE can be considered as a small perturbation.To estimate the energy scale of this perturbation, we envelop 0 and 1 Landau bands in graphite with a Hofstadter butterfly energy spectrum.The spectrum of Landau levels for a 16-layer-thick graphite without moiré perturbation calculated following [26] (also see Methods) is plotted in Fig. 4B.It shows how, with increasing B field, levels cross each other, corresponding to closure and subsequent reopening of QHE gaps. Figure 4C plots the energy levels with the moiré envelope incorporated by giving each Landau level Em of the graphite QHE the fine structure of the Hofstadter butterfly ε: Here S ≈ 0.14 meV is the scaling factor for the bandwidth of Hofstadter butterfly [38], which was estimated from the measured transport gaps.The obtained enveloped spectrum shows good agreement with the experimental data as a function of filling factor, in terms of gap sizes and positions, as plotted in Fig. 4D.
We have shown that surface states of graphite (and, potentially, other semimetals) can be efficiently tuned by moiré superlattice potential resulting in Brown-Zak oscillations and Hofstadter surface states, stemming from aligned graphite and hBN interface.The alignment provides kaleidoscopic Lifshitz transitions, which develop into Brown-Zak oscillations and Hofstadter surface states.Remarkably, at low temperatures and high magnetic fields, the periodic surface potential leads to a bulk fractal Hofstadter butterfly which penetrates the entire graphite film.Hofstadter spectrum in this bulk twistronics system offers enticing possibility to explore bulk-boundary correspondence in the QHE regime.

Device fabrication
To make Hall bar devices, graphite flakes were encapsulated by hBN through dry transfer described elsewhere [39,40].Briefly, hBN was first mechanically exfoliated onto SiO2/Si substrates.The target hBN flake was then picked up by polymethylmethacrylate (PMMA) film and released onto a graphite flake exfoliated on SiO2/Si substrate with a partially overlapped area.To make aligned hBN/graphite structures, the straight edges of hBN and graphite flakes, which are usually along their crystallographic axes, are aligned in parallel.Standard e-beam lithography, reactive ion etching, and e-beam evaporation were applied sequentially to the hBN/graphite/hBN stack to pattern top gate electrode and metal contacts to graphite flake (3 nm Cr/80 nm Au).The device was 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.To make devices in Corbino geometry, we overexposed the resist to form crosslinked PMMA bridge, to separate inner contact, top gate, and the outer contact.
Graphite capacitors used to study surface states of non-aligned heterostructures, were fabricated similarly, but with quartz substrate replacing SiO2/Si wafer to minimise the parasitic capacitance.For capacitance devices, we used graphite flakes around 50 nm thick, guaranteeing the 3D graphite electronic spectrum.
Relatively thick hBN flakes (> 40 nm) were chosen to eliminate the inhomogeneity of electrostatic potential introduced by relatively rough metal electrode.

Transport and capacitance measurements
The longitudinal and Hall voltages of the Hall bar device were recorded with lock-in amplifiers (SR830 or MFLI) while applying a small low-frequency ac current of 10 nA (except where higher current is specified).
For Corbino devices, a small ac bias (40-100 μV) was applied to the inner contact, while current was recorded from the outer one.
Differential capacitance C was measured as a function of bias voltage Vb between the metal gate and graphite using an on-chip cryogenic bridge [41], which reaches sensitivity ≈10 aF at 1 mV excitation.In general, ac excitations of 102.53 kHz frequency and opposite phases were applied to the sample and a reference capacitor of known capacitance, respectively.Output signals from these two capacitors were mixed and applied to the gate of a high electron mobility transistor (HEMT), which served as an amplifier.
The excitation voltage of the reference capacitor was modulated so that the output signal from the HEMT became zero, and the capacitance of the sample can be deduced based on the ratio of their excitation voltage at the balanced point.The typical excitation voltage applied to the samples ranged from 1 to 10 mV, depending on the thickness of the hBN dielectric layer.

Tight-binding description of surface states
To describe surface states, we focus of the evanescent modes with complex momentum kz.To this end, we solve analytically the spectrum of graphite [32], considering the boundary conditions ( = 0 for surface carbon atoms).There are no complex kz solutions at zero doping, as only real kz solutions satisfying zero boundary conditions are normalizable.Electrostatic doping of graphite surface creates an inhomogeneous in z-direction potential near the surface, which does not preserve kz momentum, allowing real kz solutions near the surface, which then turn into evanescent modes decaying into the bulk.This provides a heuristic picture for the origin of non-trivial surface state solutions, Fig. S1D,E.
To quantify the surface states, we numerically find the self-consistent potential profile of doped graphite surface, considering a finite-thickness graphite sandwiched between two gates with carrier densities nt and nb.Since graphite screens the electric field at a depth of few layers, it is equivalent to a semi-infinite bulk graphite with one surface for N ≥ 10, where 2N is the number of graphene layers.In the Hartree approximation, the potentials on the layers, Uj > 1, are related with layer electronic densities as Where ϵ = 2.6 accounts for vertical polarizability of graphene [42], c = 3.35 Å is interlayer separation.
We temporarily fix the value of U1, which plays the role of a surface chemical potential, and then selfconsistently calculate Hartree potentials and densities on all the layers.The electronic density in layer  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  ⃗ .The constant n0 is chosen to match ∑   = 0 2 =1 to provide electrical neutrality.After finding the densities on all the layers we relate U1 with nt via   = −  − ∑   2 =1 .Numerically stable self-consistent solutions of Eq. 3 and 4 for different values of U1 can also be used to calculate the quantum capacitance through Eq. 5 in Supplementary Note 1.

Tight-binding description of surface states in aligned graphite
The spectrum for graphite aligned with hBN was calculated by treating the periodic moiré potential as a perturbation applied only to the top graphene layer, we followed the standard procedure [43], using the mirror-symmetric superlattice coupling Hamiltonian applied to the two top-layer components of the graphite film wavefunction, where Pauli matrices σ operate on top layer sublattices and τ operates on valleys,   = ℛ (−1)/3 {0, 4/(3)} are 6 reciprocal lattice vectors of superlattice (where ℛ is a rotation matrix) δ = 0.018 is a lattice mismatch  = 1.42 Å is C-C distance and we use the parameters U0 = 8.5 meV, U1 = -17 meV, U3 = -14.7 meV [44,45].The results do not significantly depend on the values of superlattice couplings, and it turned out to be sufficient to restrict the momentum space to the first star of superlattice reciprocal lattice vectors to achieve convergence.

Calculation of Δσxx in high temperature mappings
To highlight Brown-Zak oscillations across a large range of magnetic fields, we calculated Δσxx by subtracting a smooth background from σxx data.In comparison to graphene-hBN systems where the background conductivity can be fitted with polynomials [21], we find even high order (>10) polynomials insufficient as many oscillatory artefacts are present.Instead, we employ a two carrier Drude model of σxx(B), σxy(B) and fit both simultaneously to yield carrier densities and mobilities nh = 2.2 ×10 12 cm −2 , µh = 24,000 cm 2 V −1 s −1 , ne = 2.8 ×10 12 cm −2 , µe = −19,000 cm 2 V −1 s −1 for zero gate bias at T = 60 K.This two carrier model fit, σxx fit (B), is then used to calculate Δσxx(nb,B) = σxx(nb,B) -σxx fit (B).Oscillations observed in Δσxx were cross examined against σxx to confirm they are not artefacts.

Tight-binding description and standing waves in graphite in UQR
The allowed energy levels of graphite bulk in the absence of moiré perturbation were calculated following Ref.[26], and using the SWMC parameterisation of [31,32] to generate the zero-field spectrum of graphite.
The Bz-field leads to Landau quantization of the electron motion in plane and forming Landau bands.The UQR is defined as the regime where only the two lowest Landau bands (0 and 1) cross the Fermi level at kF ≈ π/4c, with interlayer separation c = 0.335 nm.In thin graphite films these 0 and 1 Landau bands are further discretised by the formation of standing waves in the kz direction, where the number of states these 1D bands split into (and, hence, the energy spacing between these states) is determined by the number of layers N. Further to this, Zeeman splitting, 2μBB, is applied to the resultant states, in order to complete the plot for 16-layer graphite in Fig. 4B.

Analysis of Hofstadter spectrum and Brown-Zak oscillations
To model the conductivity of our aligned devices at high B field and low temperatures we consider the Moiré superlattice potential as a weak perturbation; each Landau-level-like energy level described above (and plotted in Fig. 4B) is split into q subbands, at a given ϕ/ϕ0 = p/q.At rational fields ϕ/ϕ0 = p/q translational symmetry is restored if we consider a magnetic unit cell of q supercells, and the q×q tight binding Hamiltonian is defined by the difference equation [37] ( where m is an integer,  =   √3/2, and the calculation has been simplified by using a normalised hopping parameter γ = 1, which yields an energy range of ε = ±3 (arb.units).Solving for energy eigenvalues generates the honeycomb lattice Hofstadter butterfly (Fig. 4A).Note that near significant flux fractions ϕ/ϕ0 = 1, 1/2, 1/3 the density of points increases, and much larger q are required to populate these regions, so to keep computation time reasonable the calculations are limited to q ≤ 100.
In order to apply this fine structure of Hofstadter butterfly to each energy level we need a relevant energy scaling factor, defined as S in the main text.Analysis of the thermal activation of gaps for device D3 at B = 13.5 T (Fig. 3E) indicate the largest fractal gaps are ≈0.1 meV (consistent between Corbino devices, see Supplementary Note 3 and Fig. S4).We then attribute this to the largest gap in the Hofstadter spectrum at the flux value ϕ/ϕ0 = 0.56 (corresponding to B = 13.5 T, dashed line in Fig. 4A), which occurs at 1.0 < ε < 1.7 (or -1.7 < ε < -1.0 due to Hofstadter butterfly symmetry) therefore resulting in S = 0.14 meV.The full spectrum is then calculated using Eq. 2 and plotted in Fig. 4C where we use the periodicity of Hofstadter butterfly (such that ε(ϕ/ϕ0 + ρ) = ε(ϕ/ϕ0) for any integer ρ) to plot states at ϕ/ϕ0 > 1.The result is not only the presence of fractal fine structure in the spectrum but also the reduction in size of the main QHE gaps, due to Hofstadter butterfly effectively broadening each energy level by a finite width of 6S.Again, the computation limit of q ≤ 100 results in a blank region at ϕ/ϕ0 = 1, which are in fact the densest regions.The overlap of these dense regions at ϕ/ϕ0 = 1 completely closes some of the main QHE gaps (e.g., ν = 0 and 4) which is in good agreement with the measured conductivity (Fig. 4D).To provide qualitative understanding of the surface states in graphite, we plot the eigenstates for homogenous bulk graphite in Fig. S1D, E, which consist of propagating (black curves, real kz) and evanescent (orange curves, complex kz) modes.There are three types of propagating modes: majority electron and hole bands with bandwidth 2γ2 and a minority carrier band near ckz = π/2, γ2 is the hopping amplitude between two non-dimer sites in the next-nearest layers.These propagating bands cross the bulk Fermi level for small in-plane momentum kx,ky (Fig. S1D, 1D metal regime in z-direction), but are spread away from the Fermi level at large kx,ky (Fig. S1E, 1D semiconductor regime in z-direction).When a potential near the surface is introduced by doping, the propagating modes in 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, as shown by red/blue arrows in Fig. S1E for hole/electron doping, respectively.The dispersion of these evanescent modes, which we denote as Type 1, crosses the Fermi level, forming a surface Fermi-line with radius larger than Fermi surface of propagating carriers, as illustrated by the yellow contour in Fig. S1C.These states are similar to surface states in doped semiconductor, with the difference that they exist only for in-plane momenta larger than the projection of bulk Fermi surface of graphite.In 1D metal regime, another type of evanescent mode, denoted as Type 2, appears for |E| > |γ2|and never crosses the Fermi level, Fig. S1D.
Returning to the Cq-Us curve in Fig. S1B, at low doping when |Us| < |γ2| (that is, |n| <  The methods of energy spectrum modelling applied to device D3 in the main text (Fig. 4) were also applied to device D2, which has a different alignment to hBN and layer parity.To be more specific, device D2 is 21layer in thickness and only aligned to one encapsulating hBN.Layer parity plays a key role in QHE in graphite.The QHE exists as a result of standing waves in kz Landau bands, and the maxima of these standing waves is in the +KH valley for even layers and -KH valley for odd layers [26].For even total number of layers (N = 2 , i.e., device D3), there are equal numbers of odd and even layers to host the states in +KH and -KH valleys, states therefore have a twofold degeneracy (∆ν = 2).However, in odd layer parity graphite films (N = 2-1, i.e., device D2), the mismatch between number of even and odd layers lifts this degeneracy and the energy gaps between states are halved.
As shown in Fig. S4, device D2 exhibits main QHE sequence gaps (Fig. S4A-D) that are significantly reduced compared to that in D3 (Fig. 3E), due to the lifted ±KH valley degeneracy (Fig. S4E) because of the odd-layer parity of graphite in D2.In Fig. S4D, we focus on the zeroth gap size measured as a function of magnetic field between two level crossings at B = 10 T and B = 16 T, with a maximal observed gap of ≈0.48 meV.Fig. S4F plots the energy levels as a function of B for 21-layer graphite, calculated in the same way as in Ref. [26].Without moiré perturbation, both the extent of the zeroth gap (8.5 T < B < 17 T) and its maximal size (1.3 meV) in the model are notably larger than in the experiment (cf.similar maximal gap size ≈1.1 meV calculated for 21 layers in Ref. [26]).Upon application of Hofstadter butterfly to each energy level (see methods) using the same scaling factor S = 0.14 meV as in Fig. 3E in the main text, the zeroth gap in the model is reduced to ≈0.6 meV (Fig. S4G), in closer agreement with our experiments.However, such effective broadening of energy levels from the Hofstadter butterfly will lead to many overlapping states and hence gap closures, which were not observed in experiment.We postulate that this is due to inadequate treatment of equal moiré perturbation to states in both even and odd layers (+KH and -KH valleys).As only one of the  odd layers is directly affected by the Moiré potential, which rapidly decreases in strength with gaining depth into graphite, the perturbation strength on even layers is expected to be significantly reduced.In Fig. S4H, we plot a revised model with S = 0.14 meV for odd layers and S = 0.04 meV for even layers, yielding less main QHE gap closures than Fig. S4G while the same zeroth gap remains.Qualitatively, this is in a closer agreement with the experiment.Future work to further improve this model may use a more detailed treatment to Landau level mixing.Experimental results are better visualized and more informative when presented as C(n,B) map, Fig. S5C.Note the branches of surface states spawn out from the neutrality point at B = 7.5 T, 3 T, 2 T, and so on.These B-fields correspond to the critical fields above which the BLBs no longer cross the Fermi level and appear only as SLLs.For instance, according to our SWMC model, at 7.5 T, BLB 2 + is just above the Fermi level (Fig. S5D).Thus, a branch of surface states spawned out around this field is labelled as S 2+ .The same happens with the electron BLB 3 + at 3 T and hole BLB 2 − at 2T, Fig. S5C.Note that above 7.5 T, only BLB 0 and 1 cross the Fermi level at charge neutrality point and graphite enters the ultraquantum regime (UQR) [48,49].The broad feature in the middle region of the map, whose position does not change notably with B is also attributed to the BLB 0, 1 [50,51].
We observed oscillations down to B ≈ 0.1 T, see Fig. S5C and Fig. S6A, which sets a lower bound for surface charge carrier mobility of ≈ 100,000 cm 2 V -1 s -1 .The high electronic quality of surface states also enables fractional features in Landau quantization of charge carriers.Another graphite capacitor device was therefore fabricated to investigate this, with a thicker hBN dielectric to reduce the inhomogeneity of electrostatic potential from the metal electrode.At high magnetic field, B = 20 T, we observe the formation of two minima on top of singly degenerate states of S 2+ (Fig. S7A).We plotted dC/dv as a function of filling factor ν at high field region (Fig. S7B).The Δν between the fractional gap is around 0.27, which is lower than expected Δν = 1/3 for fractional QHE.
To further characterize these fractional states, a graphite Hall bar device with a 6 nm thick graphite was fabricated, with graphite not aligned with hBN.In this thin graphite device, the formation of standing waves leads to bulk QHE gaps in the UQR [26].Through applying a displacement field, B-n regions can be found where the energy level of surface states locates within the bulk gap (Fig. S8).In these regions, the surface states are isolated from the bulk completely, and both vanishing longitudinal conductivity (σxx) and quantized Hall conductivity (σxy) clearly indicate the development of fractional QHE with a 1/3 degeneracy.The difference between the capacitance and transport measurements can be reconciled, when considering the negative compressibility of the fractional states: chemical potential of the surface states reduces with the injection of n [52][53][54], acquiring additional charges from the bulk.Not only S 2+ , but also S 2− shows a similar behaviour (Fig. S9), indicating a strong electron-electron interaction in the surface states.

Fig. 1 .
Fig. 1.Moiré superlattice at the interface between graphite and hexagonal boron nitride.(A) diagram of a typical heterostructure stack of graphite encapsulated in hBN with one of the interfaces aligned.Here the lattice mismatch between Grt and hBN has been exaggerated for clarity.(B, C) Optical micrographs of devices D1 (Hall bar geometry with 8.2-nm-thick graphite channel, bottom hBN aligned to graphite) and D3 (Corbino geometry with 5-nm-thick graphite channel, both top and bottom hBN flakes aligned to graphite).Scale bars are 10 μm.(D) Conductivities σxx and σxy as a function of carrier density induced by the bottom gate, nb, for device D1 at the aligned interface and non-aligned device D4, measured at T = 0.24 K and non-quantizing magnetic fields (B = 120 mT).(E) Line cuts through the dispersion relation of the SBZ in the kx ky plane, at carrier densities (bottom to top) n = -3.8,-3.6, -2.1, -2.0, 1.9, 2.3, 3.6, 3.9 × 10 12 cm -2 grouped as pairs [labels a,b,c,d correspond to the regions highlighted in (D)].Black dashed hexagon denotes the boundary of the 1 st SBZ and red (blue) denotes holes (electron) line cuts, some lines at the corners are extended into 2 nd SBZ for clarity.

Fig. 2 .
Fig. 2. Brown-Zak oscillations arising from states confined to the graphite-hBN interface.(A) Conductivity σxx as a function of B for device D2 (Corbino geometry with 7-nm-thick graphite channel) at high electron doping electrostatically induced by top (bottom) gate, which tunes the aligned (non-aligned) graphite-hBN interface.When tuning the aligned interface, peaks are evident at values of B equivalent to flux quantum 1/q per superlattice unit cell.Carrier density nt = 3.1 × 10 12 cm -2 and nb = 3.3 × 10 12 cm -2 for aligned and non-aligned interfaces, respectively.(B) Conductivity map (a smooth background subtracted, see Methods) as a function

Fig. 3 .
Fig. 3. Fractal QHE states in graphite.(A, B) Conductivity σxx as a function of n = nt + nb and B for devices D2 and D3 measured at T = 30 mK and nt = nb (i.e., displacement field D = 0).The white dashed lines denote transition from surface Landau levels to bulk ultraquantum regime (UQR).(C, D) Associated Wannier diagrams plotting the QHE (grey) and fractal QHE (purple) states in the UQR observed in (A) and (B).Horizontal axis is normalised in units of electron per supercell, with n0 = 1/A0.Below UQR, orange (brown) lines trace fractal (non-fractal) states in surface Landau levels +2 and -2.(E) Hierarchy of QHE gaps in aligned hBN/graphite/hBN heterostructure.Lower panel shows σxx(n) traces at different temperatures for device D3 at B = 13.5 T, which is used to extract gap sizes from Arrhenius activation.Upper panel is a bubble plot of QHE gaps; the area of a circle scales linearly with the gap size (ranging from 30 μeV to 1.8 meV).Colour coding is the same as in (D) and labels are integers s and t from Eq. 1, for the QHE (t only) and fractal QHE (s,t) states.

Fig. 4 .
Fig. 4. Hofstadter broadening of energy levels in graphite.(A)Hofstadter butterfly calculated for a honeycomb lattice following Ref.[37], where for simplicity the hopping parameter γ = 1.The dashed line denotes a flux ratio of ϕ/ϕ0 equivalent to B = 13.5 T in device D3, which is the field strength as in Fig.3E.(B) Allowed energy levels resulting from quantised states from 0 (1) Landau bands shown in red (grey), calculated for 16-layer-thick graphite film without a moiré perturbation[26].Zeeman splitting has been included, as indicated by light and dark lines for spin up and spin down, respectively.Integer labels in blue refer to filling factor ν. (C) Combination of panels (A) and (B) by applying the Hofstadter butterfly as a small perturbation to each energy level, Eq. 2. Labelled ν same as in (B).(D) Conductivity map for device D3, similar to that in Fig.3Bexcept plotted as a function of filling factor of main QHE sequence.Colour scale brown to yellow is from 0 to 247 μS.

Fig. S1 .
Fig. S1.DoS spectroscopy of surface states of graphite in zero magnetic field.(A) Capacitance (C) vs carrier density (n) at T = 0.3 K. Inset shows optical micrograph of one capacitor device (device D5), scale bar 20 µm.(B) Density of states (DoS) and quantum capacitance (Cq) vs surface potential (Us) from calculations based on effective-mass model (black line).Coloured symbols are experimental data from four different devices.Inset is calculated Us vs n. (C) Calculated dispersion of a 14-layer-thick graphite film for electron doping n = 6 • 10 12 cm -2 as a function of in-plane momentum kx,ky (horizontal axes) and energy E (vertical axis).Red (green) colour indicates surface states having high probability density of the wave function at the first (second) graphene bilayer.Upper outer surfaces with larger radius correspond to Type 1, and lower inner surface with a smaller radius correspond to Type 2. Blue colour indicates bulk states, and yellow contour highlights the Fermi level.(D, E) Dispersion for propagating (Im kz = 0, black lines) and evanescent modes (Im kz ≠ 0, orange lines) for bulk graphite as a function of complex kz for fixed ħνk = 0.04 eV (panel (D), 1D metal regime with Type 2 evanescent modes) and 0.15 eV (panel (E), 1D semiconductor regime with Type 1 evanescent modes), respectively.Fermi level is at 0 eV.Green/blue arrows indicate evolution of surface states from propagating modes into evanescent modes for electron/hole doping.

Fig. S4 .
Fig. S4.Gap size hierarchy in device D2. (A) Upper panel, σxx as a function of integer QHE filling factor ν at various temperatures, B = 15.2 T. Lower panel, bubble plot of extracted gap size from Arrhenius fits to the σxx minima.Gap size scales with area (60 µeV to 0.9 meV), and grey bubbles are integer gaps (labelled by ν) and purple bubbles are fractal gaps (labelled by s, t).(B-C) Examples of Arrhenius plots of ln[σxx] as a function of reciprocal temperature for integer QHE gap at ν =5 and fractal QHE gap of integers (s,t)=(-1,7), respectively (these two states are highlighted by red circles in panel A).The linear regions are fitted to yield a slope of ~½Egap.(D) Gap size extracted from Arrhenius fits as a function of B for the zeroth gap.(E) Conductivity map for device D2, same data as in Fig.3Ain the main text, except plotted as a function of filling factor ν of main QHE sequence.Colour scale brown to yellow is 0 to 210 μS.(F) Allowed energy levels resulting from quantised states from 0 (1) Landau bands shown in red (grey), calculated for 21-layer-thick graphite film without a moiré perturbation[26].Zeeman splitting has been included, as indicated by light and dark lines for spin up and spin down, respectively.(G) Combination of panels (F) and Fig.4Aby applying the Hofstadter butterfly as a small perturbation (S = 0.14 meV) to each energy level, Eq. 2. (H) Same as (G) but with S = 0.14 meV for odd layer states and S = 0.04 meV for even layer states.

Fig. S6 .
Fig. S6.Low field quantum oscillations.(A) Oscillations of Δρxx in device D7 obtained by subtracting a smooth background while sweeping the gate voltage (Vb) at 0.1 T, 0.3 K using excitation I = 1 μA.The encapsulated graphite with a thickness of 20 nm is defined to Hall bar geometry for this measurement.(B) dns/dn as a function of n at 0.6T, where ns is the carrier density injected into the surface states, n is the total electrostatically induced carrier density.It is deduced from the curve shown in Fig. S5A, based on the periodicity of the oscillations.

Fig. S8 .
Fig. S8.Fractional surface states in transport measurements.(A) Longitudinal conductivity σxx (B, n) measured at D = (nt-nb)•e/20 = 0.24 V/nm in a Hall bar device D9 (9-nm-thick-graphite), T = 0.3 K, I = 20 nA.white shaded areas are guides to the eye for the surface states.Boundaries of one such surface states are marked by white dotted curves.Logarithmic colour scale, navy to orange is 0.1 to 118.2 µS.(B) σxx cut profile (black curve) of the white dashed line in (A) and the corresponding Hall conductivity σxy (red curve).
6 • 10 11 cm -2 ), there are no surface states, and quantum capacitance is determined by electron and hole screening.Since holes have a slightly larger DoS (shallower dispersion) than electrons, we see larger CQ at hole doping.When doping reaches Us ≈ ±γ2, Type 1 and Type 2 surface states appear and contribute to quantum capacitance.The radius of surface Fermi-line for Type 1 states, grows with |n|, leading to increasing density of surface states and growth of quantum capacitance.