Topo-fermiology

The modern semiclassical theory of a Bloch electron in a magnetic field now encompasses the orbital magnetic moment and the geometric phase. These two notions are encoded in the Bohr-Sommerfeld quantization condition as a phase ($\lambda$) that is subleading in powers of the field; $\lambda$ is measurable in the phase offset of the de Haas-van Alphen oscillation, as well as of fixed-bias oscillations of the differential conductance in tunneling spectroscopy. In some solids and for certain field orientations, $\lambda/\pi$ are robustly integer-valued owing to the symmetry of the extremal orbit, i.e., they are the topological invariants of magnetotransport. Our comprehensive symmetry analysis identifies solids in any (magnetic) space group for which $\lambda$ is a topological invariant, as well as identifies the symmetry-enforced degeneracy of Landau levels. The analysis is simplified by our formulation of ten (and only ten) symmetry classes for closed, Fermi-surface orbits. Case studies are discussed for graphene, transition metal dichalchogenides, 3D Weyl and Dirac metals, and crystalline and $\mathbb{Z}_2$ topological insulators. In particular, we point out that a $\pi$ phase offset in the fundamental oscillation should \emph{not} be viewed as a smoking gun for a 3D Dirac metal.


I. INTRODUCTION
The semiclassical Peierls-Onsager-Lifshitz theory [1][2][3] connects experimentally accessible quantities in magnetic phenomena to Fermi-surface parameters of the solid at zero field.For example, field-induced quantum oscillations of magnetization [4] and resistivity [5] have become the leading method to map out the shape of the Fermi surface of normal metals [6,7] and superconductors [8]-this phenomenology has been coined "Fermiology" [6].
The semiclassical theory has been extended [9][10][11][12] to incorporate two modern notions: a wave packet orbiting in quasimomentum ðkÞ space acquires a geometric phase (ϕ B ) [13,14], as well as a second phase (ϕ R ) originating from the orbital magnetic moment of a wave packet around its center of mass [15].Further accounting for the well-known Zeeman coupling, [10,[16][17][18] to be the complete, subleading (in powers of the field) correction to the Bohr-Sommerfeld quantization rule for nondegenerate bands [2,3,19].λ is measurable as a phase offset in oscillations of the magnetization/resistivity in 3D solids, as well as in fixed-bias oscillations of the differential conductance in tunneling spectroscopy [20,21].
While it is conventionally believed that ϕ B ¼ 0 vs π distinguishes between Schrödinger and Dirac systems [22], we propose to view ϕ B =π as a continuous quantity that is sometimes fixed to an integer in certain space groups and for certain types of field-dependent orbits; moreover, while ϕ R vanishes for centrosymmetric metals without spin-orbit coupling (SOC), it plays an oft-ignored role in most other space groups.Our comprehensive symmetry analysis identifies the (magnetic) space groups in which λ=π is robustly integer valued-we will formulate λ as a topological invariant in magnetotransport, which is distinct from the traditional formulation of topological invariance in band insulators [23,24].We also extend our symmetry analysis to the multiband generalization of λ, with envisioned application to bands of arbitrary degeneracy (D); D ¼ 2 is exemplified by spin degeneracy.
Let us outline the organization of the main text.We begin in Sec.II by introducing the multiband quantization rule and describing how λ appears as the subleading phase correction.Experimental methods to extract λ are discussed in Sec.III; we present here generalized Lifshitz-Kosevich formulas for the oscillatory magnetization and density of states.These formulas extend previous works [3] in their applicability to orbits of any energy degeneracy and symmetry, including orbits in magnetic solids.In Sec.IV, we provide a general group-theoretic framework to identify solids for which λ takes only discrete values.In addition, our symmetry analysis identifies the symmetry-enforced degeneracy of Landau levels; where degeneracy is not enforced, symmetry may nevertheless constrain the possible splittings of Landau levels.We exemplify our symmetry analysis with several case studies in Sec.V, including graphene, transition metal dichalcogenides, surface states of topological insulators, and 3D Weyl and Dirac metals.In particular, we point out that a π phase offset in the fundamental oscillation should not be viewed as a smoking gun for a 3D Dirac metal.We recapitulate our main results in the concluding Sec.VI; a final remark broadens the applicability of our symmetry analysis to matrix representations of holonomy [25] in the Brillouin torus, also known as Wilson loops [26] of the Berry gauge field [13].

II. MULTIBAND BOHR-SOMMERFELD QUANTIZATION RULE
The quantization rule is derived from the effective Hamiltonian (H) that describes the low-energy dynamics of Bloch electrons in a field [9][10][11]27,28].In a basis comprising D field-modified Bloch functions at each wave vector, HðKÞ is a Weyl-symmetrized, matrix function of the kinetic quasimomentum operators, whose noncommutivity is manifest in K × K ¼ −ieB=c.H is asymptotically expandable in powers of the field: , where the leading-order term is the Peierls-Onsager Hamiltonian [1,2], while the subleading terms , respectively, encode the geometric phase, the orbital moment, and the Zeeman effect [9][10][11].
In the WKB approximation, the D-component vector wave function of H generalizes [29] the known singlecomponent solution developed by Zilberman [19] and Fischbeck [17] for a nondegenerate band.In the absence of breakdown [11,30,31], continuity of the vector wave function around a closed orbit (o) affords us the following quantization rule: We have assumed here that the field is oriented in ⃗ z, such that o is a contour of the band dispersion at fixed energy E and k z .By "closed orbit," we mean that o does not wrap around the Brillouin torus.o bounds a region in k ⊥ ≔ ðk x ; k y Þ-space with positive-definite area S. l ≔ ðℏc=ejBjÞ 1=2 above is the magnetic length, j an integer, and a ∈ Z D ≔ f0; 1; …; D − 1g.The Maslov correction ðϕ M Þ depends on the topology of the Fermi-surface orbit; e.g., it equals π for orbits that are deformable to a circle [32], but vanishes for a figure-of-eight orbit [33].To leading order in l −2 , Eq. ( 1) without λ a is a well-known result by Lifshitz and Onsager [2,3,19].The remaining term λ a is defined through the spectrum ðfe iλ a g D a¼1 Þ of the propagator (A) that is generated by H 1 over the cyclotron period.A may be expressed as a pathordered exponential (denoted exp): with o carrying a clockwise orientation; the above three one-forms represent contributions by H B 1 , H R 1 , and H Z 1 , respectively.The first one-form is the non-Abelian Berry connection [13,26] for the D-fold-degenerate subspace (henceforth denoted by P) and is defined by XðkÞ mn ≔ ihu mk j∇ k u nk i, with e ik•r u nk the Bloch function of a band labeled by n ∈ Z D .The multiband orbital magnetic moment is encoded in the second one-form with m, n ∈ Z D .ΠðkÞ ln ¼ ihu lk je −ik•r ½ Ĥ0 ; re ik•r ju nk i=ℏ are matrix elements of the velocity operator, with Ĥ0 the single-particle, translation-invariant Hamiltonian and r the position operator.For n ∈ Z D , Π nn ¼ v is the velocity of each band in P, and v ⊥ ≔ ðv 2 x þ v 2 y Þ 1=2 .X and Π in Eq. ( 3) comprise only off-block-diagonal matrix elements between P and its orthogonal complement.The third one-form in Eq. ( 2) is the well-known Zeeman coupling, with ℏσ z mn ðkÞ=2 the matrix elements of the spin operator S z , g 0 ≈ 2 the free-electron g-factor, and m the free-electron mass.While the definition of A presumes a basis choice fu nk g D n¼1 within P, one may verify that the eigenvalues of A are independent of this choice.By setting D to 1 in the above equations, the line integral of the three one-forms in Eq. ( 2) give respectively, as we have introduced in the second paragraph of the paper.
When Eq. ( 1) is viewed at fixed field, the discrete energetic solutions fE a;j g correspond to D sets of sub-Landau levels; within each set labeled by a, the difference between two adjacent levels (jE a;jþ1 − E a;j j) is approximately ℏω c ≔ 2π=ðl 2 j∂S=∂EjÞ evaluated at E a;j .When Eq. ( 1) is viewed at constant energy (e.g., the chemical potential μ), the discrete solutions (fl 2 a;j g) correspond to values of the field where Landau levels successively become equal to μ.In thermodynamic equilibrium, these are also the fields where Landau levels are suddenly depopulated with a periodicity: l 2 a;jþ1 − l 2 a;j ¼ 2π=SðμÞ, for each of a ∈ f1; …; Dg.This results in various oscillatory phenomena from which we may extract λ a .

