Films of rhombohedral graphite as two-dimensional topological semimetals

TL;DR

This paper explores thin films of rhombohedral graphite as two-dimensional topological semimetals with tunable Berry curvature and magnetic properties, revealing potential for novel electronic transport phenomena.

Contribution

It demonstrates that rhombohedral graphite films maintain 2D topological states with large Berry curvature and magnetic moments, controllable via electrostatic gating.

Findings

Large Berry curvature and magnetic moments in rhombohedral graphite films.

Electrostatic gating tunes the anomalous Hall effect.

Unique Landau level structure observed in these films.

Abstract

Topologically non-trivial states characterized by Berry curvature appear in a number of materials ranging from spin-orbit-coupling driven topological insulators to graphene. In multivalley conductors, such as mono- and bilayer graphene, despite a zero total Chern number for the entire Brillouin zone, Berry curvature with different signs concentrated in different valleys can affect the observable material's transport characteristics. Here we consider thin films of rhombohedral graphite, which appear to retain truly two-dimensional properties up to tens of layers of thickness and host two-dimensional electron states with a large Berry curvature, accompanied by a giant intrinsic magnetic moment carried by electrons. The size of Berry curvature and magnetization in the vicinity of each valley can be controlled by electrostatic gating leading to a tuneable anomalous Hall effect and a peculiar structure of the two-dimensional Landau level spectrum.

Equations

$$\hat{H}$$

$$\tilde{\Delta }$$

$${\Omega }^{(\pm )}(\mathbf{p})$$

$${\mathfrak{m}}_z(\mathbf{p})$$

$${\Omega }^{(\pm )}(\mathbf{p})\approx \pm \xi \frac{2N^2}{p^2{\gamma }_1}\frac{\tilde{\Delta }(\mathbf{p})(p/p_{\mathrm{c}}{)}^{2N}}{{\left[(\tilde{\Delta }(\mathbf{p})/{\gamma }_1{)}^2+4(p/p_{\mathrm{c}}{)}^{2N}\right]}^{3/2}},$$

$${\mathfrak{m}}_z(\mathbf{p})\approx \xi \frac{e}{\hslash }\frac{N^2}{p^2}\frac{\tilde{\Delta }(\mathbf{p})(p/p_{\mathrm{c}}{)}^{2N}}{(\tilde{\Delta }(\mathbf{p})/{\gamma }_1{)}^2+4(p/p_{\mathrm{c}}{)}^{2N}}.$$

$${\mathbf{v}}_{\pm }(\mathbf{p})={\nabla }_{\mathbf{p}}{ϵ}_{\pm }(\mathbf{p})+(e/\hslash )\mathbf{E}\times {\hat{\mathbf{e}}}_z{\Omega }^{(\pm )}(\mathbf{p}),$$

$${\sigma }_{xy}^{\mathrm{A}}=-\frac{e^2{B}_z}{{\pi }^2{\hslash }^3}\sum _i{\oint }_{{\mathcal{L}}_i(K)}\frac{{\Omega }_i^{(\pm )}{\mathfrak{m}}_{i,z}}{|{\nabla }_{\mathbf{p}}{ϵ}_i|}dp.$$

$${\sigma }_{xy}^{\mathrm{H}}=-\frac{e^3{\tau }^2}{{\pi }^2{\hslash }^2}\sum _{i,\gamma }{\oint }_{{\mathcal{L}}_i(K)}\frac{({\partial }_{p_x}{ϵ}_i)}{|{\nabla }_{\mathbf{p}}{ϵ}_i|}[{\nabla }_{\mathbf{p}}{ϵ}_i\times \mathbf{B}{]}_{\gamma }\frac{dp}{m_{\gamma y}},$$

$$\delta {\sigma }_{xy}$$

$${\hat{H}}_N=\left(\begin{array}{ccccc}{D}_1& V_{1,2}& W& 0& \cdots \\ V_{1,2}^{†}& D_2& V_{2,3}& W& \cdots \\ W^{†}& V_{2,3}^{†}& D_3& V_{3,4}& \cdots \\ 0& W^{†}& V_{3,4}^{†}& D_4& \cdots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right),$$

$$D_1$$

$$D_n$$

$$V_{n,n+1}$$

$${\kappa }_{n,n}$$

$${\kappa }_{n,n+1}$$

$${\hat{H}}_N={\gamma }_1\left(\begin{array}{cccccc}0& 1& 0& 0& 0& \cdots \\ 1& 0& 1& 0& 0& \cdots \\ 0& 1& 0& 1& 0& \cdots \\ 0& 0& 1& 0& 1& \cdots \\ 0& 0& 0& 1& 0& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right),$$

$${\Delta }_{\mathrm{bulk}}\approx 3\pi {\gamma }_1/N,$$

$$S(C)=\frac{2\pi }{{\lambda }_B^2}\left(n+\frac{1}{2}-\frac{\Gamma (C)}{2\pi }\right),$$

$$\dot{\mathbf{p}}$$

$$\dot{\mathbf{r}}$$

Full text

Introduction

