Mixing of moiré

Nature (2023)Cite this article 5061 Accesses 120 Altmetric Metrics details Van der Waals assembly enables the design of electronic states in two-dimensional (2D) materials, often by superimposing a

Nature (2023)Cite this article

5061 Accesses

120 Altmetric

Metrics details

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,2,3,4,5,6,7,8,9. This twistronics approach has resulted in numerous previously undescribed physics, including strong correlations and superconductivity in twisted bilayer graphene10,11,12, resonant excitons, charge ordering and Wigner crystallization in transition-metal chalcogenide moiré structures13,14,15,16,17,18 and Hofstadter’s butterfly spectra and Brown–Zak quantum oscillations in graphene superlattices19,20,21,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, its periodic lattice is interrupted, and surface states arise with wavefunctions exponentially decaying into the bulk of the crystal25. For example, surface charge accumulation in semiconductors leads to distinct 2D subbands tunable by electrostatic gating. By contrast, 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 such as bismuth and graphite, which have tunable surface states that are interesting but remain underexplored. Graphite films are of interest as they show both 3D and 2D electronic properties controlled by electrical doping and an external magnetic field B. Notably, graphite of a finite thickness exhibits an unusual 2.5-dimensional (2.5D) quantum Hall effect (QHE)26.

In this Article, we explore moiré engineering of highly tunable electronic states, by aligning two bulk crystals, hexagonal graphite and hexagonal boron nitride (hBN). To this end, we prepared hBN/graphite/hBN heterostructures by aligning thin graphite films (about 5–10 nm thick) on top of the hBN substrate and encapsulating the stack with another hBN crystal. Unless otherwise stated, this latter, encapsulating, hBN is intentionally misaligned (see Methods, ‘Device fabrication’ for details). As the lattice constants of hBN and graphite are close, in the 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 the moiré superlattice, the hBN encapsulation also preserves the high electronic quality of graphite films26,27,28. Figure 1a–c shows schematics and micrographs of the hBN/graphite/hBN heterostructures, fabricated into Hall bar and Corbino geometry devices. In these devices, the top and bottom electrostatic gates were used to independently control carrier densities nt and nb, at the top and bottom interfaces of the hBN/graphite/hBN heterostructure. In total, we have studied 11 graphite heterostructure devices (Extended Data Table 1).

a, Schematic of a heterostructure device with graphite (labelled Grt) encapsulated in hBN with one of the interfaces aligned. Here the lattice mismatch between graphite and hBN has been exaggerated for clarity. b,c, Optical micrographs of devices D1 (b) and D3 (c). Scale bar, 10 μm (b and c). d, Conductivities σxx and σxy as a function of the carrier density induced by the bottom gate, nb, for aligned device D1 and non-aligned device D4, measured at T = 0.24 K and non-quantizing B = 120 mT. e, Line cuts through the calculated dispersion relation in the kx–ky plane of the SBZ, at carrier densities (bottom to top) n (×1012 cm−2) = −3.8, −3.6, −2.1, −2.0, 1.9, 2.3, 3.6 and 3.9, grouped as pairs. Labels A, B, C and D correspond to the regions highlighted in d. The black dashed hexagon denotes the boundary of the first SBZ and red curves denote the hole and blue curves denote electron Fermi-surface cuts. Some lines at the corners are extended into the second SBZ for clarity.

Hexagonal graphite (Bernal stacking) is a compensated semimetal with the Fermi surface occupying only a small fraction of the Brillouin zone. The size of the Fermi surface is determined mostly by a through-layer hopping parameter, γ2 ≈ −20 meV (ref. 29). Owing to its semimetallic nature, graphite does not host surface states (evanescent modes) in the absence of dangling bonds or applied electric field. However, if an electric field above a certain value is applied perpendicular to the basal plane, tunable surface states emerge26,30 (see Methods, ‘Surface states in non-aligned graphite films in a zero-B field’ and Extended Data Fig. 1).

We find that a moiré superlattice at the surface of graphite markedly modifies its surface states, resulting in an entirely different transport behaviour observed between aligned and non-aligned devices (Fig. 1d). Devices with a non-aligned interface show a nearly featureless carrier-density dependence of longitudinal σxx(n) and transversal σxy(n) conductivities in small B. By contrast, for the aligned graphite interface, σxy(n) shows multiple zero crossings that are accompanied by peaks in σxx(n). We attribute this behaviour to the recurrence of electrostatically induced surface states occupied by electron- or hole-like charge carriers. To quantify this, we calculated Fermi-surface projections using an effective-mass model with Slonczewski–Weiss–McClure parameterization of graphite31,32 subjected to moiré superlattice potential, in combination with self-consistent Hartree analysis (see Methods, ‘Surface states in graphite films in the presence of a moiré superlattice’). Our calculations in the superlattice Brillouin zone (SBZ) of an hBN/graphite/hBN heterostructure show a multitude of surface states with numerous topological Lifshitz transitions (LTs) across a range of carrier densities (Fig. 1e). The four pairs of plots (labelled A–D in Fig. 1e) with considerable changes in the Fermi-surface topology demonstrate four LTs that correspond to the four n ranges (Fig. 1d). LTs observed at |n| ≈ 2.0 and 3.7 × 1012 cm−2 belong to two different branches of the surface states—one residing mostly on the first graphene bilayer of graphite and the other mostly on the second bilayer (Extended Data Fig. 1c–f). As B increases, the surface states in the vicinity of LTs give rise to separate branches of Landau levels. For details of evolution of σxx(n) and σxy(n) in low magnetic fields see Methods, ‘Surface states in graphite films in the presence of a moiré superlattice’ and Extended Data Fig. 2) Extended Data Fig. 2e,f provides a further comparison of aligned and non-aligned interfaces of device D1, confirming the absence of LTs in surface states at the non-aligned interface.