III. GENERALIZED LIFSHITZ-KOSEVICH FORMULAS TO EXTRACT λ
In the de Haas-van Alphen (dHvA) effect [4], each extremal orbit (o) on the Fermi surface of a 3D metal is associated to an oscillatory contribution to the longitudinal magnetization (parallel to the field in ⃗ z): which is a sum of D sets of harmonics.Being valid in the degenerate (μ ≫ kT) and semiclassical (μ ≫ ℏω c ) limits, Eq. ( 4) is our generalization of the Lifshitz-Kosevich formula [3] to orbits of any energy degeneracy (D) and symmetry, including orbits in magnetic solids.In comparison, the commonly employed Lifshitz-Kosevich formula with a "spin reduction factor" [6,34] is only applicable to twofold degenerate orbits in solids with both time-reversal and spatial-inversion symmetries.All quantities on the right-hand side of Eq. ( 4) are evaluated on o, which may be electronlike or holelike; the sign of π=4 (in the argument of the sine function) is negative (positive) for a maximal (minimal) orbit.S zz is the double derivative of S with respect to k z , and we have introduced Dingle's damping factor [35] that depends on the quasiparticle's mean free time (τ).For D ¼ 1, the field-independent phase in the argument of the fundamental harmonic is sometimes referred to as the Onsager phase: If a Fermi surface has multiple extremal orbits, each extremal orbit additively contributes a term with the same functional form as Eq. ( 4).If two extremal orbits (o i and o iþ1 ) are symmetry related, they contribute oscillatory terms that are identical in the parameters fS; S zz ; ω c ; ϕ M g, but not necessarily for the λ-phase corrections.Generally, fλ i a g D a¼1 ¼ fAEλ iþ1 a g D a¼1 (defined modulo 2π), with the sign depending on the symmetry class of the orbit, as we will elaborate in Eq. ( 13) of Sec.IV.
In 2D metals, the analogous oscillatory formula is The field dependence of μ is negligible in the semiclassical limit for 3D metals [6], as well as for 2D surface states of 3D solids [36].In strictly 2D metals with a fixed particle density, field-induced oscillations in μ render the extraction of λ implausible.An alternative method to extract λ a is to measure the temperature-broadened, 3D density of states, defined as Here, gðEÞ is the density of states of 3D Landau levels (smoothened by the Dingle factor); f 0 T ðxÞ is the derivative of the Fermi-Dirac distribution function and approaches −δðxÞ in the zero-temperature limit.The oscillatory component of G may be expressed as a harmonic expansion: This quantity may be accessed via the scanning tunneling microscope (STM) [20,21,37] or by planar tunneling junctions [38,39].The oscillatory component of the differential conductance (dI=dV), averaged over the surface of a 3D metal at fixed bias voltage (V) for the STM, is directly proportional to δGðμ þ eVÞ in the absence of surfacelocalized states [40]; for the planar tunneling junction, no spatial averaging is needed for a sufficiently large junction size.The just-mentioned proportionality presupposes that (a) the tunneling matrix elements and the density of states of the STM tip are featureless in the energy range that is accessed by the bias voltage, and that (b) the tip and sample have thermally equilibrated [38,40].Landaulevel spectroscopy via the scanning tunneling microscope has already been reported [20,21,37], but we are not aware that the phase offset of the oscillations has ever been measured.Further details on the derivation of Eqs. ( 4)-( 8) are provided in Appendix A.
Let us discuss how to extract λ from dHvA data.For simplicity in presentation, we consider the magnetization oscillations contributed by a single orbit (extremal orbit in 3D metals).We assume that l 2 S, ℏω c , and the Dingle lifetime τ have already been extracted by standard techniques [41].Let us first consider either ω c τ ≪ 1 or kT ≫ ℏω c , such that the harmonic expansion is dominated by the fundamental (r ¼ 1) harmonic.If D ¼ 1, then the experimental data may directly be fitted to a single sine function offset by the Onsager phase [cf.Eq. ( 5)].If D ¼ 2 (e.g., spin degeneracy), the sum of two fundamental harmonics produces an equifrequency harmonic proportional to [16] ; ð10Þ defined to be invariant (modulo 2π) under λ j → λ j þ 2π; these formulas will be applied to a case study of Bi 2 Se 3 in Sec.V G.If jS zz j is not otherwise measurable (for 3D metals), it is not possible to fully determine λ 1;2 , owing to our ignorance of the amplitude of the fundamental harmonic; measuring the phase offset of the fundamental harmonic merely determines Θ, a suitably defined average of λ 1;2 .Θ alone does not completely characterize the non-Abelian transport within the two-band subspace.
In the interest of measuring individual values of λ 1;2 , we propose higher-field dHvA measurements of cleaner samples (ω c τ not ≪ 1) and at lower temperatures (kT not ≫ ℏω c ).In this regime, not just the fundamental but also higher (r > 1) harmonics are needed to accurately represent the dHvA data [41].Since fλ a g D a¼1 is encoded in the interference of multiple harmonics, fλ a g D a¼1 may be extracted without knowledge of the absolute amplitude of a single harmonic.
In metals with multiple extremal orbits, there may exist field orientations where all orbits are related by symmetry, which simplifies the fitting to the Lifshitz-Kosevich formula [42].Independent of the field orientation, there exists one simplification for any nonmagnetic metal: due to time-reversal (T) symmetry, the set of all fλg comprises only pairs that are invariant under inversion about zero (λ → −λ), which effectively halves the independent parameters that require fitting.We quote this result here to exemplify the utility of a symmetry analysis of fλg, which we explore in greater generality in the next section (Sec.IV); the above-mentioned constraint by T symmetry is elaborated subsequently in Sec.V E.

