Precessing anisotropic Dirac cone and Landau subbands along a nodal spiral

We derive a three-dimensional Dirac-cone structure composed of tilted anisotropic Dirac cones around spirally located Dirac points. The Dirac points form a nodal spiral in momentum space due to accidental degeneracy, which can be realized in rhombohedral graphite. Under the interlayer electron hoppings, the Dirac cone varies in orientation and shape along this Dirac-point spiral, like a magic gyro precessing and deforming with time. The cone precession is governed by the hopping along the rhombohedral primitive unit vectors. In a perpendicular magnetic field, the electron orbits are related to different Landau level energies under the same quantization condition. The Landau subbands are thus characterized by a dispersion factor in addition to the zero-mode spectrum determined by the Dirac-point spiral.

high-symmetry axes as required by the hexagonal lattice symmetry [22]. Regarding the accidental nature, moving DPs and anisotropic Dirac cones were deduced by tuning the electron hoppings in monolayer graphene or other biparticle 2D systems [23,24]. It is therefore of our interest to analyze the effects of the electron hoppings in RG, so as to understand in general how a nodal spiral and the surrounding Dirac cones behave and lead to the LSs.
In section 2, we describe a full TB Hamiltonian for RG, based on the rhombohedral primitive unit cell, which includes the nearest neighboring vertical and non-vertical interlayer hoppings up to next-nearest neighboring layers. In sections 3 and 4, we explicitly analyze the existence and energy dispersion of the accidental DP spiral under the included interlayer hoppings, and then derive the Dirac cones from and around the DPs. Remarkably, the Dirac cone is tilted anisotropic, varying in orientation and shape along the DP spiral like a magic gyro precessing and deforming with time. The cone precession is due to the hopping along the primitive unit vectors. In comparison, a minimal model including only the nearest interlayer hopping results in a dispersionless DP spiral and all identical isotropic Dirac cones [16]. In section 5, we further figure out that the isoenergetic surface has 2D contours with various enclosed areas along the DP spiral and, therefore, the electron orbits are related to various LL energies under the same Onsager quantization condition [4]. The resulting LSs are commonly superimposed on the zero-mode LS and further characterized by a dispersion factor, reflecting the DP dispersion and the area variation, respectively. To our knowledge, this is the first derivation of the Landau quantization for RG using a full TB Hamiltonian. The calculation yields results in agreement with the numerical diagonalization [25], which can be explained neither by the minimal model nor by just considering the DP energy dispersion [18,26].

The model
The lattice of RG is depicted in figure 1(a) [14], with an atomic bond length of b = 1.42 Å, an interlayer distance of d = 3.37 Å and three primitive unit vectors a 1,2,3 associated with the c-axis 3dẑ = 3 i=1 a i . The figure shows the rhombohedral primitive cell with two atoms A and B and a three times larger (in volume) hexagonal cell, where the cyclic 1-2-3 sequence denotes the periodic ABC stacking. Figure 1(b) is a side view of the lattice with the atoms cyclically labeled. In figure 1(c), the rhombohedral Brillouin zone (RBZ) is plotted [14], together with a triply extended irreducible wedge belonging to three stacked hexagonal Brillouin zones (BZs); the inset of figure 1(c) illustrates their zone folding relation. The hexagonal BZs have six H -K -H edge extensions, which do not coincide with any high-symmetry points of the RBZ. It is noted in figures 1(a) and (c) that there are centers of symmetry in the primitive cell and the RBZ, revealing space-inversion and time-reversal symmetries of RG. Figure 1(d) is the top view of the lattice, where both the unit cells are projected onto the basal plane. For the hopping integrals represented in figure 1(b), general rules with i, j denoting the cyclic labels are given by β 0 (i = j), β 1 (i = j − 1), β 2 (i = j + 2), β 3 (i = j + 1) and β 5 (i = j − 2) between atoms A i and B j , and β 4 (| j − i| = 1) and β 5 (| j − i| = 2) between A i (B i ) and A j (B j ). In particular, β 4 refers to the hopping along a 1,2,3 .
Firstly, we make a continuum approximation for the low-energy band structure near the H -K -H hexagonal edge extension specified in figure 2(a). Here the edge only serves as a reference line, from which the in-plane wave-vector components (k x , k y ) are measured. Based on the two TB Bloch states |A and |B , the Hamiltonian is represented by a 2 × 2 matrix whose elements read where we have scaled the perpendicular wave-vector components k z by 1/d and defined 3,4,5). This chiral Hamiltonian characterizes the inversion symmetry of RG. Similar representation is easy to achieve about the H -K -H edge extension specified in figure 2(a), with φ replaced by −φ in equation (1). Referring to figure 1(d), the Hamiltonian matrix contains an interplay of those hopping processes. The hopping integral values should be consistent with firstprinciples calculation results [11], requiring β 0 = A|H|B <0 for the DP locations (described below). In both the present analysis and the following Landau quantization, the hopping integral values (in the unit of eV) transformed from the Slonczewski-Weiss-McClure (SWMcC) parameters are used [27]: β 0 = −2.73, β 1 = 0.32, β 2 = −0.0093, β 3 = 0.29, β 4 = 0.15 and β ( ) 5 = 0.0105. In equation (1), the off-diagonal elements can be written as where β 2 and β 5 are neglected in comparison with β 1 and v 3 , respectively. The anisotropic linearity of | f | in k indicates that the same chirality as in monolayer graphene while the presence of β 1 hopping dominates an offset of the DP from the hexagonal edge. Moreover, the identical diagonal elements cause the Dirac cones to be gapless, and their linearity in k indicates the possibility of band tilting.