Berry curvature berry84 ; xiao10 is a measure of the topological nature of non-trivial states appearing in materials ranging from spin-orbit-coupling driven topological insulators konig07 ; hsieh08 ; hasan10 ; asboth16 to graphene novo05 ; zhang05 ; novo06 . In multivalley conductors, such as graphene, Berry curvature with different signs concentrated in different valleys xiao07 can affect the material’s observable properties even though the Chern number for the entire Brillouin zone may be zero. Recently, there has been renewed interest in rhombohedral graphite mcclure69 due to progress in fabricating and characterizing thin films pierucci15 ; henni16 ; henck18 ; myhro18 ; lat19 ; geisenhof19 . Rhombohedral is one of the structural phases of graphite which has a specific ‘ABC’ stacking of consecutive honeycomb layers of carbon atoms such that every atom has a nearest neighbor from an adjacent layer either directly above or underneath it. The interlayer hybridization of all carbon $P_{z}$ orbitals with characteristic energy $\gamma_{1}\approx 0.38$ eV leads to a gapped electronic spectrum in the middle of a thin film of $N$ layers (the bulk gap $\approx 3\pi\gamma_{1}/N$, see Methods). However, in the surface layers of the film, half of the carbon atoms don’t have a nearest neighbor in the next layer for hybridizing their $P_{z}$ orbitals, leading, as in the Su-Schrieffer-Heeger model su79 ; asboth16 , to low-energy surface states. These low-energy states form bands that cover the entire energy range within the thin film bulk gap guinea06 ; manes07 ; min08 ; koshino09 ; heik11 ; xiao11 ; kopnin13 ; ho16 .

In this article, we theoretically model thin films of rhombohedral graphite and find that they retain their two-dimensional nature for tens of layers of thickness. The low-energy surface states give rise to a semi-metallic band structure with two bands that are almost degenerate near the valley center, but split apart at $p\sim p_{\mathrm{c}}$ where the dispersion is highly anisotropic ($p_{\mathrm{c}}=\gamma_{1}/v$ where $v$ is the Dirac velocity of electrons in graphene determined by the intralayer carbon-carbon hopping parameter $\gamma_{0}$). The presence of spatial asymmetry between the surface layers, which may be controlled using an electric displacement field applied perpendicularly to the film mcc06 ; mcc06b or which may arise from interaction-induced spontaneous symmetry breaking min08b ; jung11b , can create an insulator with an energy band gap and topologically non-trivial states represented by a giant Berry curvature and intrinsic magnetic moment of electrons. We predict that the topological nature of the surface states will be manifest in electronic transport properties including a large anomalous Hall effect and anomalous transverse photoconductivity. In addition, these features are reflected in the electronic spectra in the presence of a perpendicular magnetic field (the Landau level spectra) whereupon spatial asymmetry breaks valley degeneracy with different patterns of level crossing and hybridization in the two valleys.

Results

Semi-metallic band structure. When studied taking into account all details of intra- and interlayer carbon-carbon couplings (see Methods) in the full Slonzewski-Weiss-McClure (SWMcC) tight-binding model mcclure69 ; dressel02 ; koshino09 ; kopnin13 , the surface states in a thin film of rhombohedral graphite produce a semi-metallic band structure illustrated in Fig. 1a. The dispersion in Fig. 1a is plotted as a function of the in-plane momentum $\mathbf{p}$ counted in a zig-zag direction from the center of the valley $K$ and normalized by $p_{\mathrm{c}}=\gamma_{1}/v$. In contrast to some earlier studies ho16 ; henni16 , the dispersion in Fig. 1 takes into account both skew interlayer couplings $\gamma_{3}$ and $\gamma_{4}$. The two bands of the surface states in an ABC graphite film, below referred to as conduction (blue) and valence (red), are almost degenerate near the valley center, splitting apart in the momentum range $p\gtrsim p_{\mathrm{c}}$. The electron dispersion at $p\sim p_{\mathrm{c}}$ is highly anisotropic, due to trigonal warping effects generated by skew interlayer hoppings, with inverted orientation in valley $K^{\prime}$.

The interplay between these factors makes an undoped film of rhombohedral graphite a two-dimensional semi-metal (2DSM) as its Fermi level lies within both the conduction and valence surface bands. It becomes a two-dimensional metal (2DM) upon n- or p-doping when the Fermi level lies in only one of the conduction or valence surface bands and, eventually, a bulk metal (3DM) where the Fermi level lies within the bulk bands. In Fig. 1b, we identify parametric regimes for each of these three cases, taking into account the dependence of the spectrum on the number of layers. The parametric diagram in Fig. 1b was built by brute-force diagonalization of a hybrid $\mathbf{k}\cdot\mathbf{p}$ tight-binding approach model (HkpTB) based on the full SWMcC model, in which the intralayer hopping of electrons between carbon atoms is taken into account in a continuous description of sublattice Bloch states using $\mathbf{k}\cdot\mathbf{p}$ theory in the $K$ and $K^{\prime}$ valleys, combined with interlayer hoppings introduced in the spirit of a tight-binding model (see Methods). For a film with $N$ layers, this involves finding eigenvalues and eigenstates of a $2N\times 2N$ matrix acting in the space of the sublattice Bloch states in each valley. When studied in the presence of a vertical electric displacement field (perpendicular to the thin film), the bands responsible for the semimetallicity split by $\Delta=U_{1}-U_{N}$, where $U_{1}$ and $U_{N}$ are the onsite energies of the surface layers, so that the system may be tuned into a gap-full insulator, Fig. 1c. It may also be possible to induce a spectral gap by superconductivity kopnin13 or spontaneous symmetry breaking into a magnetic state zhang11 ; jia13 ; pamuk17 , making $\Delta$ spin and/or valley dependent as in the layer antiferromagnetic configurations discussed in the context of bilayer graphene min08b ; jung11b and experimentally observed in rhombohedral graphite for $N=3$ bao11 ; lee14 and $N=4$ myhro18 .