IV. SYMMETRY ANALYSIS OF THE λ PHASE
In certain (magnetic) space groups, λ (or P D a¼1 λ a =π for D > 1) is integer valued, owing to the symmetry of the extremal orbit.To identify these orbits and space groups, it is useful to distinguish ten symmetry classes for closed orbits-to each class we associate certain constraints for the propagator A and its spectrum (fλ a g D a¼1 ), as summarized in Table I.The goal of this section is identify the relevant symmetry class of closed orbit-for any physical system one chooses to study.Once this identification is made, the resultant symmetry constraints on fλ a g D a¼1 may be read off from the last column of Table I and verified experimentally from any of the methods detailed in Sec.III.This will be exemplified for several case studies in Sec.V.
Let us first restate the problem in simple terms: The dynamics of Bloch electrons immersed in ⃗ z are restricted to Brillouin two-tori (BT ⊥ ) of fixed k z .For a D-fold degenerate band subspace with dispersion εðkÞ, semiclassical motion occurs along (assumed) closed orbits defined by If multiple disconnected orbits exist within the same BT ⊥ , we assume they are sufficiently separated in k ⊥ -space that tunneling is negligible.Neglecting the subleading term (H 1 ) in the effective Hamiltonian, all Landau levels are at least D-fold degenerate, owing to the Onsager-Lifshitz quantization rule; here and henceforth, the "degeneracy of a Landau level" is defined in units such that the extensive degeneracy TABLE I.The first three columns distinguish among ten classes of elementary orbits.The map of k ⊥ under g is g∘k ⊥ ¼ ð−1Þ sðgÞ ǧ⊥ k ⊥ , with ǧ⊥ a 2 × 2 orthogonal matrix that represents the point-group component of g in the plane orthogonal to the field; sðgÞ ¼ 0 if g is purely a spatial transformation, and 1 if g inverts time.joj ¼ jg∘oj indicates that o is mapped to itself under g, modulo a change in orientation.uðgÞ distinguishes between proper and improper transformations on k ⊥ : ð−1Þ u ≔ det ǧ⊥ .The fourth and fifth columns describe how unitary matrices ḡ (that represent the symmetry g) constrain the propagator.Column six summarizes the constraints on λ a ; if there are none, we indicated this by Á Á Á.
of a Landau level (associated to a single spinless orbit) equals 1 [cf.Eq. ( A3)].For a subset of Landau levels, this zeroth-order degeneracy is enhanced to LD if a symmetry (g) constrains L disconnected orbits to have identical shape.We may ask if (and how) H 1 splits this LD-fold (or D-fold) degeneracy [43]; if L ¼ D ¼ 1, we ask if H 1 shifts the zeroth-order Landau-level spectrum at all.The answer to these questions depends on the class of symmetric orbit, which we proceed to analyze in full generality.

A. Tenfold classification of symmetric orbits
Since lattice translations trivially constrain A, we shall henceforth focus on symmetries (g) of the solid that correspond to nontrivial elements in the (magnetic) point group (P) of the solid; examples include (screw) rotations, (glide) reflections, spatial inversion, and time reversal.Technically speaking, magnetic point groups differ from point groups in that the former includes a symmetry that reverses time; however, this distinction is irrelevant to the following classification of symmetric orbits.Hence, we hereafter use "point group" democratically.
We are interested only in g that maps BT ⊥ to itself; such g corresponds to a subgroup (P ⊥ ) of P that generally depends on the field orientation, as well as k z .Any configuration of closed orbits in BT ⊥ may be divided into a disjoint set of elementary orbits fðg; O i Þg, where O i is defined to be the smallest, closed orbit configuration that is invariant under g.By "invariance," we mean that for every There are three topologically distinct mappings of k ⊥ ∈ o.The simplest is the identity map, where each k ⊥ in BT ⊥ (but not necessarily in the entire 3D torus) is individually invariant under g.Such mappings are labeled as class I, and all other mappings are of class II.We further distinguish between class-II mappings where g∘o is identical to o up to orientation (such mappings are labelled as class II-A), and mappings between disconnected orbits (labelled as class II-B).
There are two classes of class-I elementary orbits distinguished by whether g is purely a spatial transformation or otherwise includes a time reversal.We introduce a Z 2 index sðgÞ that equals 0 in the former and 1 in the latter.(Class I; s ¼ 0) is exemplified by BT ⊥ being a mirror/ glide-invariant plane, and ½I; s ¼ 1 by g ¼ Ti, which is the composition of time reversal (T) with spatial inversion (i); all class-I symmetries are order two.Class-II elementary orbits are likewise distinguished by whether g inverts time; they are additionally distinguished by whether g acts on k ⊥ as a two-dimensional rotation (u ¼ 0) or as a twodimensional reflection (u ¼ 1).Equivalently, given that o is clockwise oriented, uðgÞ distinguishes between symmetries that preserve (u ¼ 0) or invert (u ¼ 1) this orientation.In each of II-A and II-B, there are then four classes of elementary orbits distinguished by s, u ∈ Z 2 .This gives ten classes of elementary orbits in total, whose defining characteristics are summarized in the first three columns of Table I.
In class I and II-A, O i is composed of a single orbit o, which is self-constrained by g.In II-B, O i is composed of L disconnected orbits which are mutually constrained as g∘o i ¼ ð−1Þ u o iþ1 and o iþL ≔ o i .To clarify, o and fo j g L j¼1 are all clockwise oriented, and −o i denotes an anticlockwise-oriented orbit.L was introduced in the second paragraph of Sec.IV and is more precisely defined here as the smallest integer for which g L ∘k ⊥ ¼ k ⊥ for all k ⊥ ; generally, L divides the order (N) of g, e.g., L ¼ 3 and N ¼ 6 for the composition of T and a sixfold rotation.

B. Symmetry constraints on the propagator A
Column four summarizes how g constrains A (I,II-A) and fA j g (II-B), which are, respectively, the propagators for the self-constrained o and mutually constrained fo j g.The corresponding spectra of the propagators are denoted as fe iλ a g D a¼1 and fe iλ j a g D a¼1 .The unitary matrices ḡ that constrain these propagators form a projective representation [44,45] of the point group P ⊥ , as summarized in column five.Any g that inverts time (s ¼ 1) has the antiunitary representation ḡK, with K implementing complex conjugation; otherwise (s ¼ 0), g has the unitary representation ḡ.The relations in column five are closely analogous to the space-group relations [46] satisfied by g: Here, NðgÞ ∈ N is the smallest integer such that g N is a translation (t) by the lattice vector R, possibly composed with a 2π rotation (denoted e).Note the similarity in definition of N with L, as was defined in Sec.IVA; in general, L divides N. R is nonzero for nonsymmorphic symmetries such as screw rotations and glide reflections [46]; the translation t R is represented on Bloch functions with wave vector k by the phase factor e −ik•R (cf.column five).e in Eq. ( 11) is represented in column five by a phase factor [ð−1Þ F ] that a wave function accumulates upon a 2π rotation; F ¼ 0 (F ¼ 1) for integer-spin (half-integer-spin) representations.The former case is useful in analyzing solids with negligible spin-orbit coupling, as we will exemplify with a case study of graphene in Sec.VA.The constraints in columns four and five are derived from the symmetry transformation of the Berry connection [47] and the one-form A • dk of Eq. (3).The latter may be expressed through Hamilton's equation (ℏ _ k ¼ −jejv × B=ℏc) as −M • Bdt=ℏ, with MðkÞ an orbital moment that transforms under g like the spatial components of a (3 þ 1)dimensional pseudovector: Mj g∘k ¼ ð−1Þ sðgÞ det½ǧḡK sðgÞ ðǧMÞK sðgÞ ḡ−1 j k : ð12Þ ǧ here is a 3 × 3 orthogonal matrix representing the pointgroup component of g, and