With high B, the difference between hBN/graphite/hBN devices with aligned and non-aligned interfaces becomes even more prominent (Fig. 2a). The curves σxx(B) were measured at 60 K to suppress Landau quantization. If the aligned surface is doped away from the electron–hole compensation, σxx shows an oscillatory behaviour periodic in 1/B. Peaks in σxx appear in fields \({B}_{1/q}=\frac{1}{q}\frac{{\phi }_{0}}{{A}_{0}}\), corresponding to the integer number q of superlattice unit cells with an area A0 = √3/2λ2 that are commensurate with the magnetic flux quantum ϕ0 = h/e, where λ is the wavelength of moiré superlattice, h is the Planck’s constant and e is the elementary charge. The commensurability between ϕ0 and the magnetic flux through a moiré unit cell, ϕ = BA0, can be interpreted as a manifestation of Brown–Zak quantum oscillations at the superlattice interface, which were recently reported for aligned monolayer graphene/hBN heterostructures21,22. The formation of magnetic Bloch states leads to higher conductivity because of the straight rather than cyclotron trajectories of the surface-charge carriers21,22,33,34, as evidenced by the conductivity peaks at B1/q (Fig. 2a). Figure 2b shows some of these nb-independent conductivity peaks that were found at all distinguishable 1/q-commensurate fields. Note that not only unit fractions but also second-order fractal states (for example, B2/5 in Extended Data Fig. 3a) can be seen in Brown–Zak oscillations.

a, Conductivity σxx as a function of B for device D2 at high electron concentrations induced by the top and bottom gates that dope the aligned and non-aligned graphite–hBN interfaces, respectively. For doping of the aligned interface, peaks appear at values of B equivalent to one flux quantum per q superlattice unit cells. nt = 3.1 × 1012 cm−2 and nb = 3.1 × 1012 cm−2 for aligned and non-aligned interfaces, respectively. T = 60 K. b, Conductivity map (a smooth background subtracted, see Methods, ‘Surface states in graphite films in the presence of a moiré superlattice’) as a function of B and nb at the aligned graphite–hBN interface of device D1. Measurements were performed at 60 K to suppress Landau quantization. The right y axis denotes the inverse flux ϕ0/ϕ. c, σxx(nb,B) map for the same device at 20 K.

Because Brown–Zak oscillations stem from the translational invariance of magnetic Bloch states at rational fractions of magnetic flux ϕ/ϕ0 = p/q (where p is an integer), they are insensitive to temperature as long as electrons retain phase coherence in the area of the magnetic supercell qA0. Figure 2c shows that at intermediate temperatures (T = 20 K), states with ϕ0/ϕ up to 24 are visible (and even states with ϕ0/ϕ > 35 are distinguishable in Extended Data Fig. 3b). This provides a lower bound on the phase coherence length of greater than about 100 nm. Brown–Zak oscillations can also be interpreted as Aharonov–Bohm interference in a periodic 2D network formed by classic trajectories of electrons drifting around the Fermi contours that are joined by magnetic breakdown tunnelling in the vicinity of Van Hove singularities (see Methods, ‘Conventional interpretation of Brown–Zak oscillations’ and Extended Data Fig. 4). This interpretation enables a convenient conceptual transition into the regime of low-B fields in which we see multiple LTs of the Fermi-surface topology (Fig. 1e) and explains the disappearance of Brown–Zak oscillations for |nb| < 2 × 1012 cm−2.

In comparison, no LTs or Brown–Zak oscillations could be observed in our hBN/graphite/hBN devices if non-aligned interfaces were gated (Figs. 1d and 2a and Extended Data Fig. 2). This is not surprising, as it has previously been shown that states at the opposite surfaces of a graphite film are well screened from each other, with a screening depth of only two to three layers26. Raman measurements also do not show any qualitative difference in strain distribution or other effects of alignment for films thicker than seven to eight graphene layers at both aligned and non-aligned graphite interfaces (see Methods, ‘Raman spectroscopy of aligned graphite films’ and Extended Data Fig. 10). This conclusion is further supported by a recent report on atomic relaxation in multilayer moiré heterostructures23 that predicts a very short (one layer) penetration depth for moiré reconstruction with superlattice periodicity λ < 20 nm.

Surprisingly, if the aligned devices are cooled down to our lowest T of 30 mK and Landau fan maps are measured, we observe the development of Hofstadter’s butterfly—the fractal QHE—not just in the near-surface 2D states but across the entire graphite film (Fig. 3) as witnessed by gating either bottom (non-aligned) or top (aligned) interface. A high-B map of conductivity σxx versus n = nt + nb for device D2 (Fig. 3a) shows multiple QHE features. Figure 3c traces the observed conductivity minima on a Wannier diagram. An analogous map and Wannier diagram are presented in Fig. 3b,d for device D3. Although QHE is forbidden in 3D electronic systems, it has recently been reported for thin (up to 100 nm) graphite films26. Two main factors contribute to the observed QHE: dimensional reduction of the electronic system from a 3D semimetal to one-dimensional (1D) Landau bands in strong B and the consequent formation of standing waves in the 1D Landau bands because of a finite thickness of graphite films. Standing waves result in the quantization of the 1D Landau bands and the development of minigaps, which manifest themselves in a so-called 2.5D QHE. At high fields (above B marked by white dashed lines in Fig. 3a,b), only the two lowest Landau bands (0 and 1) cross the Fermi level and contribute to magnetotransport. In addition to being gapped by the standing waves, these two bands are split by an energy gap δ10 ≈ 0.4 meV T−1, and further spin-resolved by the Zeeman gap, 2μBB (μB is the Bohr magneton). Lifting the +KH and −KH valley degeneracy of these bands depends on the graphite layer parity26.