When subjected to an external magnetic field perpendicular to the plane of a film, the low-energy spectrum splits into interweaving electron-like and hole-like Landau levels (LLs). For a film of rhombohedral graphite, a representative example is shown in Fig. 1d computed for $N=10$ layers (see Methods), with electron-like LLs dispersing upwards for energy $E\gtrsim 20\,$meV, hole-like LLs dispersing downward for $E\lesssim 0\,$meV, and both kinds of levels present in the semimetallic range $0\lesssim E\lesssim 20\,$meV, where one can relate the electron-like and hole-like LLs to different parts of the 2D electron Fermi lines illustrated in the insets. In Fig. 1e, we show how the LL spectrum is modified by the opening of a gap (electrostatically controlled, or induced by exchange energy) which enhances some of the avoided crossings in the LL spectrum.

Two-band model. To develop intuition about the LL spectra shown in Fig. 1 as well as to anticipate the magnetoconductivity of an ABC film, we use a simplified two-band model which we derived from the full HkpTB by projecting the $2N\times 2N$ Hamiltonian onto the pair of surface states in the outermost (bottom and top) layers,

[TABLE]

Here, $m_{\star}\sim 0.4m_{0}$ (see Methods), we use a vector ${\boldsymbol{\sigma}}=(\sigma_{x},\sigma_{y},\sigma_{z})$ of Pauli matrices acting in the space of surface states and ${\bf d}=(d_{x},d_{y},d_{z})$, and incorporate the arbitrarily-oriented magnetic field with components both perpendicular, $B_{z}$, and parallel, $\mathbf{B}_{\parallel}=(B_{x},B_{y})$, to the film. In Eq. (3), the two valleys of graphite ${\bf K}=(4\pi/3a,0)$ and ${\bf K^{\prime}}=-{\bf K}$ are indexed by $\xi=1$ and $-1$ respectively. In the expression for $\tilde{\Delta}$ we take into account the effect of the Lorentz boost experienced by electrons tunneling between surface states, where $\beta_{1}=\frac{2\gamma_{4}}{\gamma_{0}}(2N-3)+\frac{2\delta}{\gamma_{1}}(N-1)$ ($\beta_{1}\propto N$ for $N\gg 1$), $a_{z}$ is the interlayer distance and $-e$ is the electronic charge. Parameter $\delta$ represents a difference in energy between the dimer sites inside the crystal and non-dimer sites $A_{1}$ and $B_{N}$ in the outer layers. In the definition of operator $\kappa$ and its Hermitian conjugate, $\kappa^{\dagger}$, we use the Landau gauge for the out-of-plane magnetic field $B_{z}$, whereas, for $B_{z}=0$, $\kappa=\xi p_{x}+ip_{y}$. The off-diagonal element, $X(\mathbf{p})=\sum\frac{(n_{1}+n_{2}+n_{3})!}{n_{1}!n_{2}!n_{3}!}\left(-\frac{\kappa^{\dagger}}{p_{\mathrm{c}}}\right)^{n_{1}}\left(-\frac{\gamma_{3}}{\gamma_{0}}\frac{\kappa}{p_{\mathrm{c}}}\right)^{n_{2}}\left(-\frac{\gamma_{2}}{2\gamma_{1}}\right)^{n_{3}}$, arises from interlayer hoppings (skew hopping between neighboring layers $\gamma_{3}$ and ‘vertical’ next-layer hopping $\gamma_{2}$ in addition to $\gamma_{1}$ mentioned earlier) passing electrons from one surface to another koshino09 , and the summation is taken over integers $n_{i}\geq 0$ such that $n_{1}+2n_{2}+3n_{3}=N$.