Dirac-point spiral
The DP location (k D , φ D ) is conveniently described with respect to the hexagon for each k z . It is determined by band degeneracy, i.e. H 12 = 0 in equation (1) near the H -K -H edge extension, for which the effects of the interplaying hoppings with different phases are revealed by From equation (2), there is only one DP deviating from the specified hexagon corner for a fixed k z . This is a result of the corner not being a high-symmetry point in the RBZ [14]. Neglecting β 2 and v 5 , the DP is given by up to first-order perturbation O(v 3 /v 0 ). The magnitude of k D is far smaller than the size of the K -line. Within the minimal model including only β 1 as the interlayer hopping, the DP spiral as a whole lies on a cylindrical surface of radius β 1 (v 0h ) −1 with the spiraling angle φ D synchronizes with −k z . Similar analysis can be performed for the H -K -H edge extension, around which φ D synchronizes with k z and the chirality of the DP spiral is reversed. Due to the interlayer hoppings other than β 1 , the spiral form persists but becomes non-cylindrical.
With respect to the RBZ (figure 1(c)), there are six, namely, three H -K -H and three H -K -H extensions. The 2D projections of the six surrounding DP spirals are depicted in figure 2(a) in the frame of the hexagons. Each of them runs inside the RBZ for −π/3 k z π/3 and into the neighboring RBZs for other k z . One DP exists at the K ( ) -line for k z = 0, another one, at the H ( ) -L lines for k z = ±π/3 and the others are located in between. We then derive the energy dispersion in terms of the polar coordinates (q, θ) measured from the DPs. The coordinate transformation runs as Obviously, θ → −θ as φ → −φ when changing between chiralities. Ignoring β 2 and v 5 , the electron energy at the DP described in equation (3) is given by The same result can be derived for all the hexagonal egdes. Equation (4) reveals the DP energy dispersion with k z , indicating that RG is a semimetal. The DP spiral is surrounded by a sausagelink Fermi surface [16]. Near k z = 0 (±π/3), there is an electron (hole) pocket; at k z = ±π/6, the pockets shrink to the DP, where the Fermi energy E F = E D (±π/6) = 0 is determined.