a,b, Conductivity σxx as a function of n = nt + nb and B for devices D2 (a) and D3 (b), T = 30 mK and nt = nb (that is, zero displacement field). The white dashed curves indicate the transition from surface Landau levels to the bulk UQR. Black arrows point to the threshold filling factors that bound the region of bulk in which fractal states are observed (ν = −9 and 12). c,d, Associated Wannier diagrams for panels a (c) and b (d): 2D QHE (grey) and fractal 2.5D QHE (purple) states in the UQR. The x axis is in units of n0 = 1/A0. Below UQR, orange lines trace fractal states and brown lines trace non-fractal states in the surface Landau levels +2 and −2. e, Hierarchy of 2.5D QHE gaps in aligned hBN/graphite/hBN. Bottom, σxx(n) traces at different T for device D3 at B = 13.5 T, which are used to extract gap sizes from Arrhenius activation. Top, bubble plot of the QHE gaps in which the area of circles scales linearly with the found gap sizes (ranging from 30 μeV to 1.8 meV). Grey and purple colour coding is the same as in d and labels are integers s and t from equation (1) for standard QHE (t only) and fractal QHE (s, t) states.

To show how the QHE states of Hofstadter’s butterfly penetrate through the entire graphite bulk, we have also measured σxx as a function of n using both top and bottom gates at fixed B (see Extended Data Fig. 5a,c and the corresponding Wannier diagrams in Extended Data Fig. 5b,d). 2.5D QHE gaps appear as diagonal features because these Landau levels can be filled equivalently by either nb or nt and, therefore, the states extend throughout the bulk. This is the case for both standard QHE and Hofstadter’s butterfly gaps, which shows that, in the ultraquantum regime (UQR), the moiré surface potential affects the entire bulk of graphite. Conductivity maps σxx(nt, nb) at the high-B field for doubly aligned device D3 is generally similar to that of singly aligned device D2, with QHE and Hofstadter’s butterfly gaps following both nt and nb (Extended Data Fig. 5). One notable difference attributed to the symmetry breaking of the interface alignment is that σxx(nt, nb) is asymmetric for singly aligned device D2, but symmetric for doubly aligned device D3 (Extended Data Fig. 5).

Hofstadter’s butterfly35—a fractal set of energy eigenvalues for magnetic fluxes ϕ/ϕ0 = p/q—is shown in Fig. 4a for a honeycomb lattice36, matching the geometry of our moiré superlattice. In electron transport measurements, this fractal pattern manifests itself in a Landau fan diagram and its Wannier representation, which are described by the Diophantine equation:

where integer t is the Landau filling factor ν = nh/eB and integer s is the superlattice Bloch band-filling index; n0 = 1/A0 is the density of one electron per superlattice unit cell. For s = 0, equation (1) corresponds to the conventional Landau fan with t ≡ ν (Fig. 3c,d, grey lines), whereas for s ≠ 0, it traces Hofstadter states (Fig. 3c,d, purple lines) emanating from magnetic fields satisfying ϕ/ϕ0 = p/q.

a, Hofstadter’s butterfly calculated for a honeycomb lattice following ref. 36, with a normalized energy scale. The dashed line marks ϕ/ϕ0 equivalent to B = 13.5 T, the field strength as in Fig. 3e. b, Landau levels resulting from quantized states from 0 Landau bands are shown in red and 1 Landau bands are shown in grey and calculated for a 16-layer-thick graphite film without a moiré perturbation26. Zeeman splitting is included, as indicated by lighter and darker curves for the spin up and down, respectively. Labels in black refer to the filling factor ν. c, Expected spectra by applying Hofstadter’s butterfly in a as a small perturbation to each Landau level in b. Same labelling as in b. d, Conductivity map replotted from Fig. 3b as a function of ν. Two prominent gap closures around ϕ/ϕ0 = 1 have been labelled by the Landau band origin and spin of the levels that are crossing; the same colour coding as in b,c. Colour scale: brown to yellow, 0 μS to 115 μS.

Figure 3e shows a hierarchy in the observed QHE gaps: those measured at integer filling factors (s = 0) are one order of magnitude larger than the gaps because of Hofstadter’s butterfly (s ≠ 0). This suggests that the effect of the moiré superlattice on the QHE can be considered a small perturbation. To model the impact of this perturbation, we envelop the standing waves of 0 and 1 Landau bands in graphite with Hofstadter’s butterfly energy spectrum. Figure 4b shows the Landau-level spectrum calculated for a 16-layer-thick graphite film without taking into account moiré perturbation (see Methods, ‘Bulk states in graphite films in the presence of a surface moiré superlattice’), where Landau levels cross each other with increasing B, which corresponds to closure and subsequent reopening of gaps in the 2.5D QHE. Figure 4c plots the same 16-layer-graphite spectrum but each Landau level, Em, is now augmented by the Hofstadter’s butterfly spectrum, ε, using

Here S ≈ 0.42 meV is the scaling factor for the bandwidth of Hofstadter’s butterfly37, which was estimated from the measured transport gaps. The obtained \({E}_{{\rm{m}}}^{{\rm{moir \acute{{\rm{e}}} }}}\) spectrum shows good agreement with our experimental data in terms of both sizes and positions of the gaps (Fig. 4d) (gap closures at ϕ/ϕ0 = 1 are labelled with crossings of the corresponding \({E}_{{\rm{m}}}^{{\rm{moir \acute{{\rm{e}}} }}}\) states).