Qualitatively, the spectrum of $\hat{H}$, $\epsilon_{\pm}({\bf p})=p^{2}/(2m_{\star})\pm|{\bf d}|$, reproduces the exact multilayer $2N\times 2N$ solutions shown in Fig. 1, and it is particularly useful to discuss the features of the LL spectra (for a detailed comparison, see Methods). Without symmetry breaking, $\tilde{\Delta}=0$, the spectrum of $\hat{H}$ is valley degenerate, leading to a $4$-fold degeneracy (spin and valley) of each of the LLs in Fig. 1d. In the low magnetic field range, one can also identify closely-packed groups of 3 LLs (i.e. 12-fold degenerate) whose origin we trace to three mini-valleys forming at $p\sim p_{\mathrm{c}}$ sketched in Fig. 1b. At high magnetic fields, $\hat{H}$ generates $N$ separate groups of 4-fold degenerate low-energy LLs, which were considered to be degenerate in previous studies min08 ; henni16 where the effects of $\gamma_{4}$ and $\delta$ were neglected, and $X(\mathbf{p})\approx(-\kappa^{\dagger}/p_{\mathrm{c}})^{N}$ resulted in a Berry phase $N\pi$ singularity at $p=0$. At high fields, this group of $N$ LLs clearly separates from electron- and hole- like levels in the spectrum, whereas, at low fields, their mixing with hole-like dispersive LLs leads to an additional two-fold degeneracy for all but the $\epsilon=0$ level.

Topological properties. The other notable feature of the low-energy ‘non-dispersive’ LLs is that their states in opposite valleys are located on opposite surfaces of the film. Consequently, when asymmetry, $\Delta$, between the bottom and top surfaces is introduced by, e.g, a displacement field, $E_{z}$, these LLs split apart, lifting the valley degeneracy. This evolution can be traced in the LL spectrum shown in Fig. 1e where $K$ and $K^{\prime}$ valley states are marked using different colors. The first noticeable difference is that, now, groups of $N$ lowest-energy LLs in opposite valleys can be traced to $\pm\Delta/2$ convergence points marked by arrows in Fig. 1e. To highlight the difference in the LLs spectra in the opposite valleys, we plot them separately in Fig. 2a, b and point out that groups of LLs that converge towards the top of the valence band (highlighted in black) slightly differ in valleys $K$ and $K^{\prime}$. As these states originate from the valence band maxima, where the electron dispersion is approximately parabolic, the difference between LL spectra reflects topological properties of electron states as characterized by Berry curvature $\Omega^{(\pm)}(\mathbf{p})$ and associated magnetic moment $\mathfrak{m}_{z}(\mathbf{p})$ berry84 ; xiao10 ,

[TABLE]

where u_{\pm}^{T}={\cal N}_{\pm}\left(\begin{array}[]{cc}d_{z}\pm d,&d_{x}+id_{y}\\ \end{array}\right) is the sublattice Bloch spinor in the $\epsilon_{\pm}$ (conduction/valence) band, ${\cal N}_{\pm}=1/\sqrt{2d(d\pm d_{z})}$, $d=|{\bf d}|$, and $\nabla_{\mathbf{p}}=(\partial_{p_{x}},\partial_{p_{y}})$. Note that $\Omega^{(-)}(\mathbf{p})=-\Omega^{(+)}(\mathbf{p})$, whereas $\mathfrak{m}_{z}(\mathbf{p})$ is the same for the $\epsilon_{+}$ and $\epsilon_{-}$ bands.

Using the two-component model (3) for $N\gg 1$ and neglecting trigonal warping (parameters $\gamma_{2}$ and $\gamma_{3}$), the Berry curvature $\Omega^{(\pm)}$ for conduction/valence bands is given by zhang11

[TABLE]

which generalizes the result for finite $\Delta$ determined in monolayer and bilayer graphene xiao07 ; xiao10 . The orbital magnetic moment for both bands is

[TABLE]

We note the strong $N^{2}$ layer dependence and the fact that in-plane magnetic field $\mathbf{B}_{\parallel}$ can contribute to the gap $\tilde{\Delta}(\mathbf{p})$.

Whereas $\Omega^{(\pm)}(\mathbf{p})$ and $\mathfrak{m}_{z}(\mathbf{p})$ are peaked at $p\approx 0$ in monolayer and bilayer, for rhombohedral graphene with $N\gg 1$ they are peaked at $p\approx p_{\mathrm{c}}$. Moreover, trigonal warping (parameters $\gamma_{2}$ and $\gamma_{3}$) introduces anisotropy such that $\Omega^{(\pm)}(\mathbf{p})$ and $\mathfrak{m}_{z}(\mathbf{p})$ are peaked in the vicinity of the valence band maxima. This is demonstrated in Fig. 2c where we plot the distribution of magnetic moment $\mathfrak{m}_{z}$ across momentum space in both valleys for $\Delta=10\,$meV. This magnetic moment leads to a valley splitting, $2\mathfrak{m}_{z}B_{z}$, (e.g. of the LL converging towards the valence band edge, Fig. 2a), that can be characterized by a factor $g_{\mathrm{v}}=\mathfrak{m}_{z}/\mu_{\mathrm{B}}$ which can be as large as $g_{\mathrm{v}}\sim 100$, Fig. 3a. An application of in-plane magnetic field would additionally lift the degeneracy between the three valence band maxima, as illustrated using the LL fan in Fig. 2b.

