Films of rhombohedral graphite as two-dimensional topological semimetals

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.


Introduction
Berry curvature [1,2] is a measure of the topological nature of non-trivial states appearing in materials ranging from spin-orbit-coupling driven topological insulators [3][4][5][6] to graphene [7][8][9]. In multivalley conductors, such as graphene, Berry curvature with different signs concentrated in different valleys [10] 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 [11] due to progress in fabricating and characterizing thin films [12][13][14][15][16][17]. 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 γ 1 ≈ 0.38 eV leads to a gapped electronic spectrum in the middle of a thin film of N layers (the bulk gap ≈ 3πγ 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 [6,18], to low-energy surface states.
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 ∼ p c where the dispersion is highly anisotropic (p c = γ 1 /v where v is the Dirac velocity of electrons in graphene determined by the intralayer carbon-carbon hopping parameter γ 0 ). The presence of spatial asymmetry between the surface layers, which may be controlled using an electric displacement field applied perpendicularly to the film [27,28] or which may arise from interaction-induced spontaneous symmetry breaking [29,30], 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 intraand interlayer carbon-carbon couplings (see Methods) in the full Slonzewski-Weiss-McClure (SWMcC) tight-binding model [11,22,25,31], 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 p counted in a zig-zag direction from the center of the valley K and normalized by p c = γ 1 /v. In contrast to some earlier studies [13,26], the dispersion in Fig. 1 takes into account both skew interlayer couplings γ 3 and γ 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 p c . The electron dispersion at p ∼ p c is highly anisotropic, due to trigonal warping effects generated by skew interlayer hoppings, with inverted orientation in valley K .
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 bruteforce diagonalization of a hybrid k · 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 k · p theory in the K and K valleys, combined with interlayer hoppings introduced in the spirit of a tightbinding model (see Methods). For a film with N layers, this involves finding eigenvalues and eigenstates of a 2N × 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 ∆ = 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 [25] or spontaneous symmetry breaking into a magnetic state [32][33][34], making ∆ spin and/or valley dependent as in the layer antiferromagnetic configurations discussed in the context of bilayer graphene [29,30] and experimentally observed in rhombohedral graphite for N = 3 [35,36] and N = 4 [15].
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 Here, m ∼ 0.
, a z is the interlayer distance and −e is the electronic charge. Parameter δ 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 κ and its Hermitian conjugate, κ † , we use the Landau gauge for the out-ofplane magnetic field B z , whereas, for B z = 0, κ = ξp x + ip y . The off-diagonal element, , arises from interlayer hoppings (skew hopping between neighboring layers γ 3 and 'vertical' next-layer hopping γ 2 in addition to γ 1 mentioned earlier) passing electrons from one surface to another [22], and the summation is taken over integers n i ≥ 0 such that n 1 + 2n 2 + 3n 3 = N .
Qualitatively, the spectrum ofĤ, ± (p) = p 2 /(2m )±|d|, reproduces the exact multilayer 2N × 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,∆ = 0, the spectrum ofĤ 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 closelypacked groups of 3 LLs (i.e. 12-fold degenerate) whose origin we trace to three mini-valleys forming at p ∼ p c sketched in Fig. 1b. At high magnetic fields,Ĥ generates N separate groups of 4-fold degenerate low-energy LLs, which were considered to be degenerate in previous studies [13,21] where the effects of γ 4 and δ were neglected, and X(p) resulted in a Berry phase N π 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 = 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, ∆, 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 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 ±∆/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 . 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 Ω (±) (p) and associated magnetic moment m z (p) [1,2], , whereas m z (p) is the same for the + and − bands.
Using the two-component model (1) for N 1 and neglecting trigonal warping (parameters γ 2 and γ 3 ), the Berry curvature Ω (±) for conduction/valence bands is given by which generalizes the result for finite ∆ determined in monolayer and bilayer graphene [2,10].
The orbital magnetic moment for both bands is We note the strong N 2 layer dependence and the fact that in-plane magnetic field B can contribute to the gap∆(p).
Whereas Ω (±) (p) and m z (p) are peaked at p ≈ 0 in monolayer and bilayer, for rhombohedral graphene with N 1 they are peaked at p ≈ p c . Moreover, trigonal warping (parameters γ 2 and γ 3 ) introduces anisotropy such that Ω (±) (p) and m z (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 m z across momentum space in both valleys for ∆ = 10 meV.
This magnetic moment leads to a valley splitting, 2m z B z , (e.g. of the LL converging towards the valence band edge, Fig. 2a), that can be characterized by a factor g v = m z /µ B which can be as large as g v ∼ 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 the presence of in-plane electric field E, carriers in bands with Berry curvature experience a drift [2,37], v ± (p) = ∇ p ± (p) + (e/ )E ×ê z Ω (±) (p), so that, together, Berry curvature and valley polarization produce an anomalous contribution to the Hall effect, Here i lists all Fermi lines L i in valley K for both ± bands, and the linear integral is taken in the anticlockwise direction along each Fermi line, L i . When combined with the classical kinetic Hall coefficient,  Fig. 3a, both tenand twenty-layer films with ∆ = 40 meV are gapful 2D semiconductors (2DSC), so that the states with the maximium m z Ω near the valence band edge (Fig. 1c) are reached at a relatively small p-doping, resulting in a hump in σ xy (n e ) at n 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 m z Ω is the largest is set above the Fermi level for undoped materials, shifting the hump in σ xy to n e > 0. A double-peak structure is caused by the convolution of the m z Ω product and the density of states with a pronounced Van Hove singularity (pointed by down arrows) in Eq.(4). At higher electron dopings, the Fermi level reaches the maximum of m z Ω 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, δσ 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 ω and absorption coefficient g(ω), plotted in Fig. 3d, where ζ = ±1 stands for left/right handed circular polarization of the pump (radiation is approaching from the top) with power density W , Ω − K is the Berry curvature of the valence band at valley K, and τ rec is the lifetime of photo-excited holes at the top, p 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 [38] 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/ √ 2 and a small asymmetry gap, determined by δ.
Semiconducting ABC graphitic films may be also used to create topologically-protected edge modes at interface regions with an inverted sign of ∆, controlled,e.g., by an oppositelydirected displacement field. Similarly to a ±∆ domain wall in bilayer graphene [39], a ±∆ 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 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 [40,41], topological features of bands in N -layer rhombohedral graphite readŝ where we use 2 × 2 blocks Here a z is the interlayer spacing, v = ( 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 [42].
Semiclassical quantization in a magnetic field and magnetic breakdown. For a given cyclotron orbit C with area S(C) in reciprocal space [2], the semiclassical quantization condition is and Γ(C) is the Berry phase of orbit C defined as the E = (p) − B · m(p) contour. We find that for the K valley and ∆ > 0, the Landau levels are described by Eq. (8)  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 u v (p)|∇ p |u c (p) 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 [2,47]:ṗ = −e[v × n z ]B z − eE 1 + eΩB z / (9) r = v + eE × n z Ω 1 + eΩB z / where v(p) = ∇ (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 [48]), 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. (9) due to eΩ v B z / > 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 ∆ indicate that magnetic breakdown starts when eΩ v B z / 0.6 at the valence band Fermi line near the breakdown region indicated with green arrows. Despite many years of studies [48], we are not aware of any similar cases studied before. Since Eqs. (9, 10) were derived in the limit |eΩB z / | 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.