Bulk fractal states observed here are different from those of graphene-based 2D electronic systems, and our mixed moiré system demonstrates plenty of additional, non-trivial physics, inaccessible in 2D systems. First, the conductivity of bulk graphite can be efficiently tuned using interface alignment—σxx is increased more than two times at zero B field and up to an order of magnitude at high B in aligned versus non-aligned devices (Extended Data Fig. 4c,d). Second, B-field dependence of the amplitude of Brown–Zak oscillations observed in aligned graphite films differs from that of graphene, showing a non-monotonic dependence of the amplitude of oscillations on 1/B and emphasizing the richness of our 3D twistronics system (Extended Data Fig. 4e). Third, the absence of bulk fractal states beyond threshold filling factors (ν = −9 and 12; Fig. 3a,b and Extended Data Fig. 5) shows that the mixing of moiré surface and bulk states can be controlled electrostatically by changing the screening of moiré potential to the bulk bands.

In summary, we have shown that surface states in graphite (and, potentially, other semimetals and doped semiconductors) can be strongly modified by a moiré superlattice potential. The alignment between hBN and graphite provides a kaleidoscope of LTs that develop into Brown–Zak oscillations and Hofstadter surface states. Remarkably, moiré surface states in high magnetic fields also affect the entire electronic spectrum of these graphite films, which results in fractal Hofstadter butterflies that can be referred to as 2.5D in analogy with the 2.5D QHE in graphite. Our approach thus offers a possibility to explore mixed-dimensionality effects arising because of surface superlattices extending their influence deep into the bulk of 3D electronic systems.

To make hBN/graphite/hBN heterostructures, graphite flakes were encapsulated by hBN through dry transfer as described elsewhere38,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 SiO2/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 SiO2/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.

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(ro/ri)G, where G is the measured conductance and ro is the outer radius and ri 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 ϕ/ϕ0 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 B0 for each sequence of oscillations, where B0 is the magnetic field at which ϕ0 = B0A0 with magnetic flux quantum ϕ0 and superlattice unit cell area A0. Using this value of B0, the moiré wavelength λ is calculated as λ = \(\sqrt{2{A}_{0}/\sqrt{3}}\). Furthermore, B0 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 Vb between the metal gate and graphite using an on-chip cryogenic bridge40, 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.

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