Anomalous Hall effect. In the small magnetic field range where electron dynamics can be described classically (rather than using Landau quantization), the splitting of magnetic moments at the valence band edges in the opposite valleys would shrink the size of the Fermi pockets for the holes in one valley and expand them in the other. The resulting valley polarization of holes would manifest itself in transport characteristics. This is because, in the presence of in-plane electric field ${\bf E}$, carriers in bands with Berry curvature experience a drift  xiao10 ; nagaosa10 ,

[TABLE]

so that, together, Berry curvature and valley polarization produce an anomalous contribution to the Hall effect,

[TABLE]

Here $i$ lists all Fermi lines ${\cal L}_{i}$ in valley $K$ for both $\epsilon_{\pm}$ bands, and the linear integral is taken in the anticlockwise direction along each Fermi line, ${\cal L}_{i}$. When combined with the classical kinetic Hall coefficient,

[TABLE]

computed for the same band structure ($m_{\alpha\beta}^{-1}=\tfrac{\partial^{2}\epsilon_{i}}{\partial p_{\alpha}\partial p_{\beta}}$) and elastic scattering time $\tau=5\times 10^{-14}\,$s, this yields the overall Hall conductivity, $\sigma_{xy}$, displayed in Fig. 3b. The anomalous contribution is most pronounced for small scattering times, $\tau\sim 10^{-13}\,$s, otherwise the classical kinetic Hall conductivity dominates, as shown in Fig. 3c, where we plot the ratio $\sigma_{xy}/(B\tau^{2})$ for $\tau=10^{-13},2\times 10^{-13},10^{-12}\,\rm s$ [note that the value obtained for $\tau=10^{-12}\,\rm s$ is indistinguishable from $\sigma^{\mathrm{H}}_{xy}/(B\tau^{2})$].

The doping density dependence of the overall Hall conductivity, shown in Fig. 3b for ten- and twenty-layer films with $\Delta=10$ meV and $40$ meV, carries certain features reflecting the distribution of Berry curvature and magnetic moment (Fig. 2c) across the electron dispersion. According to the parametric plot shown on the basal plane of Fig. 3a, both ten- and twenty-layer films with $\Delta=40$ meV are gapful 2D semiconductors (2DSC), so that the states with the maximium $\mathfrak{m}_{z}\Omega$ near the valence band edge (Fig. 1c) are reached at a relatively small p-doping, resulting in a hump in $\sigma_{xy}(n_{\mathrm{e}})$ at $n_{\mathrm{e}}<0$ indicated by a star. For a smaller gap, the films appear to be 2D semimetals (2DSM), so that the part of the band where $\mathfrak{m}_{z}\Omega$ is the largest is set above the Fermi level for undoped materials, shifting the hump in $\sigma_{xy}$ to $n_{\mathrm{e}}>0$. A double-peak structure is caused by the convolution of the ${\mathfrak{m}}_{z}\Omega$ product and the density of states with a pronounced Van Hove singularity (pointed by down arrows) in Eq.(6). At higher electron dopings, the Fermi level reaches the maximum of ${\mathfrak{m}}_{z}\Omega$ for the conduction band, causing a negative peak in Hall coefficient (pointed by up arrows).

Pumping inter-band transitions with circularly polarized photons induces a partial valley polarization in graphene and graphite (thus breaking the time-inversion symmetry). Combined with the Berry curvature effect, this would produce a Hall-like drift current (at $B_{z}=0$) perpendicular to the static in-plane electric field, which can be characterized by an anomalous transverse photoconductivity, $\delta\sigma_{xy}$, which appears to be especially pronounced for ABC graphitic films in the 2DSC regime. This effect is determined by the resonant inter-band absorption of circularly-polarized photons with energy $\omega$ and absorption coefficient $g(\omega)$, plotted in Fig. 3d,

[TABLE]

where $\zeta=\pm 1$ stands for left/right handed circular polarization of the pump (radiation is approaching from the top) with power density $W$, $\Omega^{-}_{K}$ is the Berry curvature of the valence band at valley $K$, and $\tau_{\mathrm{rec}}$ is the lifetime of photo-excited holes at the top, $p_{\mathrm{max}}$, of the valence band in the photo-activated valley.

Discussion

We have shown that thin films of rhombohedral graphite with up to tens of layers of thickness host two-dimensional electron states characterized by a large Berry curvature and a giant intrinsic magnetic moment. Note that stacking faults  garciaruiz19 in an rhombohedral film give rise to an additional strongly dispersing band that simply overlays the surface states spectrum. For example, in Fig. 4 we show the spectrum of an eleven layer rhombohedral (ABC) film with an ABA stacking fault at the top surface, which features a graphene-like Dirac band with velocity $v/\sqrt{2}$ and a small asymmetry gap, determined by $\delta$.