C. Symmetry constraints on λ
From taking the determinant of each equation in column four, we derive constraints on λ a that are summarized in column six.
Three of six classes in [I,II-A] are characterized by the reality condition e i P a λ a ∈ R.This implies λ ¼ 0 or π for a nondegenerate band; i.e., the orbit, respectively, encircles an even or odd number of Dirac points.λ ¼ π is exemplified by the Dirac surface state [48] of the Z 2 topological insulator Bi 2 Se 3 [49][50][51][52]; an orbit on the Dirac cone is selfconstrained by time-reversal (T) symmetry (class II-A, u ¼ 0, s ¼ 1).λ ¼ π may result solely from crystalline symmetry, as exemplified by a mirror-symmetric orbit (class II-A, u ¼ 1, s ¼ 0) that encircles a surface Dirac cone of the topological crystalline insulator SnTe [53] (further elaborated in Sec.V B).Finally, λ ¼ π can be protected by a composition of T and crystalline symmetry (class I, s ¼ 1), as exemplified by Weyl metals having the same space group as WTe 2 ; the robustness of λ depends sensitively on the field orientation with respect to the crystallographic symmetry axis, as elaborated in Sec.V C.
For degenerate bands, the reality condition fixes P D a¼1 λ a to 0 or π, but not λ a individually.This will be exemplified by our case studies of 3D Dirac metals and Z 2 topological insulators, in Secs.V F and V G, respectively.In Sec.V D, we further demonstrate that the reality constraint may be further strengthened for spin-degenerate orbits (D ¼ 2) that are self-constrained by T symmetryinto a zero-sum rule: λ 1 þ λ 2 ¼ 0 mod 2π.Moreover, if one considers the set of fλg contributed by all closed orbits in a T-symmetric solid, we propose that fλg comprises only pairs of λ, such that each pair individually satisfies a zerosum rule, as elaborated in Sec.V E. This global constraint on fλg might be viewed as an analog of the Nielsen-Ninomiya fermion-doubling theorem [54], which states that there are always an equal number of left-and righthanded Weyl fermions on a lattice.
The four II-B classes are characterized by Let us discuss the implications of the above equation for the three following cases (i)-(iii): (i) For s ≠ u and even L (which necessarily holds if u ¼ 1), the minimal Landau-level degeneracy is L=2.As illustration, the two valley-centered Fermi surfaces (L ¼ 2) in the transition-metal dichalcogenide WSe 2 are mutually constrained by T symmetry (class II-B, u ¼ 0, s ¼ 1); the Landau levels are nondegenerate but exhibit a symmetric splitting, as elaborated in Sec.VA.
(ii) For s ≠ u and odd L, the minimal Landau-level degeneracy is L; if D is also odd, it is necessary that fλg contain either 0 or π.This is exemplified by three of four disconnected Fermi pockets (L ¼ 3) per valley of bilayer graphene (with trigonal warping) [55,56]; we refer to the three Fermi pockets that are mutually constrained by the composition of T and sixfold rotational symmetry c 6z (class II-B, u ¼ 0, s ¼ 1).For pedagogy, it is instructive to consider a model of bilayer graphene with spinless electrons; hence, the orbit associated to each Fermi pocket is nondegenerate (D ¼ 1).The threefold degenerate λ ¼ π reflects that each spinless pocket encircles a Dirac point [55,56].(iii) We would demonstrate that the same set of Fermi pockets, if equipped with multiple point-group symmetries, can belong to multiple symmetry classes in Table I.For example, the above-mentioned three pockets are also invariant under the threefold rotational symmetry c 3z (u ¼ s ¼ 0).We may therefore apply Eq. ( 13) with u ¼ s, which generally implies that the minimal Landau-level degeneracy is L. Again, L ¼ 3 in the present example.We remind the reader that this threefold degeneracy constraint was consistently implied by Tc 6z ; in addition, Tc 6z implies a stronger constraint that fλg contains either 0 or π.

V. CASE STUDIES
The utility of Table I is illustrated in the following case studies of existing conventional and topological metals, which were introduced in the previous section (Sec.IV C).

A. Orbits mutually constrained by time-reversal symmetry: Application to graphene and transition metal dichalcogenides
Our first study encompasses materials with two timereversal-related valleys in their band dispersion [57], as exemplified by graphene and monolayer WSe 2 .We will demonstrate (i) how orbits in time-reversal-invariant (T) solids can nevertheless develop a nonzero orbital magnetic moment, and (ii) the role of point-group symmetry in discretizing the Berry phase of valley-centered orbits.
To explain (i), we point out that T symmetry relates the magnetic moment of wave packets at k and −k [cf.Eq. ( 12)]; this mapping in k-space distinguishes the symmetry transformation of magnetic moments in solids from that in atoms.This allows for valley-centered orbits that are separated in k-space to individually develop a magnetic moment-since time reversal relates one valley to the other, the net magnetic moment must vanish.
In more detail, let us consider a finite chemical potential (as measured from the Dirac point) where valley-centered orbits (o i ) are disconnected; we introduce here a valley index i ∈ f1; 2g.The two orbits are mutually constrained as T∘o 1 ¼ o 2 (class II-B, u ¼ 0, s ¼ 1); each of o i is also self-constrained as Ti∘o i ¼ o i (class I, s ¼ 1).Ti imposes reality for the orbital component of wave functions at each k, leading to M ¼ 0 [cf.Eq. ( 12) with D ¼ s ¼ 1 and ǧ ¼ −I].In analyzing graphene (which has negligible spinorbit coupling), it is useful to first neglect the spinor structure of its wave function and then account for the Zeeman effect after-this was implicit in our previous assumption of D ¼ 1.Such "spinless" wave functions transform in an integer-spin representation (F ¼ 0 in Table I); the corresponding A in Table I should be interpreted as Eq. ( 2) without the Zeeman term.
The Berry phase of graphene is π, and therefore λi ≔ as is consistent with the reality constraint in (class I, s ¼ 1) of Table I.We have added an accent to λi to remind ourselves that it is the purely orbital contribution to λ i .Further accounting for the Zeeman contribution, with λi ¼ π, AE distinguishing two spin species, and m c the cyclotron mass.The symmetric splitting of λ AE about π implies an invariance under inversion: T symmetry imposes the mutual constraint: which follows from (class II-B, u ¼ 0, s ¼ 1) of Table I.Equations ( 14) and ( 15) together imply To recapitulate, we have reproduced the well-known fact that graphene's Landau levels are valley degenerate but spin split by the Zeeman effect.
The valley degeneracy of the Landau levels may be split by breaking spatial-inversion (i) symmetry, e.g., with a hexagonal BN [58] or SiC [59] substrate.In zero field, the i asymmetry is predicted to produce a band gap of 53 and 260 meV, respectively-this may be interpreted as a nonzero Semenoff mass for the Dirac fermion [60].Since each valley-centered orbit is no longer self-constrained by Ti, it develops a nonzero orbital moment (as was first noted in Ref. [61]), as well as a noninteger ϕ i B =π. Consequently, Eq. ( 14) remains valid with λi deviating from π, as we confirm with a first-principles calculation in Fig. 1(a) [62].While fλ i þ ; λ i − g (that is associated to one valley) is no longer invariant under inversion, we remark that the invariance persists for the full set: − g, owing to the mutual constraint by T symmetry [cf.Eq. ( 15)].
The predicted Semenoff masses induced by substrates tend to be small-it is instructive to compare graphene to WSe 2 , a transition metal dichalcogenide with a large Semenoff mass due to the natural absence of inversion symmetry in its space group.WSe 2 is similar to graphene in that its low-energy bands (at zero field) are also centered at two valleys, but it differs from graphene in that its bands are spin split by spin-orbit coupling.Figure 1(b) illustrates our calculated λ 1 as a function of energy, with i ¼ 1, 2 again a valley index, and λ 1 ¼ −λ 2 , owing to T symmetry; at 0.2 eV below the band maximum, λ 1 − λ 2 ≈ π=2.