where ε = 2.6 accounts for the vertical polarizability of graphene41, c = 3.35 Å is the interlayer separation and 2N is the number of graphene layers. We temporarily fix the value of U1, 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 n0 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 U1 to nt 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 systems40,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 CG = εε0A/d and quantum capacitance CQ in series, 1/C = 1/CG + 1/CQ, where A is the device area, ε0 is the vacuum permittivity and d and ε are the thickness and relative permittivity of the hBN dielectric layer43. The quantum capacitance reflects the DOS = dn/dU1 on the surface of graphite: CQ = Ae2dn/dU1, 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 (γ0 = 3.16 eV, γ1 = 0.39 eV, γ2 = −17 meV, γ3 = −0.315 eV, γ4 = 44 meV, γ5 = 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 (|U1| < |γ2|, that is, |n| < 6 × 1011 cm−2), 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 CQ at hole doping. When doping reaches U1 ≈ ±γ2, 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 CQ. 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 × 1012 cm−2 to 6 × 1012 cm−2, 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 graphite32, 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 kz) and an evanescent mode (orange curves, complex kz). There are no complex kz solutions at zero doping, as only real kz 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 kz momentum, allowing real kz 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γ2 and a minority carrier band near ckz = π/2. These propagating bands cross the bulk Fermi level at a small in-plane momentum kx,ky (Extended Data Fig. 1h, 1D metal regime in the z direction) but are spread away from the Fermi level at large kx,ky (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| > |γ2| and never crosses the Fermi level (Extended Data Fig. 1h).

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 procedure19, 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 U0 = 8.5 meV, U1 = −17 meV and U3 = −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 |nb| > 1012 cm−2, and all QHE states can be traced back to nb ≈ 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 × 1012 cm−2 resulting in the diamond-like features in σxx occurring at flux fractions ϕ/ϕ0 = 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 |nb| > 2 × 1012 cm−2, independent of nt 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 polynomials21, 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 nh = 2.2 × 1012 cm−2, µh = 24,000 cm2 V−1 s−1, ne = 2.8 × 1012 cm−2 and µe = −19,000 cm2 V−1 s−1 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.

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 ϕ/ϕ0 = p/q. Levels in Fig. 4b were calculated from the tight-binding description of Landau bands in graphite from ref. 26 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 fields46,47. Hofstadter’s butterfly for a honeycomb lattice (Fig. 4a) was calculated from the finite-difference equation in ref. 36, 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 ϕ/ϕ0 = 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 ΔEfractal ≈ 0.1 meV. We assign ΔEfractal to the largest gap in Hofstadter’s butterfly at the flux value ϕ/ϕ0 = 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 ε(ϕ/ϕ0 + ρ) = ε(ϕ/ϕ0) for any integer ρ) to plot states at ϕ/ϕ0 > 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 graphite26 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).

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(U1) using U1 = eVb − e2n/CG (ref. 40). 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 S2+. 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 cm2 V−1 s−1 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 S2+ (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 = (nt − nb)e/2ε0. At D = 0.24 V nm−1, 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 bulk48,49,50.

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 C3 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 captured51,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 breakdown55 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 ELT, 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 × 1012 cm−2. 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.

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 studied56, 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 system57, 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 E2g hBN peak at 1,363 cm−1. The FWHM is 21 cm−1 and 35 cm−1 for non-aligned and aligned regions of the MLG, respectively, which agrees well with the results in ref. 56. 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.

All data are available from the corresponding authors upon reasonable request.

The codes generated in this work to compute surface states of graphite multilayers with and without coupling to aligned hBN are available at https://github.com/slizovskiy/GraphitehBN.

Mak, K. F. & Shan, J. Semiconductor moiré materials. Nat. Nanotechnol. 17, 686–695 (2022).

Article ADS CAS PubMed Google Scholar

Lau, C. N., Bockrath, M. W., Mak, K. F. & Zhang, F. Reproducibility in the fabrication and physics of moire materials. Nature 602, 41–50 (2022).

Article ADS CAS PubMed Google Scholar

Ciarrocchi, A., Tagarelli, F., Avsar, A. & Kis, A. Excitonic devices with van der Waals heterostructures: valleytronics meets twistronics. Nat. Rev. Mater. 7, 449–464 (2022).

Article ADS Google Scholar

Liu, Y. et al. Moiré superlattices and related moiré excitons in twisted van der Waals heterostructures. Chem. Soc. Rev. 50, 6401–6422 (2021).

Article CAS PubMed Google Scholar

Kennes, D. M. et al. Moiré heterostructures as a condensed-matter quantum simulator. Nat. Phys. 17, 155–163 (2021).

Article CAS Google Scholar

He, F. et al. Moiré patterns in 2D materials: a review. ACS Nano 15, 5944–5958 (2021).

Article CAS PubMed Google Scholar

Carr, S., Fang, S. & Kaxiras, E. Electronic-structure methods for twisted moiré layers. Nat. Rev. Mater. 5, 748–763 (2020).

Article ADS CAS Google Scholar

Yankowitz, M., Ma, Q., Jarillo-Herrero, P. & LeRoy, B. J. van der Waals heterostructures combining graphene and hexagonal boron nitride. Nat. Rev. Phys. 1, 112–125 (2019).

Article CAS Google Scholar

Andrei, E. Y. et al. The marvels of moiré materials. Nat. Rev. Mater. 6, 201–206 (2021).

Article ADS CAS Google Scholar

Cao, Y. et al. Correlated insulator behaviour at half-filling in magic-angle graphene superlattices. Nature 556, 80–84 (2018).

Article ADS CAS PubMed Google Scholar

Cao, Y. et al. Unconventional superconductivity in magic-angle graphene superlattices. Nature 556, 43–50 (2018).

Article ADS CAS PubMed Google Scholar

Chen, G. et al. Tunable correlated Chern insulator and ferromagnetism in a moiré superlattice. Nature 579, 56–61 (2020).

Alexeev, E. M. et al. Resonantly hybridized excitons in moiré superlattices in van der Waals heterostructures. Nature 567, 81–86 (2019).

Article ADS CAS PubMed Google Scholar

Regan, E. C. et al. Mott and generalized Wigner crystal states in WSe2/WS2 moiré superlattices. Nature 579, 359–363 (2020).

Article ADS CAS PubMed Google Scholar

Li, H. et al. Imaging two-dimensional generalized Wigner crystals. Nature 597, 650–654 (2021).

Article ADS CAS PubMed Google Scholar

Xu, Y. et al. Correlated insulating states at fractional fillings of moiré superlattices. Nature 587, 214–218 (2020).

Article ADS CAS PubMed Google Scholar

Liu, E. et al. Signatures of moiré trions in WSe2/MoSe2 heterobilayers. Nature 594, 46–50 (2021).

Article ADS CAS PubMed Google Scholar

Zhou, Y. et al. Bilayer Wigner crystals in a transition metal dichalcogenide heterostructure. Nature 595, 48–52 (2021).

Article ADS CAS PubMed Google Scholar

Ponomarenko, L. A. et al. Cloning of Dirac fermions in graphene superlattices. Nature 497, 594–597 (2013).

Article ADS CAS PubMed Google Scholar

Hunt, B. et al. Massive Dirac fermions and Hofstadter butterfly in a van der Waals heterostructure. Science 340, 1427–1430 (2013).

Article ADS CAS PubMed Google Scholar

Krishna Kumar, R. et al. High-temperature quantum oscillations caused by recurring Bloch states in graphene superlattices. Science 357, 181–184 (2017).

Article ADS CAS PubMed Google Scholar

Krishna Kumar, R. et al. High-order fractal states in graphene superlattices. Proc. Natl Acad. Sci. USA 115, 5135–5139 (2018).

Article ADS CAS PubMed PubMed Central Google Scholar

Halbertal, D. et al. Multilayered atomic relaxation in van der Waals heterostructures. Phys. Rev. X. 13, 011026 (2023).

Mandelli, D., Ouyang, W., Urbakh, M. & Hod, O. The princess and the nanoscale pea: long-range penetration of surface distortions into layered materials stacks. ACS Nano 13, 7603–7609 (2019).

Article CAS PubMed Google Scholar

Davison, S. G. & Steslicka, M. Basic Theory of Surface States (Oxford Univ. Press, 1996).

Yin, J. et al. Dimensional reduction, quantum Hall effect and layer parity in graphite films. Nat. Phys. 15, 437–442 (2019).

Article CAS Google Scholar

Shi, Y. et al. Electronic phase separation in multilayer rhombohedral graphite. Nature 584, 210–214 (2020).

Article ADS CAS PubMed Google Scholar

Yang, Y. et al. Stacking order in graphite films controlled by van der Waals technology. Nano Lett. 19, 8526–8532 (2019).

Article ADS CAS PubMed Google Scholar

Dresselhaus, M. S. & Dresselhaus, G. Intercalation compounds of graphite. Adv. Phys. 51, 1–186 (2002).

Article ADS CAS Google Scholar

Morozov, S. V. et al. Two-dimensional electron and hole gases at the surface of graphite. Phys. Rev. B 72, 201401 (2005).

Article ADS Google Scholar

Slonczewski, J. C. & Weiss, P. R. Band structure of graphite. Phys. Rev. 109, 272–279 (1958).

Article ADS CAS Google Scholar

McClure, J. W. Band structure of graphite and de Haas-van Alphen effect. Phys. Rev. 108, 612–618 (1957).

Article ADS CAS Google Scholar

Brown, E. Bloch electrons in a uniform magnetic field. Phys. Rev. 133, A1038–A1044 (1964).

Article ADS MathSciNet Google Scholar

Zak, J. Magnetic translation group. Phys. Rev. 134, A1602–A1606 (1964).

Article ADS MATH Google Scholar

Hofstadter, D. R. Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields. Phys. Rev. 14, 2239–2249 (1976).

Article ADS CAS Google Scholar

Rammal, R. Landau level spectrum of Bloch electrons in a honeycomb lattice. J. Phys. 46, 1345–1354 (1985).

Article MathSciNet Google Scholar

Huber, R. et al. Band conductivity oscillations in a gate-tunable graphene superlattice. Nat. Commun. 13, 2856 (2022).

Article ADS CAS PubMed PubMed Central Google Scholar

Mishchenko, A. et al. Twist-controlled resonant tunnelling in graphene/boron nitride/graphene heterostructures. Nat. Nanotechnol. 9, 808–813 (2014).

Article ADS CAS PubMed Google Scholar

Wang, L. et al. One-dimensional electrical contact to a two-dimensional material. Science 342, 614–617 (2013).

Article ADS CAS PubMed Google Scholar

Yu, G. L. et al. Hierarchy of Hofstadter states and replica quantum Hall ferromagnetism in graphene superlattices. Nat. Phys. 10, 525–529 (2014).

Article CAS Google Scholar

Slizovskiy, S. et al. Out-of-plane dielectric susceptibility of graphene in twistronic and Bernal bilayers. Nano Lett. 21, 6678–6683 (2021).

Article ADS CAS PubMed PubMed Central Google Scholar

Yu, G. L. et al. Interaction phenomena in graphene seen through quantum capacitance. Proc. Natl Acad. Sci. USA 110, 3282–3286 (2013).

Article ADS CAS PubMed PubMed Central Google Scholar

Luryi, S. Quantum capacitance devices. Appl. Phys. Lett. 52, 501–503 (1988).

Article ADS Google Scholar

Wallbank, J. R. et al. Excess resistivity in graphene superlattices caused by umklapp electron–electron scattering. Nat. Phys. 15, 32–36 (2019).

Article CAS Google Scholar

Lee, M. et al. Ballistic miniband conduction in a graphene superlattice. Science 353, 1526–1529 (2016).

Article ADS CAS PubMed Google Scholar

Sugihara, K. Charge-density wave and magnetoresistance anomaly in graphite. Phys. Rev. B 29, 6722–6731 (1984).

Article ADS CAS Google Scholar

Arnold, F. et al. Charge density waves in graphite: towards the magnetic ultraquantum limit. Phys. Rev. Lett. 119, 136601 (2017).

Article ADS CAS PubMed Google Scholar

Eisenstein, J. P., Pfeiffer, L. N. & West, K. W. Negative compressibility of interacting two-dimensional electron and quasiparticle gases. Phys. Rev. Lett. 68, 674–677 (1992).

Article ADS CAS PubMed Google Scholar

Eisenstein, J. P., Pfeiffer, L. N. & West, K. W. Compressibility of the two-dimensional electron gas: measurements of the zero-field exchange energy and fractional quantum Hall gap. Phys. Rev. B 50, 1760–1778 (1994).

Article ADS CAS Google Scholar

Li, L. et al. Very large capacitance enhancement in a two-dimensional electron system. Science 332, 825–828 (2011).

Article ADS CAS PubMed Google Scholar

Wilkinson, M. Critical properties of electron eigenstates in incommensurate systems. Proc. R. Soc. Lond. A 391, 305–350 (1997).

ADS MathSciNet Google Scholar

Cohen, M. H. & Falicov, L. M. Magnetic breakdown in crystals. Phys. Rev. Lett. 7, 231–233 (1961).

Article ADS CAS Google Scholar

Davis, L. C. & Liu, S. H. Landau spectrum and line broadening in real metals. Phys. Rev. 158, 689–697 (1967).

Article ADS CAS Google Scholar

Berry, M. V. & Mount, K. E. Semiclassical approximations in wave mechanics. Rep. Prog. Phys. 35, 315 (1972).

Article ADS CAS Google Scholar

Alexandradinata, A. & Glazman, L. Semiclassical theory of Landau levels and magnetic breakdown in topological metals. Phys. Rev. B 97, 144422 (2018).

Article ADS CAS Google Scholar

Eckmann, A. et al. Raman fingerprint of aligned graphene/h-BN superlattices. Nano Lett. 13, 5242–5246 (2013).

Article ADS CAS PubMed Google Scholar

Cheng, B. et al. Raman spectroscopy measurement of bilayer graphene’s twist angle to boron nitride. Appl. Phys. Lett. 107, 033101 (2015).

Article ADS Google Scholar

Download references

This research was supported by the European Research Council (ERC) under the Horizon 2020 research and innovation programme of the European Union (grant agreement no. 865590), the Royal Society International Exchanges 2019 Cost Share Program (IEC\R2\192001), the Research Council UK (BB/X003736/1), EC-FET Core 3 European Graphene Flagship Project, Quantum Flagship Project 2D-SIPC, EPSRC grants EP/S030719/1 and EP/V007033/1, and the Lloyd Register Foundation Nanotechnology Grant. J.Y., S.X. and Y.Y. acknowledge the support of the National Natural Science Foundation of China (12150002 to J.Y., 12274354 to S.X., 12104505 to Y.Y.). Q.Y. acknowledges the funding from Leverhulme Early Career Fellowship (ECF-2019-612), the Royal Society University Research Fellowship (URF\R1\221096) and the Research Council UK (EP/X017575/1). K.S.N. is grateful to the Ministry of Education, Singapore, for the support under its Research Centre of Excellence award to the Institute for Functional Intelligent Materials (I-FIM, project no. EDUNC-33-18-279-V12) and to the Royal Society (UK, grant no. RSRP\R\190000). C.M. acknowledges the funding from Graphene NOWNANO CDT EP/L01548X/1.

These authors contributed equally: Ciaran Mullan, Sergey Slizovskiy, Jun Yin

Department of Physics and Astronomy, University of Manchester, Manchester, UK

Ciaran Mullan, Sergey Slizovskiy, Jun Yin, Ziwei Wang, Qian Yang, Yaping Yang, Kostya S. Novoselov, A. K. Geim, Vladimir I. Fal’ko & Artem Mishchenko

National Graphene Institute, University of Manchester, Manchester, UK

Sergey Slizovskiy, Qian Yang, Shuigang Xu, Yaping Yang, Sheng Hu, Kostya S. Novoselov, A. K. Geim, Vladimir I. Fal’ko & Artem Mishchenko

State Key Laboratory of Mechanics and Control for Aerospace Structures, Key Laboratory for Intelligent Nano Materials and Devices of Ministry of Education, Institute for Frontier Science, Nanjing University of Aeronautics and Astronautics, Nanjing, China

Jun Yin

Key Laboratory for Quantum Materials of Zhejiang Province, Department of Physics, School of Science, Westlake University, Hangzhou, China

Shuigang Xu

Laboratoire National des Champs Magnétiques Intenses (LNCMI), CNRS Université Grenoble Alpes, Université Toulouse 3, INSA Toulouse, EMFL, Grenoble, France

Benjamin A. Piot

State Key Laboratory of Physical Chemistry of Solid Surfaces, Collaborative Innovation Center of Chemistry for Energy Materials (iChEM), College of Chemistry and Chemical Engineering, Xiamen University, Xiamen, China

Sheng Hu

National Institute for Materials Science, Tsukuba, Japan

Takashi Taniguchi & Kenji Watanabe

Institute for Functional Intelligent Materials, National University of Singapore, Singapore, Singapore

Kostya S. Novoselov

Henry Royce Institute for Advanced Materials, Manchester, UK

Vladimir I. Fal’ko

You can also search for this author in PubMed Google Scholar

J.Y., S.X., Y.Y. and S.H. designed and fabricated the devices. C.M., J.Y. and B.A.P. performed transport measurements. C.M., J.Y., S.S. and A.M. analysed the data. S.S., C.M. and V.I.F. developed the theory and performed the calculations. Z.W. performed Raman measurements and analysis. T.T. and K.W. provided hBN crystals. A.M., A.K.G., C.M., S.S., J.Y., Z.W., K.S.N. and Q.Y. contributed to writing the paper. All authors discussed the results and commented on the paper.

Correspondence to Jun Yin, Vladimir I. Fal’ko or Artem Mishchenko.

The authors declare no competing interests.

Nature thanks the anonymous reviewers for their contribution to the peer review of this work.

is available for this paper at https://doi.org/10.1038/s41586-023-06264-5.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

a, Capacitance (C) vs carrier density (n) for device D6 at T = 0.3 K. Inset shows optical micrograph of a typical capacitor device, scale bar 20 µm. b, Density of states (DoS) and quantum capacitance (Cq) vs surface potential (U1) from calculations based on effective-mass model (black line). Coloured symbols are experimental data from four different devices (D5, D6, D8, D11). Inset is calculated U1 vs n. c–f, Calculated dispersions of a 20-layer-thick graphite film for hole/electron dopings n = −6, −4, 4, 6 × 1012 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 wavefunction 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. g, same as c–f, calculated for zero doping where no distinct surface states are observed. h,i, 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 h, 1D metal regime with Type 2 evanescent modes) and 0.15 eV (panel i, 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.

a, c, σxx,σxy as a function of B and nb for aligned device D1 at T = 0.24 K. Landau level features emanate from LTs (|n| ≈ 2.0 and 3.7 × 1012 cm−2) and from zero doping. b, d, σxx,σxy as a function of B and nb for non-aligned device D4 at T = 0.22 K. e, Comparison of low field σxx line traces for the two interfaces of the single aligned D1, where aligned interface (nb) hosts many features related to the LTs which are notably absent in non-aligned interface (nt) response. f, Mapping of σxx for D1 as a function of nb and nt, highlighting that the vertical features related to LTs are independent of nt throughout the measured range. B = 0.3 T, T = 1.8 K, and the colour scale is black to white, 5 to 18 mS.

a, Extended range plot of Fig. 2b in the main text for device D1; ∆σxx (conductivity minus a smooth background) as a function of B and nb at T = 60 K. Inset, σxx(ϕ0/ϕ) trace at nb = 3.53 × 1012 cm−2 highlighting higher order Brown-Zak oscillation at ϕ0/ϕ = 5/2. b, Extended range plot of Fig. 2c in the main text; σxx as a function of B and nb at T = 20 K. Inset, conductivity averaged across a carrier density range nb = 2.00 to 3.84 × 1012 cm−2, ‹σxx(nb)›, plotted as a function of ϕ0/ϕ showing oscillations continuing down to experimental mapping resolution (Bstep = 1 mT). c–h High temperature mappings of σxx(nt, B) (panels c & g) and σxx(nb, B) (panels d & h) for doubly aligned device D3 (c–d) and singly aligned device D2 (g–h). Measurements were conducted at T = 60 K, colour scale black to white is 80 to 190 µS for c & d, 49 to 154 µS for g and 64 to 97 µS for h. Panels e & f show 1st and 2nd derivatives, respectively, calculated from d. Horizontal dashed lines show the significant flux fractions \(\phi /{\phi }_{0}\) from 1/2 to 1/8 fitted to the data. Colour scale for e is blue to red −40 to 40 µS T−1 and in f is black to white −100 to 100 µS T−2.

a, Dispersion in the SBZ plotted up to Fermi energy for doping 2.1 × 1012 cm−2 (left) and 3.8 × 1012 cm−2 (right) respectively. b, networks of classical trajectories accessible to surface state electrons at the same dopings in a. Examples of shortest interfering paths of equal length are shown by green and brown arrows. The area enclosed between these paths equals ABZ/(eB)2 irrespective of energy. c, Conductivity enhancement at high B field in aligned Corbino devices; σxx as a function of single gate induced carrier density for non-aligned device (D10, black curve, bottom gate tuned), singly aligned device (D2, red curve, top gate tuned) and doubly aligned device (D3, blue curve, bottom gate tuned) all in Corbino geometry, measured at T = 0.3 K and B = 18 T. Shaded regions highlight surface Landau band features. d, Conductivity at zero B field in aligned (D1) vs non-aligned (D9) Hall bar devices. Increased scattering due to alignment results in reduced conductivity in the bulk at B = 0 T. Here the top (non-aligned interface) gate is swept for D1. T = 0.3 K. e, Brown-Zak oscillations in σxx differ in amplitude non-monotonically as a function of ϕ0/ϕ in aligned Hall bar device D1. Measured at nb = 2.6 × 1012 cm−2, where T = 60 K curve is from the same dataset used to generate Fig. 2b and Extended Data Fig. 2a.

a, σxx(nt, nb) of device D2 as a function of top and bottom gates, nb and nt, measured at high field B = 15.6 T, T = 0.3 K, colour scale: brown to yellow, 0 to 59 μS. b, Wannier diagram depicting the QHE and fractal QHE gaps in the bulk as diagonal grey and purple lines, respectively. High doping for top (bottom) gate results in horizontal (vertical) features that correspond to surface states accessing the +2 and −2 Landau bands from the opposite graphite surfaces (highlighted by orange shading). c, σxx(nt, nb) of device D3, measured at high field B = 13.6 T, T = 0.3 K, colour scale: brown to yellow, 0 to 98 μS. d, The same as b but for D3.

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. 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 ν = 0 gap. e, Conductivity map for device D2, same data as in Fig. 3a in the main text, except plotted as a function of filling factor ν of main QHE sequence. Colour scale: brown to yellow, 0 to 80 μ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é perturbation26. 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. 4a in the main text by applying Hofstadter’s butterfly as a small perturbation (S = 0.42 meV) to each energy level, Eq. 2. h, Same as g but with S = 0.42 meV for odd layer states and S = 0.12 meV for even layer states.

a, Typical C(n) at 0.6 T (bottom black curve), 5.1 T (middle blue curve) and 12 T (top red curve), respectively, in capacitor device D6 at T = 0.3 K. Bottom insets are magnifications of the marked areas. b, DoS vs U1 at 5.1 T, 0.3 K. c, Surface Landau fan diagram C (n, B) at 0.3 K. Colour scale: navy to white, 254.7 fF to 270.0 fF. Dashed line marks the critical field B = 7.5 T, where bulk Landau band 2+ leaves the Fermi level and graphite enters the UQR. d, Dispersion relation for bulk Landau bands calculated using the SWMC model at B = 7.5 T26.

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.6 T, 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 Extended Data Fig. 7a, based on the periodicity of the oscillations.

a,b, Fractional surface states measured in capacitance (device D8). a, C(n) at B = 20 T, T = 0.3 K. Inset magnifies the encircled region, but plots as C(ν). b, dC/dν (ν, B) of S2+ in high B region. Colour scale: blue to red, −6 to 6 fF. Right-top inset shows the corresponding C (n, B) map. Colour scale: orange to red, 254.3 to 260 fF. c, Longitudinal conductivity σxx (B, n) measured at D = 0.24 V/nm in a Hall bar device D9, T = 0.3 K, I = 20 nA. The 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, 0.1 to 118.2 µS. d, σxx cut profile (black curve) of the white dashed line in c and the corresponding Hall conductivity σxy (red curve).

a, Optical image and b, AFM profile of graphite flake used in Raman measurement (area shown in the black box in a). The flake contains both regions of MLG and graphite with thickness around 10 layers. Scale bar 10 um. c, Optical image of the stack fabricated for Raman measurements. The aligned and non-aligned regions are marked by blue and red dashed lines. Scale bar 10 um. d, Raman map of full width half maximum (FWHM) of 2D peak. The MLG and graphite regions are colour-coded in grayscale and red-to-purple, respectively. e, Comparison of 2D peak between aligned and non-aligned regions in MLG (up) and graphite (down). Each spectrum shown here is averaged over ten spectra at different spots. The intensity is normalized by that of hBN peak at 1363 cm−1.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

Mullan, C., Slizovskiy, S., Yin, J. et al. Mixing of moiré-surface and bulk states in graphite. Nature (2023). https://doi.org/10.1038/s41586-023-06264-5

Download citation

Received: 04 December 2022

Accepted: 25 May 2023

Published: 19 July 2023

DOI: https://doi.org/10.1038/s41586-023-06264-5

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.

Lithium Energy advances metallurgical testwork and PFS on downstream graphite processing for battery industry

Gemcorp Supports Walkabout’s Tanzanian Graphite Mine

News

Mixing of moiré