Semiconducting ABC graphitic films may be also used to create topologically-protected edge modes at interface regions with an inverted sign of $\Delta$, controlled,e.g., by an oppositely-directed displacement field. Similarly to a $\pm\Delta$ domain wall in bilayer graphene martin08 , a $\pm\Delta$ interface in $N$ layer ABC graphite would host $N$ co-propagating one-dimensional bands inside the spectral gap (with opposite direction of propagation in $K$ and $K^{\prime}$ valleys). Also, an atomic step in the film thickness may produce an isolated pair of edge states inside the gap (one in each valley with opposite directions of drift), but this feature will depend on the crystallographic orientation of the edge: it would be best developed for a zigzag termination of the additional layer, and it would be suppressed for the armchair edge due to valley mixing. Finally, as the gapful spectrum of an ABC film may be the result of many body effects leading to spontaneous spin/valley symmetry breaking lemonik10 ; jung11 , topological features of ABC graphite, enhanced by the $N\gg 1$ number of layers, would produce various possibilities for gapless edge modes in the system. For valley-symmetric magnetic phases, this would result in valley-current carrying domain walls, hosting $N$ one-dimensional channels with the opposite direction of drift in valleys $K$ and $K^{\prime}$. For phases with a valley antisymmetric order parameter, this would result in $2N$ electrical-current carrying edge modes at the physical edge of the sample. Overall, multilayer rhombohedral graphite films offer an attractive playground for studying giant topological effects in their electronic transport characteristics.

Methods

Hybrid $\mathbf{k}\cdot\mathbf{p}$ theory and tight-binding model. Our calculations are based on the hybrid $\mathbf{k}\cdot\mathbf{p}$ theory and tight-binding model (HkpTB). HkpTB combines the expansion of the electronic Hamiltonian for a single graphene layer in the lowest relevant order of momentum counted from the center of the $K$ and $K^{\prime}$ valleys, which we use in the form $\kappa=\xi p_{x}+ip_{y}+iexB_{z}/\hbar$, with interlayer hopping that takes into account the lattice arrangement for rhombohedral graphite displayed in Fig. 5a. Here we use a $2N$ component basis $(\psi_{\mathrm{A1}},\psi_{\mathrm{B1}},\psi_{\mathrm{A2}},\psi_{\mathrm{B2}},\cdots,\psi_{\mathrm{AN}},\psi_{\mathrm{BN}})$ of $A$ and $B$ sublattice Bloch states for each valley and take into account all hoppings in the Slonczewski-Weiss-McClure parametrization of graphite mcclure69 ; dressel02 ; koshino09 , as marked in Fig. 5a. The Hamiltonian that determines the electronic bands in $N$-layer rhombohedral graphite reads

[TABLE]

where we use $2\times 2$ blocks

[TABLE]

Here $a_{z}$ is the interlayer spacing, $v=(\sqrt{3}/2)a\gamma_{0}/\hbar$, $v_{3}=(\sqrt{3}/2)a\gamma_{3}/\hbar$, $v_{4}=(\sqrt{3}/2)a\gamma_{4}/\hbar$ and $a$ is the in-plane lattice constant. The intra ($D_{n}$) and interlayer ($V_{n,n+1}$) elements take into account the in-plane components $B_{x}$, $B_{y}$ of an arbitrarily-oriented magnetic field via the Peierls substitution, generalising an approach applied previously to bilayer graphene kheirabadi16 . Matrix element $W$ describes next-neighboring-layer hopping in the vertical direction and $\delta$ represents a difference in energy between the dimer sites inside the crystal and non-dimer sites $A_{1}$ and $B_{N}$ in the outer layers.

The minimal tight-binding model, consisting of only nearest-neighbour intra- and interlayer hopping $\gamma_{0}$ and $\gamma_{1}$, approximates $X({\bf p})\approx(-\kappa^{\dagger}/p_{\mathrm{c}})^{N}$ in Eq. (3) yielding flat, isotropic surface bands $\epsilon_{\pm}\approx(p/p_{\mathrm{c}})^{N}$. Skew interlayer hopping $\gamma_{3}$ and next-nearest-layer hopping $\gamma_{2}$ produce trigonal warping as shown in the dispersion Fig. 1a and the Fermi surface plots in Fig. 1b. Skew interlayer hopping $\gamma_{4}$ and parameter $\delta$ break electron-hole symmetry in the energy spectrum. For numerical diagonalization of Eq. 8 we use the following values of tight-binding parameters kuz09 $\gamma_{0}=3.16\,$eV, $\gamma_{1}=0.381\,$eV, $\gamma_{3}=0.38\,$eV, $\gamma_{4}=0.14\,$eV, $\gamma_{2}=-0.020\,$eV dressel02 and $\delta=0.0\,$eV.

The band structure near the valley center is plotted in Fig. 5b, obtained by diagonalization of the Hamiltonian in Eq. 8 for $N=10$ layers. There are $N-1$ conduction/valence bands with apex at $p=0$ and energy $\pm\gamma_{1}$ arising from the interlayer hybridization of carbon orbitals within the middle of the graphite film. Additionally, there are two flat bands at the conduction/valence band edge with wave functions made of $p_{z}$ carbon orbitals on sites $A_{1}$ and $B_{N}$ in the surface layers (the ones that don’t find a nearest neighbor in the next layer to hybridize into a dimer state). Zooming into the dispersion of the flat bands at low energy reveals semi-metallic behavior in Fig. 1a and parametric regimes identified in Fig. 1b. The full band model was used to compute the anomalous Hall coefficient using the formal definitions of $\Omega^{(\pm)}(\mathbf{p})$ and $\mathfrak{m}_{z}(\mathbf{p})$ with the results in Fig. 5c; Fig. 5d indicates how features in the low-energy dispersion are related to peaks in the anomalous Hall effect. Fig. 5e, f shows the anomalous Hall effect for a smaller number of layers with different asymmetry values.