B. Orbits self-constrained by mirror/glide symmetry:
Application to topological crystalline insulators λ=π ∈ Z may originate from crystalline symmetry alone.This occurs for orbits that are self-constrained by a mirror/ glide symmetry (IIA, u ¼ 1, s ¼ 0), which may arise in time-reversal-invariant and magnetic solids.Any circular orbit in (IIA, 1, 0) intersects a mirror/glide-invariant line at two points denoted by k a and k b ; at each point, the assumednondegenerate band transforms in either the even or odd representation of mirror/glide.The reality condition for e iλ in (IIA, 1, 0) of Table I has a simple interpretation: λ ¼ 0 if the representations at k a and k b are identical, and π if the two representations are distinct [47,65,66].The former is exemplified by a band that is nondegenerate at all k ⊥ bounded by o, which implies that o is contractible to a point-due to continuity of the representation along the high-symmetry line, the representations at k a and k b must be identical.λ ¼ π occurs iff there is an odd number of linear band touchings along the segment of the mirror line within o-at each touching (a Dirac point), the mirror/glide representation flips discontinuously.λ ¼ π is exemplified by the mirrorsymmetric surface states of the SnTe class [53] of topological crystalline insulators, as well as by glide-symmetric Dirac cones in band-inverted, nonsymmorphic metals [65].

C. Effect of field orientation on the crystalline symmetry of extremal orbits: Application to 3D Weyl fermion
The previous examples demonstrate that deforming the crystal structure offers a way to tune λ.For 3D solids, we may continuously tune between integer and noninteger λ=π without explicitly breaking any symmetry-by modifying the orientation of the field with respect to the crystal structure, we effectively modify the symmetry of the extremal orbit.For specific orientations, the extremal orbit may be invariant under a point-group symmetry that stabilizes integer-valued λ=π.We shall illustrate this with a 3D topological metal whose type-I Weyl points lie on a high-symmetry plane that is invariant under Ts 2z (the composition of time reversal with a twofold screw rotation), as exemplified by strained WTe 2 [67].When the field is aligned parallel to the screw axis, the maximal orbit of a Weyl point lies within the Ts 2z -invariant plane (class I, s ¼ 1) and is characterized robustly by λ ¼ π.As illustrated in Fig. 2, λ deviates from π when the field is tilted away, owing to the absence of any symmetry for the tilted maximal orbit.

D. Orbits self-constrained by time-reversal symmetry
Let us consider orbits (o) that are self-constrained as T∘o ¼ o (II-A, u ¼ 0, s ¼ 1); these are orbits that lie in T-invariant cross sections of the Brillouin torus and encircle Kramers-degenerate points.The contexts in which selfconstrained orbits will be discussed include (a) T-invariant solids with negligible SOC, (b) SO-coupled solids with both T and spatial-inversion (i) symmetries, and (c) SO-coupled solids with T but without i symmetry.Cases (a) and (b) correspond to spin-degenerate bands.In the above cases, we would show, respectively, that (a) The lack of SOC allows us to constrain the purely orbital component ( λ) of λ AE , where AE distinguishes the two spin species; recall that λ and λ AE differ only by the Zeeman splitting, as shown in Eq. ( 14).The orbitaveraged orbital moment vanishes, and λ ¼ 0 reflects the trivial Berry phase of band extrema at T-invariant wave vectors.Combining λ ¼ 0 with Eq. ( 14), we obtain the following zero-sum rule: This zero-sum rule is also satisfied for spin-degenerate orbits in spin-orbit-coupled solids: λ 1 þ λ 2 ¼ 0 mod 2π.(c) Both the orbital moment and Zeeman effect average out, and λ ¼ π reflects the nontrivial Berry phase associated to the Kramers-degenerate Dirac point at T-invariant wave vectors.