Three-dimensional Dirac-cone structure
The eigenvalues of equation (1) near the specified H -K -H edge extension are obtained by E = H 11 ± |H 12 |, which describe the low-energy band structure of RG [16]. In terms of (q, θ ), the Dirac cone immediately approaches up to O(v 3 /v 0 ), where the ± signs refer to the upper and lower half Dirac cones, respectively, and v 5 is neglected compared with v 4 . When the v 3 and v 4 terms are ignored, all the Dirac cones are isotropic linear and identical to that of monolayer graphene with a Fermi velocity v 0 except that they are centered about the dispersing DPs. The presence of the v 3 and v 4 terms causes k z -dependent anisotropy. The 2D isoenergetic contours around the DPs are illustrated in figure 2(b) for k z = 0, π/6 and π/3. They have energies = E − E D according to equation (5).
The phase factor of the v 3 term does not change under a phase shift θ → θ + π . This term causes a rotation of the Dirac cone, shown in figure 2(b), as k z varies and the rotation angle is given by k z /2. However, the phase factor of the ±v 4 term has a reversed sign under that phase shift. This term tilts the Dirac cones with respect to the k z -axis. With both the v 3 and v 4 terms in equation (5), tilted anisotropic Dirac cones are obtained along the DP spiral. The -contours are outlined in figure 2(b). Furthermore, the 2D projections of the cone axes are marked by the dashed-green lines with the angles θ t = −k z + π since the extremals of the ±v 4 term lie at θ = −k z + π and −k z . The cone axis, angled from the k z -axis by α t , precesses clockwise with θ t varying from 2π to 0 in pace with k z from −π to π, as shown in figure 2(c). Remember that ±v 4 is attributed to β 4 , the electron hopping along the rhombohedral primitive unit vectors. Like a magic gyro precessing and deforming with time, the Dirac cone varies in orientation and shape with k z along the DP spiral. It can shown that the chirality of the precessing tilted Dirac cone is reversed when changing the chirality of the DP spiral. The area S( , k z ) enclosed by the -contour is calculated from 2π 0 dθ Q 0 q dq, with Q determined by : with θ = θ − k z /2. In principle, equation (6) can be integrated by the trigonometric substitution tan θ /2 = t because the integrand is a rational function of sin θ and cos θ . Alternatively, we expand equation (6) in powers of v 4 and obtain S( The integrations are performed through certain trigonometric substitutions to obtain where symmetries of the involved functions are called for simplifying the calculations of I 2 and I 3 . Accordingly, equation (6) is approximated by with F(k z ) = 6(v 4 /v 0 ) 2 [1 + 2(v 3 /v 0 )cos 3k z ], which yields nearly the same results as the direct numerical integration does. Equation (7) states that the enclosed area of the -contour is independent of the chirality of the cone precession. S( , k z ) decreases as k z increases from 0 to π/3 for an isoenergetic surface with . That is, grows with increasing k z for a constant area S( , k z ). The variation of S( , k z ), as illustrated in figure 2(b), provides a physical picture of the following Landau quantization.

Onsager quantization
In a uniform magnetic field B 0ẑ , the Onsager quantization rule is used [4,28]. It is given by with n being the LL index. The zero phase shift in equation (8) for RG results from the same electron chirality and Berry phase as in monolayer graphene, regardless of the anisotropy of the Dirac cones [4,29]. Ignoring v 3 and v 4 in equations (7) and (8), the quantized energy ± (0) (n, B 0 ) = ±(2ehn B 0 ) 1/2 v 0 is independent of k z and identical to the LL spectrum in monolayer graphene with the Fermi velocity v 0 , as if RG were a quasi-2D stack of independent graphene layers. The consequent LS spectrum reads where the signs ± refer to the unoccupied and the occupied LSs, respectively. The LS dispersions arise from only the superimposition on the zero mode (n = 0), which amounts to the DP energy E D (k z ).
Introducing v 3 and v 4 gives the quantized energy The dispersion factor contained in the braces arises from the increase of the LL energy with increasing k z for a constant quantized S( ). The LS spectrum is then obtained by superimposing Equation (10) in combination with equations (4) and (9) shows that the interlayer hoppings bring about characteristics beyond the minimal model. The electron-hole asymmetry about E F = 0 is induced by the DP dispersion. For field strengths allowing LS bulk gaps, charge compensation is simply achieved because only the n = 0 LS intersects E F = 0. The intersection occurs at k z = π/6, about which the cos 3k z function contained in E D (k z ) is odd. The renormalized Fermi velocity is determined at k z = π/6, given by which is smaller than v 0 of monolayer graphene. Equation (11) also indicates that v 3 and v 4 reduce the LS spacings. The dispersion factor in equation (9) leads to opposite trends for the occupied versus the unoccupied LSs, with the former dispersing in the same trend as the n = 0 LS (equation (4)) and the latter dispersing oppositely. Moreover, the k z dispersion is enhanced with increasing n. Due to this dispersion enhancement, the bulk gap between two neighboring higher-mode LSs can be closed.

Discussion and conclusions
The LSs are calculated from equation (10) for B 0 = 30 T and plotted in figures 3(a)-(c). Also plotted are the numerical diagonalization results for which β 2 and β 5 are additionally taken into account [25]. In contrast to both AA-and AB-stacked graphites [30,31], much weakly dispersive low-lying LSs are obtained and only the n = 0 LS intersects with E F = 0. Figures 3(a) and (b) depict the LSs ranging from zero energy to ∼±0.2|β 0 | ±550 meV. The Onsager quantization and the numerical diagonalization agree with each other very well from n = 0 to 4 and their agreement holds approximately up to ∼±0.27|β 0 | ±730 meV, as shown in figure 3(c). The difference grows slightly in the increase of n as a result of the continuum approximation in use. With the current field strength, the bulk gaps vanish for LSs as high as n = 25 (27) in the occupied (unoccupied) spectrum (not shown). Figures 3(a)-(c) also include the LSs obtained from the minimal model, which are dispersionless and have a larger Fermi velocity. In figure 3(d), the quantized ± (n, B 0 , k z ) spectrum is plotted within the same energy range as in figure 3(c), where the dispersion enhancement and opposite trends due to the dispersion factor are clearly seen. The main characteristics of the LSs are helpful in understanding the magnetoelectronic properties. The LS spectrum makes up a one-dimensional magneto-electronic structure since they disperse only in the k z -direction. The width w D = 2E D (k z = 0) 10.5 meV of the n = 0 LS (figure 3(a)) can be signified by the low-temperature electronic specific heat, which majorly comes from the intra-LS thermal excitations [32]. All the LSs could lead to the peculiar magneto-optical spectra. The optical transition between the LSs of m and n must obey a specific selection rule, |m − n| = 1, for the linear Dirac cone at a fixed k z . The spectral intensity, which depends on k z , is expected to exhibit absorption peaks for k z = 0 and π/3 with the highest density of states. The magnetoabsorption measurements are available [33][34][35], for identifying the characteristics of interband transitions. A representative transition due to the n = 13 and 14 LSs is shown in figure 3(d), where double peaks with a separation ∼9.5 meV can be resolved from the infrared spectroscopy [33]. Besides, probing the cyclotron resonance near E F = 0 is feasible by measuring the far-infrared spectroscopy [34,35]. Moreover, RG has a LS bulk gap due to the n = 0 and 1 LSs in very weak magnetic fields ∼ B 0 w 2 D /2eh ( 0.11 T). This is in contrast to the strongly dispersive AA-and AB-stacked graphites, for which unattainable field strengths are required to open the bulk gap due to the lowest LSs [26,36]. RG is a candidate system for the existence of the 3D QHE [13]. In conclusion, we have derived a 3D Dirac-cone structure in which the Dirac cone is tilted anisotropic and, like a magic gyro, precesses and deforms with k z along the DP spiral. In association with the DP spirals, there are tilted Dirac cones precessing with opposite chiralities in reciprocal space. The electron orbits for different k z are related to different LL energies under the same quantization condition, giving rise to a dispersion factor in the LSs in addition to the zero-mode spectrum determined by the DP spiral. RG is a semimetal with the DP spiral appearing accidentally in the bulk momentum space. We have shown that the interlayer electron hoppings included in a full TB Hamiltonian do not lift the DP spirals while they bring about significant characteristics; in particular, the β 4 hopping along the rhombohedral primitive unit vectors causes the cone precession. The analysis provides a general approach to investigating how the interlayer hoppings can affect 3D layered systems that possess nodal lines such as the DP spirals.