The simplified $2\times 2$ (two band) Hamiltonian in Eq. 1 was derived using a Schrieffer-Wolff transformation which eliminated all high energy conduction and valence subbands projecting Hamiltonian $\hat{H}_{N}$ onto a basis of $(\psi_{\mathrm{A1}},\psi_{\mathrm{BN}})$ Bloch states using perturbation theory in intra- and interlayer hoppings up to the $N$th order with effective mass $m_{\star}^{-1}=(2v^{2}/\gamma_{1})(2\gamma_{4}/\gamma_{0}+\delta/\gamma_{1})\sim(0.4m_{0})^{-1}$ where $m_{0}$ is the free electron mass. The energy spectrum predicted by the two band model Eq. 3 is compared with that of the full HkpTB Eq. 8 in Fig. 6a for zero magnetic field and zero asymmetry, in Fig. 6b as a function of perpendicular magnetic field $B_{z}$ with zero asymmetry, and in Fig. 6c for zero magnetic field and finite asymmetry $\Delta=40\,\rm meV$.

To study the Landau level spectra in a perpendicular magnetic field ${\bf B}=(0,0,B_{z})$ we use the Landau gauge for the vector potential ${\bf A}=(0,B_{z}x,0)$. Then, $\kappa$ and $\kappa^{\dagger}$ transform into raising and lowering operators for magnetic oscillator states $\phi_{n}$ in valley $K$ mcc06 ; koshino09 (vice versa in valley $K^{\prime}$) with $\kappa\phi_{n}=i(\hbar/\lambda_{B})\sqrt{2(n+1)}\phi_{n+1}$ and $\kappa^{\dagger}\phi_{n}=-i(\hbar/\lambda_{B})\sqrt{2n}\phi_{n-1}$ with $\lambda_{B}=\sqrt{\hbar/(eB_{z})}$. Numerical diagonalization of the Hamiltonian in Eq. 8 is performed in a basis with a series of oscillator states for each sublattice component, $(\phi_{0},\phi_{1},\cdots,\phi_{N_{0}+n-1})$ for $A_{n}$ and $(\phi_{0},\phi_{1},\cdots,\phi_{N_{0}+n})$ for $B_{n}$, where $N_{0}=250$ is a cutoff, sufficient for convergence for $B>0.5\,$T. The result of the exact diagonalization is shown in Fig. 1d and e, indicating that the two-band model is sufficient to catch the qualitative features in the behavior of LLs at small magnetic fields, but with substantial quantitative deviations developing at $B_{z}>8\,$T.

Estimate of the bulk band gap. For an infinite number of layers (3D rhombohedral graphite), the bulk bands are gapless at $p=p_{\mathrm{c}}$ (the position of the Dirac point follows a continuous spiral as a function of the out-of-plane wave vector heik11 ; ho13 ; ho14 ; ho16 ; hyart18 ). In order to estimate the magnitude of the bulk gap for finite $N$, we consider the minimal model (including only $\gamma_{0}$ and $\gamma_{1}$) with in-plane momentum of magnitude $p=p_{\mathrm{c}}$ and directed in the $x$ direction only, i.e. $\mathbf{p}=(p_{\mathrm{c}},0)$. Then, the $2N\times 2N$ Hamiltonian [Eq. 8 with $\gamma_{2}=\gamma_{3}=\gamma_{4}=\delta=U_{i}=|\mathbf{B}|=0$] may be written as

[TABLE]

This is simply the tight-binding model of a linear chain with $2N$ sites and nearest-neighbor hopping $\gamma_{1}$ between every site. The solutions are $E=2\gamma_{1}\cos(\pi j/(2N+1))$ for $j=1,2,\ldots,2N$ which describes an electron-hole symmetric series of energies. The positive energy closest to zero ($j=N$) which describes the surface state is $E_{0}=2\gamma_{1}\cos(\pi N/(2N+1))=2\gamma_{1}\sin(\pi/[2(2N+1)])$ and $E_{0}\approx\pi\gamma_{1}/(2N)$ for $N\gg 1$. The positive energy which is the next closest to zero ($j=N-1$) which represents the lowest-energy bulk band is $E_{1}=2\gamma_{1}\cos(\pi(N-1)/(2N+1))=2\gamma_{1}\sin(3\pi/[2(2N+1)])$ and $E_{1}\approx 3\pi\gamma_{1}/(2N)$ for $N\gg 1$. Thus, we estimate the bulk band gap $\Delta_{\mathrm{bulk}}=2E_{1}$ to be

[TABLE]