Demonstration
In all cases, the propagator A½o satisfies with TK an antiunitary representation of T (cf.row 4 of Table I with The second equation may be contrasted with the usual antiunitary representation (denoted TK) of T, which satisfies T TÃ ¼ ð−1Þ F ; the additional factor of A −1 in Eq. ( 16) indicates that Eq. ( 16) represents an extension of the magnetic point group by contractible loop translations (represented by A) in k-space [29], which generalizes a previous work on noncontractible loop translations [44].The first equation in Eq. ( 16) implies det A ¼ In case (c), we contract o T-symmetrically to an infinitesimal loop encircling the Kramers-degenerate Dirac point [69]-since the velocity is finite in this limit, the nongeometric contribution to λ vanishes; what remains is the π-Berry-phase contribution.This completes the demonstration.
One implication of (a) and (b) for spin-degenerate orbits that are self-constrained by T symmetry: the net phase offset Θ [cf.Eq. ( 10)] of two time-reversal-related fundamental harmonics can only be 0 or π.The former occurs if the two values of fλg are closer to 0 than to π and vice versa.

E. Global constraint on fλg for time-reversal-symmetric solids
Let us impose a global constraint on the set (denoted fλg) of all λ i a that are contributed by closed orbits of bulk states in a T-symmetric solid.To clarify, for a ddimensional solid, a "bulk state" is spatially extended in d directions.By combining our results for orbits that are self-constrained (cf.Sec.V D) and mutually constrained (cf.Sec.VA) by time-reversal (T) symmetry, we would show that fλg are comprised only of inversion-invariant pairs, i.e., pairs that are symmetric about zero.In our counting of "pairs," each closed orbit with energy degeneracy D contributes D values of λ, independent of whether any of these D values are mutually degenerate, or degenerate with λ from a distinct orbit.A corollary of this result is the global sum rule: P λ ¼ 0 mod 2π.Let us first prove this claim for spin-degenerate orbits.All orbits in a T-symmetric solid are either self-constrained or mutually constrained by T symmetry.As proven in Sec.V D, each self-constrained orbit contributes an inversion-invariant pair fλ; −λg.Utilizing (class II-B, u ¼ 0, s ¼ 1) of Table I, the net contribution of any T-related pair of spin-degenerate orbits is an inversion-invariant quartet fλ 1 ; λ 2 ; −λ 1 ; −λ 2 g; this was exemplified by our case study of graphene in Sec.VA.
We would next prove the global constraint for spin-orbitcoupled solids lacking spatial-inversion symmetry [case (c) in Sec.V D].By similar reasoning as in the previous paragraph, each T-related pair of spin-nondegenerate orbits contributes fλ; −λg.A self-constrained orbit necessarily contributes λ ¼ π, as proven in Sec.V D. Furthermore, self-constrained orbits always come in pairs; this follows because each self-constrained orbit encircles a T-invariant wave vector, and there are always 2 d such wave vectors in a d-dimensional lattice.This completes the proof.
It is instructive to draw analogies between our result and the Nielsen-Ninomiya fermion-doubling theorem (as applied to T-invariant solids) [54].In full generality, the theorem states that there are always an equal number of left-and right-handed Weyl fermions on a lattice; a simple implication is that Weyl fermions come in pairs.If Weyl fermions exist at wave vectors that are not T-invariant, they necessarily come in pairs due to T symmetry mapping k to −k.In the absence of point-group symmetry in the space group of the lattice, bands always come in pairs that touch at isolated T-invariant wave vectors-such Kramers-degenerate points are also Weyl points.That there are an even number of Weyl fermions on a lattice thus complements our previous observation that there are an even number of Kramers-degenerate points.If the Fermi surface encircles Kramers-degenerate Weyl points instead of intersecting them, we recover our previous claim that self-constrained orbits come in pairs.
To recapitulate, the global zero-sum rule for fλg is exemplified by all case studies of bulk orbits in this work; more generally it constrains the bulk oscillatory phenomena of all T-symmetric solids.
Let us extend our discussion to closed orbits contributed by surface states of a 3D T-symmetric solid; these orbits lie in a 2D surface Brillouin zone.By a "surface state," we mean a state localized spatially to a surface that is translation invariant in two directions; this surface lies at the interface between vacuum and a bulk that is semi-infinite in the direction orthogonal to the surface.If the solid is spinorbit coupled, surface bands are generically nondegenerate, and we would use nearly the same argument (given above) for spin-nondegenerate, bulk orbits; the sole difference is that the fermion-doubling theorem does not apply, and it is possible to have an odd number of self-constrained orbits (associated to an odd number of surface Dirac fermions).One implication is that χ s ≔ P λ (summed over surface orbits) may equal 0 or π.If the bulk of the solid is insulating, χ s may be viewed as a Z 2 index that classifies insulators from the perspective of its surface magnetotransport; this is equivalent to the Z 2 strong classification [49][50][51][52] of 3D insulators in Wigner-Dyson class AII.If the solid has negligible spin-orbit coupling, then the surface orbits satisfy the same type of global constraint as for bulk orbits-to demonstrate this, one may utilize the above argument for spin-degenerate bulk orbits.
F. Spin-degenerate orbits in spin-orbit-coupled, centrosymmetric metals: Application to the 3D Dirac metal Na 3 Bi A zero-sum rule also applies individually to each Ti-invariant orbit [o in (class I, s ¼ 1)]; this rule applies whether or not the orbit is self-constrained by T symmetry.We have previously discussed such orbits in the context of graphene; here, we extend our discussion to Ti-invariant, spin-degenerate orbits in SO-coupled solids.The reality of det A and the contractibility of o together imply Let us apply this result to a 3D massless Dirac fermiona four-band touching point between conically dispersing bands which are spin degenerate at generic wave vectors (owing to Ti symmetry).We shall assume that the 3D Dirac point is not centered at a T-invariant point.Hence, an orbit at finite chemical potential (as measured from the Dirac point) is Ti invariant but not T invariant; e.g., such Dirac points in the topological metal Na 3 Bi are stabilized by threefold rotational symmetry [70].For a field aligned along the rotational (trigonal) axis, the equidistant splitting of λ a is illustrated in Fig. 3(a) for constant-k z orbits on a surface of constant energy (0.08 eV below the Dirac point); λ 1;2 ≈ AEπ=4 for the maximal orbit (occurring at wave vector kz , as indicated by a dashed line on the plot).In Fig. 3(b), we further plot λ a ðμ; kz Þ for a range of μ below the Dirac point.The ab initio calculation of λ a is detailed in Appendix D.
Let us discuss the experimental implications for Na 3 Bi.A single Dirac point contributes to the oscillatory magnetization a term of the form of Eq. ( 4), with D ¼ 2 and λ 1;2 ðμ; kz Þ plotted in Fig. 3 We remind the reader that λ i 1 ¼ −λ i 2 owing to Ti symmetry; hence, in combination, fλ 1  1 ; λ 1 2 g ¼ fλ 2 1 ; λ 2 2 g.To reiterate, the two Dirac points contribute identically to the oscillatory magnetization, which we plot in Fig. 4 for a temperature of 0.37 K and three representative Fermi energies.These plots are obtained from Eq. ( 4), which receive as inputs the ab initio calculated band parameters fS; S zz ; m c ; λ 1;2 g for the maximal orbit, as well as the Dingle lifetime τ ¼ 10 −12 s.For the range of B plotted, kT=ℏω c ∼ 10 −2 and ω c τ ∼ 1, which implies that the Dingle factor (rather than temperature smearing) is the limiting factor in observing higher harmonics; for the parameters considered, the third harmonic does not appreciably modify the plots.The second harmonic is especially transparent at μ c ≔ −69 meV [Fig.4(b)], where λ 1 ¼ π=2 ¼ −λ 2 [cf.Fig. 3(b)] leads to the complete destructive interference of all time-reversed pairs of odd harmonics (including the fundamental harmonic) [cf.Eq. ( 9)]; i.e., the dHvA period is effectively halved.The point of destructive interference intermediates two energy intervals where the phase offset Θ [cf.Eq. ( 10)] in the fundamental harmonic differs by π: Θ ¼ π for μ > μ c and closer to the Dirac point [cf.green shaded region in Fig. 3(b)]; Θ ¼ 0 for μ ≤ μ c (blue shaded region).Our case study demonstrates that the experimental tunability [37,71,72] of μ in Na 3 Bi should expose a wide range of T-symmetric interference patterns.
For the parameter ranges in Figs. 3 and 4, the Fermi surfaces enclosing each Dirac point are indeed disconnected.As μ is increased beyond ∼10 meV, the two Fermi surfaces intersect at a saddle point and merge into a dumbbell Fermi surface.For a field parallel to the trigonal axis, this merged Fermi surface has a single minimal orbit o 3 (encircling the saddle point), in addition to the aforementioned two maximal orbits.At slightly higher energy (approximately 20 meV), the dumbbell transforms into an ellipsoid, and correspondingly, the three extremal orbits merge into a single maximal orbit (denoted o 4 ).o 3 (and also o 4 ) is T invariant [under (class II-A, u ¼ 0, s ¼ 1)] and Ti invariant (class I, s ¼ 1); this exemplifies how a single orbit having multiple symmetries may fall under different classes in Table I.Both classes impose (consistently) a zero-sum-rule constraint on fλg.

Comment on magnetotransport experiments of 3D Dirac metals
We conclude this section by commenting on the interpretation of several magnetotransport experiments on 3D Dirac metals [73][74][75][76][77][78].The experimental data in these works have only been fitted to the fundamental harmonic (r ¼ 1) in the traditional Lifshitz-Kosevich formula [6].The phase offset in this fitted fundamental harmonic is interpreted as ϕ B − ϕ M AE π=4; ϕ B ¼ π that is inferred from such a fitting is viewed as a smoking gun for a 3D Dirac metal.These interpretations have been justified [73][74][75][76] by appealing to the one-band Onsager-Lifshitz quantization rule [Eq.( 1) We propose an alternative interpretation that nevertheless makes sense of the robustness of their measurements, across different experimental groups and slightly varying samples-this robustness is solely a consequence of T symmetry, rather than anything intrinsic to a 3D-Dirac Fermi surface.In other words, Θ ¼ π should not be viewed as a smoking gun for a 3D Dirac metal.
As described in Eqs. ( 4)- (10), the phase offset (Θ) of the summed fundamental harmonic, after subtracting the Maslov and Lifshitz-Kosevich corrections, is a suitably defined average of λ 1;2 .As proven in Sec.V D, Θ is restricted to 0 or π for orbits that are self-constrained by T symmetry.We believe that previous measurements [73][74][75][76][77][78] of Θ have been misinterpreted as a measurement of the Berry phase.The convergence in measurements of several groups indicate the quantization of Θ due to T symmetry.However, Θ ¼ π is not an intrinsic property of 3D Dirac metals, as we have exemplified by our case study of Na 3 Bi.Moreover, we point out that Θ ¼ π has been measured in conventional (i.e., nontopological) metals, even in the first metal (3D, single-crystal bismuth) for which quantum oscillations were reported [4,5].As far back as 1954, it was remarked by Dhillon and Shoenberg that "For Bismuth good agreement was found by the theoretical formula except that the signs of the fundamental and odd harmonics had to be reversed" [41].Here, the "theoretical formula" refers to the traditional Lifshitz-Kosevich formula [3] without the Θ correction to the phase.In a later work, the quantization of Θ was implicit in the form of Shoenberg's spin "reduction factor" [6,34], which involves a phenomenological g-factor in the presence of spin-orbit coupling.
One contribution of this work is to present an analytic formula for the effective g-factor through the spectrum of Eq. ( 2) [80].Our proposed formula is calculable from ab initio wave functions (as elaborated in Appendix D), which allows for a quantitative comparison between theory and experiments.Another contribution of this work is the recognition that the quantization of Θ originates from T symmetry; hence, in magnetic materials, Shoenberg's spin reduction factor is not generally applicable, but our generalized Lifshitz-Kosevich formulas in Sec.III remain valid.
We further emphasize that Θ is an incomplete characterization of the Fermi-surface orbit; a more complete characterization involves measuring the individual values of λ 1;2 , which carry the complete non-Abelian information about the Fermi-surface wave function.Two methods are available to measure λ 1;2 : (i) Where either of ω c τ ≪ 1 or kT ≫ ℏω c applies, one needs to measure both the amplitude and the phase offset of the fundamental harmonic; λ 1;2 is then extracted from the r ¼ 1 terms in Eq. ( 4) and Eqs. ( 9)- (10) [82].(ii) Where neither of the above inequalities apply, the experimental data should be fitted to the full harmonic expansion of Eq. ( 4) and λ 1;2 extracted accordingly.Note that the zero-sum rule for λ 1;2 implies that there is only one parameter to be fitted.
Going beyond general symmetry constraints, one may ask if Θ and λ 1;2 may be analytically calculated in a caseby-case basis.Given material-specific information about the magnetic point group of the 3D Dirac point, as well as the relevant symmetry representations that span the lowenergy Hilbert space at the Dirac point, one may indeed determine Θ from a k • p analysis, as exemplified in the next section.
G. Calculating λ 1;2 from k • p analysis: Application to 3D massive Dirac fermions in Z 2 topological insulator Bi 2 Se 3 In certain symmetry classes, and for sufficiently small (but nonzero) chemical potentials (μ) relative to the Dirac point, 3D massless Dirac fermions may exhibit λ 1 ≈ λ 2 ≈ π; this would imply from Eq. ( 10) that Θ ¼ π.Precisely, we restrict jμj to be small enough that a linearized k • p Hamiltonian is a good approximation-how small jμj needs to be to fulfill the above criterion depends on materialspecific band parameters.As far as the linear approximation is valid, λ 1;2 that follows from the subsequent k • p analysis is approximately π (hence, our use of ≈ throughout this section).λ 1;2 should manifest in the dHvA oscillations for a sufficiently weak field, such that the semiclassical approximation remains valid, i.e., l 2 SðμÞ ≫ 1 [83].
Let us demonstrate how λ 1;2 ≈ π arises for 3D massless Dirac fermions centered at a wave vector (Γ) with a magnetic point group (P) that combines T symmetry with the D 3d point group [84].λ 1;2 ≈ π relies not just on P, but also on symmetries that are absent in P; these additional symmetries emerge only at long wavelength (k → 0), where the averaged-out crystal field has higher symmetry than the magnetic point group.
Our study of D 3d is motivated by the identical symmetry class of Bi 2 Se 3 , a well-known Z 2 topological insulator [49][50][51][52].For the symmetry representations of Bi 2 Se 3 , the k • p Hamiltonian around Γ (of the 3D Brillouin torus) assumes the form [85,86] to quadratic accuracy in k.Here, σ j and τ j are Pauli matrices spanning spin and orbital (labeled P1 þ z and P2 − z ) [85] subspaces, respectively, and , with Δ, the Dirac mass, vanishing at the critical point between trivial and topological insulator.
Assuming Δ ¼ 0 for now, ½H 0 ðkÞ; I 2×2 ⊗ τ 1 ¼ Oðk 2 Þ is an emergent conservation law of the linearized H 0 .Blockdiagonalizing H 0 with respect to I 2×2 ⊗ τ 1 ¼ AE1, we derive two decoupled 3D Weyl Hamiltonians H AE with opposite chirality, each satisfying the time-reversal constraint TH AE ðkÞT −1 ¼ H AE ð−kÞ, with T ¼ iσ 2 K squaring to −I in a half-integer-spin representation (F ¼ 1).Independent of the field orientation, the spin-degenerate extremal orbit effectively decouples to two nondegenerate (D ¼ 1) extremal orbits that are individually T invariantthis implies from a previous demonstration [case (c) of Sec.V D] that λ 1;2 ≈ π, which may be viewed as the Berry phase (ϕ B ) for each decoupled orbit.This would effectively imply a single set of dHvA harmonics [indexed by r in Eq. ( 4)] with twice the usual amplitude and a Lifshitz-Kosevich phase correction of −π=4 for a maximal orbit.
Away from the critical point, we prove in Appendix C that the degeneracy is Zeeman split as λ 1;2 ≈ π AE πΔ=mv 2 , with m the free-electron mass, for a field parallel to the trigonal axis (k z ).It is worth remarking that ϕ B reduces from π, and the phase (ϕ R ) associated to the orbital moment increases from 0; however, their sum (ϕ B þ ϕ R ) remains fixed to π [87].
In fact, naturally occurring Bi 2 Se 3 does not lie at the critical point; i.e., the corresponding 3D Dirac fermion is massive: we estimate Δ=mv 2 ≈ 0.13 utilizing fitted parameters from an ab initio calculation in Ref. [85].Utilizing Eq. ( 10) with jΔ=mv 2 j < 0.5, we further derive Θ ¼ π.This splitting of λ manifests as two sets of harmonics in the dHvA oscillation, which should be measurable utilizing techniques outlined in Sec.V F. The tunability of Δ by doping [88] or strain [89][90][91] suggests that the dHvA oscillations are correspondingly tunable-in particular, jΔ=mv 2 j ¼ 0.5 is the point of complete destructive interference for the fundamental harmonic; this point separates two regimes where Θ ¼ 0 vs π.
There exists galvanomagnetic evidence that supports the quantization of Θ.The fundamental Shubnikov-de Haas harmonic [5,92,93] in the transverse resistivity of Bi 2 Se 3 has been fitted [94] as with the higher harmonics suppressed by temperature.The above proportionality may be derived from (a) ρ xx ≈ σ xx =σ 2 xy , which is valid for a large Hall angle σ xy ≫ σ xx , and (b) the proportionality between the oscillatory components of the transverse conductivity and the longitudinal magnetic susceptibility: Δσ xx ∝ m c Δχ [95,96].The cyclotron mass (m c ) is positive for the lone electronlike orbit in Bi 2 Se 3 [94], and the oscillatory susceptibility (Δχ) is obtained from differentiating Eq. ( 4) (with D ¼ 2) with respect to B and keeping only the fastest oscillatory terms.The proportionality stated in (b) is valid for metals having only a single extremal orbit, and for sufficiently weak fields such that μ ≫ ℏω c [97].The fit in Eq. ( 18) implies Θ ¼ 0.4-1=2 þ 1=8 ≈ 0 to the accuracy of their fit.It was further remarked in Ref. [94], without explanation, that the measured Θ ≈ 0 was independent of the field orientation.To elaborate, four measurements, carried out at angles 0, π=4, π=3, and π=2 relative to the trigonal axis in Bi 2 Se 3 , produced the same value for Θ.Our explanation for this robustness is that, independent of the field orientation, the extremal orbit is invariant under timereversal symmetry, which quantizes Θ to 0 or π.In contrast, crystalline symmetries of an extremal orbit depend sensitively on the field orientation, as explained in Sec.V C.

VI. DISCUSSION
In fermiology, the shape of the Fermi surface is deducible from the period of dHvA oscillations; here, we propose that the topology of the Fermi-surface wave functions is deducible from the phase offset (λ a ) of dHvA oscillations in 3D solids, as well as in fixed-bias oscillations of the differential conductance in scanningtunneling microscopy.λ a manifests as the complete, subleading [Oð1Þ] corrections to the Bohr-Sommerfeld quantization rules for closed orbits without breakdown [cf.Eq. ( 1)]; we have formulated λ a as eigenphases of a propagator A defined in Eq. (2).
In certain solids and for certain field orientations with respect to a crystal axis, P D a¼1 λ a (with D the degeneracy of the low-energy band subspace) are topologically invariant under deformations of the zero-field Hamiltonian that preserve the symmetry and global shape of the orbit.Precisely, globally equivalent orbits correspond to a graph with a homotopy equivalence defined in Ref. [[29]].To identify orbits with robustly integer-valued P D a¼1 λ a =π, as well as identify the degeneracy of Landau levels, we classify symmetric orbits into ten (and only ten) classes summarized in the first three columns of Table I; the rest of the table describes the corresponding constraints on the propagator A.
The results of this table remain valid if we substitute A → W and λ a → ϕ B;a ; here, W is defined as the unitary matrix generated by the Berry connection and may be viewed as the purely geometric component of A; ϕ B;a are defined as the eigenphases of W and may be viewed as the non-Abelian generalization [26] of the Berry phase.W provides a purely geometric characterization of bands that is intimately related to the topology of wave functions over the Brillouin torus [44,98,99].While W (ϕ B;a ) generically differs from A (λ a ) due to the orbital moment and the Zeeman effect, W and A transform identically under symmetry, and therefore, they satisfy the same constraints.In particular, if we define fϕ B g as the set of ϕ B contributed by all bulk orbits in a T-symmetric solid, then fϕ B g comprises only pairs of ϕ B that individually satisfy a zero-sum rule, in close analogy with the global constraint on fλg that is described in Sec.V E.
Higher-order (in jBj) corrections to the Bohr-Sommerfeld quantization rule become relevant in higherfield experiments that intermediate the semiclassical and quantum limits [16,17]; these corrections may be interpreted as zero-field magnetic response functions [18].These corrections are accounted for in the generalized Lifshitz-Kosevich formulas [cf.Eqs. ( 4)- (8)] by the simple replacement l 2 SðμÞ þ λ a ðμÞ → l 2 SðμÞ þ λ a ðμ; jBjÞ, with λ a expandable asymptotically in powers of jBj; λ a ðμ; jBj ¼ 0Þ ¼ λ a ðμÞ is obtained from diagonalizing the propagator A [cf. Eq. ( 2)].In addition, the field dependence of the chemical potential μ due to fixed particle density [6] cannot be neglected at the order of accuracy in which we consider the field dependence of λ.It would be interesting to extend our symmetry analysis to the field-dependent component of the phase offset.The linearized k • p Hamiltonian from Eq. ( 17) is with a nonzero mass Δ.This Hamiltonian has an emergent reflection symmetry: which is broken by terms that are cubic in k [86].If the field is oriented parallel to the trigonal axis ðk z Þ, the maximal orbit lies at k z ¼ 0, owing to the just-mentioned symmetry.
It is here that ½Hðk ⊥ ; 0Þ; ȓz ¼ 0, and we may blockdiagonalize H with respect to even and odd representations of ȓz : H AE describe two 2D massive Dirac fermions with opposite chirality.It is known for each of H AE that ϕ R þ ϕ B ¼ π is independent of the Semenoff mass (Δ) [63]; we rederive this result in our language in Appendix C, where we further clarify the applicability of their result to more general twoband Hamiltonians.What remains is to calculate the Zeeman contribution to λ. Employing that the spin operator S z is represented ℏγ z =2, and the Zeeman phase is Employing the parameters from Ref. [85], Δ ¼ 0.28 eV, v ¼ 6.2 × 10 5 ms −1 , we obtain ϕ Z ≈ 0.13π.

APPENDIX C: REALITY OF THE ORBITAL COMPONENT OF e iλ
For a general dispersion ϵðkÞ, the area of constant energy contour is a function of ϵ: SðϵÞ.For a general scalar-valued function fðkÞ, we define a FðϵÞ ≔ R SðϵÞ fðkÞd 2 k.Then the average of f on the contour is defined by where M is orbital magnetization.The sign change comes from a change in orientation in path.For the two-band Hamiltonian in the form of energy dispersion is symmetric with respect to zero.Here, r is real and Δ is k independent.In this case, orbital magnetization is simply related to Berry curvature by M ¼ ½e=ðℏcÞΩϵ, where Berry curvature is defined by ½ðdϵϕ B Þ=ðdϵÞ is evaluated to be AEπ R ↻ dk • ∇ k θ for the conduction band and valence band, respectively, which is fixed to 0 or π (mod 2π).
Any two-band Hamiltonian that can be transformed into Eq.(C3) (i.e., with at least one k-independent multiplicative coefficient of a Pauli matrix) is also characterized by real e iðϕ R þϕ B Þ .We propose a sufficient condition for the existence of this transformation, which we denote by U in what follows.
a ¼ AE1, with D ¼ 2 and 1 for cases (b) and (c), respectively; in case (a), we should interpret det A ¼ e i λ ¼ AE1.For cases (a) and (b), we may further restrict det A to þ1 by the following argument: while preserving det A, contract o T-symmetrically to an infinitesimal loop that encircles the T-invariant point.In the classes of (a) and (b), the band dispersion is extremized at T-invariant points; hence the band velocity vðkÞ vanishes along the infinitesimal loop [68].det A is thus solely contributed by the nongeometric one-forms that depend inversely on the velocity.Further applying the identity det exp R VðkÞjdkj ¼ exp R TrVðkÞjdkj for any matrix VðkÞ and the time-reversal constraints TrMðkÞ ¼ −TrMð−kÞ [cf.Eq. (12) with g ¼ T] and TrσðkÞ ¼ −Trσð−kÞ, we derive the desired result.The above proof required path-ordering and matrix traces in case (b), where A, M, and σ are 2 × 2 matrices; such matrix operations are not needed in case (a).

FIG. 2 .
FIG. 2. (a) Fermi surface of a Weyl fermion centered on a generic wave vector in a high-symmetry plane; red circle indicates an extremal orbit that has no symmetry.(b) For a model of the Weyl fermion, we plot λ as a function of the field orientation, which we parametrize by spherical coordinates illustrated in (a).
(b).Na 3 Bi actually has two Trelated Dirac points, which both lie on the rotationinvariant line through Γ [70].Assuming the Fermi surfaces of the two Dirac points are disconnected, the effect of the second Dirac point is simply to double the magnitude of the oscillations.This follows because the maximal orbits (o 1 , o 2 ) on disconnected Fermi surfaces are mutually constrained by T symmetry; hence, from (class II-B,

FIG. 3 .
FIG. 3. Characterization of the 3D Dirac metal Na 3 Bi.(a) Left: Surface of fixed energy ( Ē ¼ −0.08 eV measured from the 3D Dirac point); contours of fixed energy and k z are illustrated by red circles.Right: Plot of fλ a ð Ē; k z Þg a¼1;2 for a field parallel to the trigonal axis (k z ).k z is measured in units where the reciprocal period (0.64 Å −1 ) is unity; the minimum (maximum) k z on the plot corresponds to the south (north) poles of the constant-energy surface.The extremal orbit occurs at k z ¼ kz ð ĒÞ, indicated by a dashed line.(b) λ 1;2 ½E; kz ðEÞ vs E, for a range of E below the Dirac point.
oscillatory terms, i.e., only the terms derived from D acting on the cosine function, πD sin πD cosðΓÞ ¼ πkTdΓ=dμ sinhðπkTdΓ=dμÞ cosðΓÞ; ðA8Þ we arrive at the desired result.APPENDIX B: 3D MASSIVE DIRAC FERMION WITH D 3D SYMMETRY