for $N\gg 1$. We estimate $\Delta_{\mathrm{bulk}}\approx 120\,$meV for $N=30$ which seems to be in agreement with the numerical calculation of Henni et al henni16 , and we estimate $\Delta_{\mathrm{bulk}}\approx 36\,$meV for $N=100$.

Semiclassical quantization in a magnetic field and magnetic breakdown. For a given cyclotron orbit $C$ with area $S(C)$ in reciprocal space xiao10 , the semiclassical quantization condition is

[TABLE]

and $\Gamma(C)$ is the Berry phase of orbit $C$ defined as the $E=\epsilon(\mathbf{p})-\mathbf{B}\cdot\mathfrak{m}(\mathbf{p})$ contour. We find that for the $K$ valley and $\Delta>0$, the Landau levels are described by Eq.(10) with $n=0,1,2...$, while for $\Delta<0$ the levels are described by Eq.(10) with $n=1,2...$, see Fig. 2a. Similarly, taking into account that parallel magnetic field creates an effective top-bottom asymmetry gap $\sim[\mathbf{p}\times\mathbf{B}_{\parallel}]_{z}$ (see Eq.(3)), for $B_{y}=20\,$T, Landau levels of the $p_{x}>0$ conduction band pocket are described in the K valley by Eq.(10) with $n=1,2...$, while the two pockets with $p_{x}<0$ are described by Eq.(10) with $n=0,1,2,..$. These rules are swapped in the other valley (or, for $B_{y}=-B_{y}$ ), see Fig. 2b.

As is evident from Fig. 2b for $\Delta=10\,\rm meV$ and Fig. 7a for $\Delta=5\,\rm meV$, there is a significant mixing of LLs of the valence and conduction bands for valley $K$, while such mixing is almost absent in valley $K^{\prime}$ (this is reversed for $\Delta\to-\Delta$). At zero magnetic field, the states in the valence and conduction bands are orthogonal for a given quasimomentum, while a semiclassical wave-packet drifting along the Fermi contour at non-zero magnetic field has non-zero spread of quasimomenta around the Fermi contour, leading to non-zero projection onto another band. Thus, there are two main factors that contribute to interband breakdown: (i) the extent to which the states in the two bands are non-orthogonal when taken at nearby points in momentum space, and (ii) the time that the wave-packet spends in the part of the trajectory where breakdown is most probable.

To study the degree of non-orthogonality of valence and conduction band states at different momentum points, we plot the direction and amplitude of the maximal gradient ${\langle}u_{\mathrm{v}}(p)|\nabla_{p}|u_{\mathrm{c}}(p){\rangle}$ with blue vectors in Fig. 7b. The potential breakdown region corresponds to large gradients connecting the two Fermi surfaces (as illustrated by green arrows in Fig. 7b). The second aspect, which appears to be crucial for the magnetic breakdown, is related to the value and the sign of Berry curvature. This can be related to a notable decrease in the momentum-space velocity of the wave-packet as it passes an area of large Berry curvature. Indeed, the semiclassical equations of motion for an electron wave-packet in band $n$ are xiao10 ; niu99 :

[TABLE]

where ${\bf v}(p)={\bf\nabla}\epsilon(p)$ and the energy of the wave-packet is offset by the orbital magnetic moment. As seen in Fig. 7b, there is a large Berry curvature and small velocity near the two points of every valence band Fermi contour (these points are near the saddle point of the dispersion, similarly to the discussion in glazman18 ), while the conduction band Fermi contours have higher velocity and are located at low Berry curvature regions. Looking at the sign of Berry curvature of the valence band, we see that electrons in the valence band have higher probability to be in the breakdown-prone region for the $K$ valley due to lower momentum space velocity (larger denominator in Eq. (11) due to $e\Omega^{\mathrm{v}}B_{z}/\hbar>0$). This qualitatively explains the difference of breakdown patterns in the two valleys. Empirical evidence from numerical studies of LL spectra at different layer numbers and $\Delta$ indicate that magnetic breakdown starts when $e\Omega^{\mathrm{v}}B_{z}/\hbar\gtrsim 0.6$ at the valence band Fermi line near the breakdown region indicated with green arrows. Despite many years of studies glazman18 , we are not aware of any similar cases studied before. Since Eqs. (11, 12) were derived in the limit $|e\Omega B_{z}/\hbar|\ll 1$, which is no longer valid near the breakdown region, further theoretical study of this phenomemon is warranted.

**Data availability

**The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments

**We thank Servet Ozdemir, Vladimir Enaldiev, Yanmeng Shi, Jun Yin, and Artem Mishchenko for useful discussions. We acknowledge support from EU Graphene Flagship Project, EPSRC grants EP/S019367/1, EP/P026850/1 and EP/N010345/1, EC Project 2D-SIPC, and ERC Synergy Grant Hetero2D. M.K. acknowledges the financial support of JSPS KAKENHI Grant Number JP17K05496.

**Author contributions

**S.S., E.M., M.K. and V.I.F. contributed to the formulation of the problem and development of the theory. S.S. performed the numerical calculations and prepared the figures. S.S., E.M. and V.I.F. wrote the manuscript.

**Competing interests

**The authors declare no competing interests.