RRevealing the topology of Fermi-surface wavefunctions from magnetic quantumoscillations
A. Alexandradinata, Chong Wang, Wenhui Duan,
2, 3, 4 and Leonid Glazman Department of Physics, Yale University, New Haven, Connecticut 06520, USA Institute for Advanced Study, Tsinghua University, Beijing 100084, China Department of Physics and State Key Laboratory of Low-DimensionalQuantum Physics, Tsinghua University, Beijing 100084, China Collaborative Innovation Center of Quantum Matter, Tsinghua University, Beijing 100084, China (Dated: January 10, 2018)The modern semiclassical theory of a Bloch electron in a magnetic field now encompasses theorbital magnetic moment and the geometric phase. These two notions are encoded in the Bohr-Sommerfeld quantization condition as a phase ( λ ) that is subleading in powers of the field; λ is measurable in the phase offset of the de Haas-van Alphen oscillation, as well as of fixed-biasoscillations of the differential conductance in tunneling spectroscopy. In some solids and for certainfield orientations, λ/π are robustly integer-valued owing to the symmetry of the extremal orbit,i.e., they are the topological invariants of magnetotransport. Our comprehensive symmetry analysisidentifies solids in any (magnetic) space group for which λ is a topological invariant, as well asidentifies the symmetry-enforced degeneracy of Landau levels. The analysis is simplified by ourformulation of ten (and only ten) symmetry classes for closed, Fermi-surface orbits. Case studiesare discussed for graphene, transition metal dichalchogenides, 3D Weyl and Dirac metals, andcrystalline and Z topological insulators. In particular, we point out that a π phase offset in thefundamental oscillation should not be viewed as a smoking gun for a 3D Dirac metal. I. INTRODUCTION
The semiclassical Peierls-Onsager-Lifshitz theory connects experimentally-accessible quantities in mag-netic phenomena to Fermi-surface parameters of the solidat zero field. For example, field-induced quantum oscil-lations of magnetization and resistivity have becomethe leading method to map out the shape of the Fermisurface of normal metals and superconductors – thisphenomenology has been coined ‘Fermiology’. The semiclassical theory has been extended to in-corporate two modern notions: a wavepacket orbitingin quasimomentum ( k ) space acquires a geometric phase( φ B ), as well as a second phase ( φ R ) originating fromthe orbital magnetic moment of a wavepacket around itscenter of mass. Further accounting for the well-knownZeeman coupling, λ := φ B + φ R + φ Z is known to bethe complete, subleading (in powers of the field) correc-tion to the Bohr-Sommerfeld quantization rule for nonde-generate bands. λ is measurable as a phase offset inoscillations of the magnetization/resistivity in 3D solids,as well as in fixed-bias oscillations of the differential con-ductance in tunneling spectroscopy. While it is conventionally believed that φ B =0 vs π dis-tinguishes between Schr¨odinger and Dirac systems, wepropose to view φ B /π as a continuous quantity that issometimes fixed to an integer in certain space groupsand for certain types of field-dependent orbits; moreover,while φ R vanishes for centrosymmetric metals withoutspin-orbit coupling (SOC), it plays an oft-ignored rolein most other space groups. Our comprehensive sym-metry analysis identifies the (magnetic) space groups inwhich λ/π is robustly integer-valued – we will formulate λ as a topological invariant in magnetotransport, which is distinct from the traditional formulation of topolog-ical invariance in band insulators. We also extendour symmetry analysis to the multi-band generalizationof λ , with envisioned application to bands of arbitrarydegeneracy ( D ); D =2 is exemplified by spin degeneracy.Let us outline the organization of the main text. Webegin in Sec. II by introducing the multi-band quantiza-tion rule, and describing how λ appears as the sublead-ing phase correction. Experimental methods to extract λ are discussed in Sec. III; we present here generalizedLifshitz-Kosevich formulae for the oscillatory magnetiza-tion and density of states. These formulae extend previ-ous works in their applicability to orbits of any energydegeneracy and symmetry, including orbits in magneticsolids. In Sec. IV, we provide a general group-theoreticframework to identify solids for which λ takes only dis-crete values. In addition, our symmetry analysis identi-fies the symmetry-enforced degeneracy of Landau levels;where degeneracy is not enforced, symmetry may never-theless constrain the possible splittings of Landau levels.We exemplify our symmetry analysis with several casestudies in Sec. V, including: graphene, transition metaldichalchogenides, surface states of topological insulators,and 3D Weyl and Dirac metals. In particular, we pointout that a π phase offset in the fundamental oscillationshould not be viewed as a smoking gun for a 3D Diracmetal. We recapitulate our main results in the conclud-ing Sec. VI; a final remark broadens the applicabilityof our symmetry analysis to matrix representations ofholonomy in the Brillouin torus, also known as Wilsonloops of the Berry gauge field. a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n II. MULTI-BAND BOHR-SOMMERFELDQUANTIZATION RULE
The quantization rule is derived from the effectiveHamiltonian ( H ) that describes the low-energy dynamicsof Bloch electrons in a field. In a basis compris-ing D field-modified Bloch functions at each wavevector, H ( K ) is a Weyl-symmetrized, matrix function of the ki-netic quasimomentum operators, whose noncommutivityis manifest in K × K = − ie B /c . H is asympotically ex-pandable in powers of the field: H + H + . . . , where theleading-order term is the Peierls-Onsager Hamiltonian, while the subleading terms H := H B + H R + H Z respec-tively encode the geometric phase, the orbital moment,and the Zeeman effect. In the WKB approximation, the D -component vec-tor wavefunction of H generalizes the known single-component solution developed by Zilberman andFischbeck for a nondegenerate band. In the absenceof breakdown, continuity of the vector wavefunc-tion around a closed orbit ( o ) affords us the followingquantization rule: l S ( E, k z ) + λ a ( E, k z ) = 2 πj + φ M . (1)We have assumed here that the field is oriented in (cid:126)z ,such that o is a contour of the band dispersion at fixedenergy E and k z . By ‘closed orbit’, we mean that o does not wrap around the Brillouin torus. o bounds aregion in k ⊥ :=( k x , k y )-space with positive-definite area S . l :=( (cid:126) c/e | B | ) / above is the magnetic length, j aninteger, a ∈ Z D := { , , . . . , D − } . The Maslov correc-tion ( φ M ) depends on the topology of the Fermi-surfaceorbit, e.g., it equals π for orbits that are deformable toa circle, but vanishes for a figure-of-eight orbit. Toleading order in l -2 , Eq. (1) without λ a is a well-knownresult by Lifshitz and Onsager. The remaining term λ a is defined through the spec-trum ( { e iλa } Da =1 ) of the propagator ( A ) that is generatedby H over the cyclotron period. A may be expressed asa path-ordered exponential (denoted exp): A [ o ] = exp (cid:104) i (cid:72) o (cid:110) ( X + A ) · d k + g (cid:126) mv ⊥ σ z | d k | (cid:111)(cid:105) , (2)with o carrying a clockwise orientation; the above threeone-forms represent contributions by H B , H R and H Z respectively. The first one-form is the non-abelianBerry connection for the D -fold-degenerate sub-space (henceforth denoted by P ), and is defined by X ( k ) mn := i (cid:10) u m k (cid:12)(cid:12) ∇ k u n k (cid:11) , with e i k · r u n k the Bloch func-tion of a band labelled by n ∈ Z D . The multi-band orbitalmagnetic moment is encoded in the second one-form A mn · d k = (cid:88) l/ ∈ Z D X xml Π yln dk x / v y + ( x ↔ y ) , (3)with m, n ∈ Z D . Π ( k ) ln = i (cid:10) u l k (cid:12)(cid:12) e − i k · ˆ r [ ˆ H , ˆ r ] e i k · ˆ r (cid:12)(cid:12) u n k (cid:11) / (cid:126) are matrix elements of the velocity operator, with ˆ H thesingle-particle, translation-invariant Hamiltonian and ˆ r the position operator. For n ∈ Z D , Π nn = v is the veloc-ity of each band in P , and v ⊥ :=( v x + v y ) / . X and Π in Eq. (3) comprise only off-block-diagonal matrix ele-ments between P and its orthogonal complement. Thethird one-form in Eq. (2) is the well-known Zeeman cou-pling, with (cid:126) σ zmn ( k ) / S z , g ≈ m thefree-electron mass. While the definition of A presumesa basis choice { u n k } Dn =1 within P , one may verify thatthe eigenvalues of A are independent of this choice. Bysetting D to 1 in the above equations, the line integral ofthe three one-forms in Eq. (2) give λ := λ = φ B + φ R + φ Z respectively, as we have introduced in the second para-graph of the paper.When Eq. (1) is viewed at fixed field, the discrete ener-getic solutions { E a,j } correspond to D sets of sub-Landaulevels; within each set labelled by a , the difference be-tween two adjacent levels ( | E a,j +1 − E a,j | ) is approxi-mately (cid:126) ω c :=2 π/ ( l | ∂S/∂E | ) evaluated at E a,j . WhenEq. (1) is viewed at constant energy (e.g., the chemicalpotential µ ), the discrete solutions ( { l a,j } ) correspond tovalues of the field where Landau levels successively be-come equal to µ . In thermodynamic equilibrium, theseare also the fields where Landau levels are suddenly de-populated with a periodicity: l a,j +1 − l a,j =2 π/S ( µ ), foreach of a ∈{ , . . . , D } . This results in various oscillatoryphenomena from which we may extract λ a . III. GENERALIZED LIFSHITZ-KOSEVICHFORMULAE TO EXTRACT λ In the de Haas-van Alphen effect, each extremal orbit( o ) on the Fermi surface of a 3D metal is associated toan oscillatory contribution to the longitudinal magneti-zation (parallel to the field in (cid:126)z ): δ M = − π ) / kT | B | Sl | S zz | / × D (cid:88) a =1 ∞ (cid:88) r =1 e − rπω c τ sin (cid:2) r (cid:0) l S + λ a − φ M (cid:1) ± π/ (cid:3) r / sinh (2 π rkT / (cid:126) ω c ) , (4)which is a sum of D sets of harmonics. Being valid inthe degenerate ( µ (cid:29) kT ) and semiclassical ( µ (cid:29) (cid:126) ω c ) lim-its, Eq. (4) is our generalization of the Lifshitz-Kosevichformula to orbits of any energy degeneracy ( D ) andsymmetry, including orbits in magnetic solids. In com-parison, the commonly-employed Lifshitz-Kosevich for-mula with a ‘spin reduction factor’ is only applicableto two-fold degenerate orbits in solids with both time-reversal and spatial-inversion symmetries. All quanti-ties on the right-hand side of Eq. (4) are evaluated on o , which may be electron- or hole-like; the sign of π/ S zz is thedouble derivative of S with respect to k z , and we have in-troduced Dingle’s damping factor that depends on thequasiparticle’s mean free time ( τ ).For D =1, the field-independent phase in the argumentof the fundamental harmonic is sometimes referred to asthe Onsager phase: − πγ := λ − φ M ± π/ . (5)If a Fermi surface has multiple extremal orbits, eachextremal orbit additively contributes a term with thesame functional form as Eq. (4). If two extremal or-bits ( o i and o i +1 ) are symmetry-related, they contributeoscillatory terms that are identical in the parameters { S, S zz , ω c , φ M } , but not necessarily for the λ -phase cor-rections. Generally, { λ ia } Da =1 = {± λ i +1 a } Da =1 (defined mod-ulo 2 π ), with the sign depending on the symmetry classof the orbit, as we will elaborate in Eq. (13) of Sec. IV.In 2D metals, the analogous oscillatory formula is δ M = − π kT | B | S D (cid:88) a =1 ∞ (cid:88) r =1 e − rπω c τ sin[ r ( l S + λ a − φ M )]sinh (2 π rkT / (cid:126) ω c ) (cid:12)(cid:12)(cid:12)(cid:12) µ . (6)The field dependence of µ is negligible in the semiclassi-cal limit for 3D metals, as well as for 2D surface statesof 3D solids. In strictly-2D metals with a fixed particledensity, field-induced oscillations in µ render the extrac-tion of λ implausible.An alternative method to extract λ a is to measure thetemperature-broadened, 3D density of states, defined as G T ( E + µ ) := − (cid:90) ∞−∞ dε f (cid:48) T ( ε − µ − E ) g ( ε ) . (7)Here, g ( E ) is the density of states of 3D Landau levels(smoothened by the Dingle factor); f (cid:48) T ( x ) is the deriva-tive of the Fermi-Dirac distribution function, and ap-proaches − δ ( x ) in the zero-temperature limit. The oscil-latory component of G may be expressed as a harmonicexpansion: δ G ( E ) = √ π kT ( (cid:126) ω c ) l | S zz | / × D (cid:88) a =1 ∞ (cid:88) r =1 r / e − rπω c τ cos (cid:2) r ( l S + λ a − φ M ) ± π/ (cid:3) sinh (2 π rkT / (cid:126) ω c ) (cid:12)(cid:12)(cid:12)(cid:12) E, ¯ k z . (8)This quantity may be accessed via the scanning tun-neling microscope (STM) or by planar tunnelingjunctions. The oscillatory component of the differ-ential conductance ( dI/dV ), averaged over the surfaceof a 3D metal at fixed bias voltage ( V ) for the STM,is directly proportional to δ G ( µ + eV ) in the absence ofsurface-localized states; for the planar tunneling junc-tion, no spatial averaging is needed for a sufficiently largejunction size. The just-mentioned proportionality pre-supposes that: (a) the tunneling matrix elements and thedensity of states of the STM tip are featureless in the en-ergy range that is accessed by the bias voltage, and that(b) the tip and sample have thermally equilibrated. Landau-level spectroscopy via the scanning tunneling mi-croscope has already been reported, but we are notaware that the phase offset of the oscillations has everbeen measured. Further details on the derivation of Eqs.(4-8) are provided in App. A.Let us discuss how to extract λ from dHvA data. Forsimplicity in presentation, we consider the magnetiza-tion oscillations contributed by a single orbit (extremalorbit in 3D metals). We assume that l S , (cid:126) ω c and theDingle lifetime τ have already been extracted by stan-dard techniques. Let us first consider either ω c τ (cid:28) kT (cid:29) (cid:126) ω c , such that the harmonic expansion is domi-nated by the fundamental ( r =1) harmonic. If D =1, thenthe experimental data may directly be fitted to a singlesine function offset by the Onsager phase [cf. Eq. (5)].If D =2 (e.g., spin degeneracy), the sum of two funda-mental harmonics produces an equi-frequency harmonicproportional to (cid:12)(cid:12)(cid:12)(cid:12) cos (cid:18) λ − λ (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) sin (cid:16) l S +Θ − φ M ± π (cid:17) (9)withΘ := λ + λ π (cid:18) − sign (cid:20) cos (cid:18) λ − λ (cid:19)(cid:21)(cid:19) (10)defined to be invariant (modulo 2 π ) under λ j → λ j +2 π ;these formulae will be applied to a case study of Bi Se in Sec. V G. If | S zz | is not otherwise measurable (for 3Dmetals), it is not possible to fully determine λ , owingto our ignorance of the amplitude of the fundamentalharmonic; measuring the phase offset of the fundamentalharmonic merely determines Θ, a suitably-defined aver-age of λ , . Θ alone does not completely characterize thenon-abelian transport within the two-band subspace.In the interest of measuring individual values of λ , ,we propose higher-field dHvA measurements of cleanersamples ( ω c τ not (cid:28)
1) and at lower temperatures ( kT not (cid:29) (cid:126) ω c ). In this regime, not just the fundamentalbut also higher ( r>
1) harmonics are needed to accuratelyrepresent the dHvA data. Since { λ a } Da =1 is encoded inthe interference of multiple harmonics, { λ a } Da =1 may beextracted without knowledge of the absolute amplitudeof a single harmonic.In metals with multiple extremal orbits, there mayexist field orientations where all orbits are related bysymmetry, which simplifies the fitting to the Lifshitz-Kosevich formula. Independent of the field orienta-tion, there exists one simplification for any non-magneticmetal: due to time-reversal ( T ) symmetry, the set of all { λ } comprises only pairs that are invariant under inver-sion about zero ( λ →− λ ), which effectively halves the in-dependent parameters that require fitting. We quote thisresult here to exemplify the utility of a symmetry analy-sis of { λ } , which we explore in greater generality in thenext 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 (cid:80) Da =1 λ a /π for D>
1] is integer-valued owing to the symmetry of the ex-tremal orbit. To identify these orbits and space groups,it is useful to distinguish ten symmetry classes for closedorbits – to each class we associate certain constraints forthe propagator A and its spectrum [ { λ a } Da =1 ], as summa-rized in Tab. I. The goal of this section is identify therelevant symmetry class of closed orbit – for any physi-cal system one chooses to study. Once this identificationis made, the resultant symmetry constraints on { λ a } Da =1 may be read off from the last column of Tab. I, and ver-ified experimentally from any of the methods detailed inSec. III. This will be exemplified for several case studiesin Sec. V.Let us first restate the problem in simple terms: the dy-namics of Bloch electrons immersed in (cid:126)z are restricted toBrillouin two-tori ( BT ⊥ ) of fixed k z . For a D -fold degen-erate band subspace with dispersion ε ( k ), semiclassicalmotion occurs along (assumed) closed orbits defined by ε ( k ⊥ , k z )= E , with k ⊥ parametrizing BT ⊥ ( k z ). If multi-ple disconnected orbits exist within the same BT ⊥ , weassume they are sufficiently separated in k ⊥ -space thattunneling is negligible. Neglecting the subleading term( H ) in the effective Hamiltonian, all Landau levels areat least D -fold degenerate owing to the Onsager-Lifshitzquantization rule; here and henceforth, the ‘degeneracyof a Landau level’ is defined in units where the exten-sive degeneracy of a Landau level (associated to a singlespinless orbit) equals 1 [cf. Eq. (A3)]. For a subset ofLandau levels, this zeroth-order degeneracy is enhancedto LD if a symmetry ( g ) constrains L disconnected or-bits to have identical shape. We may ask if (and how) H splits this LD -fold (or D -fold) degeneracy; if L = D =1,we ask if H shifts the zeroth-order Landau-level spec-trum at all. The answer to these questions depends onthe class of symmetric orbit, which we proceed to analyzein full generality. A. Tenfold classification of symmetric orbits
Since lattice translations trivially constrain A , we shallhenceforth focus on symmetries ( g ) of the solid that cor-respond to nontrivial elements in the (magnetic) pointgroup ( P ) of the solid; examples include (screw) rota-tions, (glide) reflections, spatial inversion and time rever-sal. Technically speaking, magnetic point groups differfrom point groups in that the former includes a symmetrythat reverses time; however this distinction is irrelevantto the following classification of symmetric orbits, hencewe hereafter use ‘point group’ democratically.We are interested only in g that maps BT ⊥ to itself;such g correspond to a subgroup ( P ⊥ ) of P that gener-ally depends on the field orientation as well as k z . Anyconfiguration of closed orbits in BT ⊥ may be divided intoa disjoint set of elementary orbits { ( g, O i ) } , where O i is defined to be the smallest, closed orbit configuration thatis invariant under g . By ‘invariance’, we mean that forevery k ⊥ ∈O i , the map of k ⊥ under g (denoted as g ◦ k ⊥ )belongs also in O i . Similarly, if a closed orbit o ∈O i , sowould g ◦ o ∈O i .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 in-dividually invariant under g . Such mappings are labelledas class I, and all other mappings are of class II. We fur-ther distinguish between class-II mappings where g ◦ o isidentical to o up to orientation [class II-A], or they aredisconnected orbits [class II-B].There are two classes of class-I elementary orbits dis-tinguished by whether g is purely a spatial transforma-tion, or otherwise includes a time reversal. We introducea Z index s ( g ) which equals 0 in the former, and 1 inthe latter. Class [ I, s =0] is exemplified by BT ⊥ being amirror/glide-invariant plane, and [ I, s =1] by g = T i , whichis the composition of time reversal ( T ) with spatial in-version ( i ); all class-I symmetries are order two. Class-IIelementary orbits are likewise distinguished by whether g inverts time; they are additionally distinguished bywhether g acts on k ⊥ as a two-dimensional rotation( u =0), or as a two-dimensional reflection ( u =1). Equiv-alently, given that o is clockwise-oriented, u ( g ) distin-guishes between symmetries that preserve ( u =0) or invert( u =1) this orientation. In each of II-A and II-B, thereare then four classes of elementary orbits distinguishedby s, u ∈ Z . This gives ten classes of elementary orbits intotal, whose defining characteristics are summarized inthe first three columns of Tab. I. u s Symmetry constraints λ (I) ∀ k ⊥ , A =¯ g A ¯ g -1 ¯ g = e iπFµ - i k · R − k ⊥ = g ◦ k ⊥ A =¯ g A ∗ ¯ g -1 (¯ gK ) = e iπFµ - i k · R e i (cid:80) aλa ∈ R (II-A) 0 0 A =¯ g A ¯ g -1 ¯ g N = A ± N/L e iπFµ − k ⊥ ∈ o , A =¯ g A ∗ ¯ g -1 (¯ gK ) N = A ± N/L e iπFµ e i (cid:80) aλa ∈ R | o | = | g ◦ o | A =¯ g A -1 ¯ g -1 ¯ g N = e iπFµ - i k · R e i (cid:80) aλa ∈ R A =¯ g A t ¯ g -1 (¯ gK ) N = e iπFµ - i k · R − (II-B) 0 0 A i +1 =¯ g i A i ¯ g -1 i ¯ g N . . . ¯ g = e iπFµ - i k · R { λ i +1 a } = { λ ia } k ⊥ ∈ o , A i +1 =¯ g i A ∗ i ¯ g -1 i ¯ g N K . . . ¯ g K = e iπFµ - i k · R { λ i +1 a } = { - λ ia }| o |(cid:54) = | g ◦ o | A i +1 =¯ g i A -1 i ¯ g -1 i ¯ g N . . . ¯ g = e iπFµ - i k · R { λ i +1 a } = { - λ ia } A i +1 =¯ g i A t i ¯ g -1 i ¯ g N K . . . ¯ g K = e iπFµ - i k · R { λ i +1 a } = { λ ia } TABLE I. The first three columns distinguish between tenclasses of elementary orbits. The map of k ⊥ under g is g ◦ k ⊥ =( − s ( g ) ˇ g ⊥ k ⊥ , with ˇ g ⊥ a two-by-two orthogonal ma-trix that represents the point-group component of g in theplane orthogonal to the field; s ( g )=0 if g is purely a spatialtransformation, and =1 if g inverts time. | o | = | g ◦ o | indicatesthat o is mapped to itself under g , modulo a change in orienta-tion. u ( g ) distinguishes between proper and improper trans-formations on k ⊥ : ( − u :=det ˇ g ⊥ . Fourth and fifth columnsdescribe how unitary matrices ¯ g (that represent the symme-try g ) constrain the propagator. Column six summarizes theconstraints on λ a ; if there are none, we indicate this by − . 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 com-posed of L disconnected orbits which are mutually con-strained as g ◦ o i =( − u o i +1 and o i + L := o i . To clarify, o and { o j } Lj =1 are all clockwise-oriented, and − o i denotesan anticlockwise-oriented orbit. L was introduced in thesecond paragraph of Sec. IV, and is more precisely de-fined 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 six-foldrotation. B. Symmetry constraints on the propagator A Column four summarizes how g constrains A [I,II-A]and {A j } [II-B], which are respectively the propagatorsfor the self-constrained o and mutually-constrained { o j } .The corresponding spectra of the propagators are de-noted as { e iλ a } Da =1 and { e iλ ja } Da =1 .The unitary matrices ¯ g that constrain these propaga-tors form a projective representation of the pointgroup P ⊥ , as summarized in column five. Any g that in-verts time [ s =1] has the antiunitary representation ¯ gK ,with K implementing complex conjugation; otherwise[ s =0], g has the unitary representation ¯ g . The relationsin column five are closely analogous to the space-grouprelations satisfied by g : g N = e µ t R , with µ ( g ) ∈{ , } . (11)Here, N ( g ) ∈ N is the smallest integer such that g N is atranslation ( t ) by the lattice vector R , possibly composedwith a 2 π rotation (denoted e ). Note the similarity in def-inition of N with L , as was defined in Sec. IV A; in gen-eral, L divides N . R is nonzero for nonsymmorphic sym-metries such as screw rotations and glide reflections; the translation t R is represented on Bloch functions withwavevector k by the phase factor e − i k · R [cf. column five]. e in Eq. (11) is represented in column five by a phasefactor [( − F ] that a wavefunction accumulates upon a2 π rotation; F =0 (resp. =1) for integer-spin (resp. half-integer-spin) representations. The former case is usefulin analyzing solids with negligible spin-orbit coupling, aswe will exemplify with a case study of graphene in Sec.V A.The constraints in columns four and five are de-rived from the symmetry transformation of the Berryconnection and the one-form A · d k of Eq. (3). Thelatter may be expressed through Hamilton’s equation[ (cid:126) ˙ k = −| e | v × B / (cid:126) c ] as − M · B dt/ (cid:126) , with M ( k ) an orbitalmoment that transforms under g like the spatial compo-nents of a (3+1)-dimensional pseudovector M (cid:12)(cid:12) g ◦ k = ( − s ( g ) det[ˇ g ] ¯ g K s ( g ) (ˇ g M ) K s ( g ) ¯ g − (cid:12)(cid:12) k . (12)ˇ g here is a three-by-three orthogonal matrix representingthe point-group component of g , and K s M K s = M ∗ if s =1. C. Symmetry constraints on λ From taking the determinant of each equation in col-umn four, we derive constraints on λ a which are summa-rized in column six.Three of six classes in [I,II-A] are characterized by thereality condition e i (cid:80) aλa ∈ R . This implies λ =0 or π fora nondegenerate band, i.e., the orbit respectively encir-cles an even or odd number of Dirac points. λ = π isexemplified by the Dirac surface state of the Z topo-logical insulator Bi Se ; an orbit on the Dirac coneis self-constrained by time-reversal ( T ) symmetry [classII-A, u =0 , s =1]. λ = π may result solely from crystallinesymmetry, as exemplified by a mirror-symmetric orbit[class II-A, u =1 , s =0] that encircles a surface Dirac coneof the topological crystalline insulator SnTe [furtherelaborated in Sec. V B]. Finally, λ = π can be protectedby a composition of T and crystalline symmetry [classI, s =1], as exemplified by Weyl metals having the samespace group as WTe ; the robustness of λ depends sen-sitively on the field orientation with respect to the crys-tallographic symmetry axis, as elaborated in Sec. V C.For degenerate bands, the reality condition fixes (cid:80) Da =1 λ a to 0 or π , but not λ a individually. This willbe exemplified by our case studies of 3D Dirac metalsand Z topological insulators, in Sec. V F and Sec. V Grespectively. In Sec. V D, we further demonstrate thereality constraint may be further strengthened for spin-degenerate orbits ( D =2) that are self-constrained by T symmetry – into a zero-sum rule: λ + λ =0 mod 2 π .Moreover, if one considers the set of { λ } contributed byall closed orbits in a T -symmetric solid, we propose that { λ } comprise only pairs of λ , such that each pair individ-ually satisfies a zero-sum rule, as elaborated in Sec. V E.This global constraint on { λ } might be viewed as an ana-log of the Nielsen-Ninomiya fermion-doubling theorem, which states that there are always an equal number ofleft- and right-handed Weyl fermions on a lattice.The four II-B classes are characterized by { λ i +1 a } D − a =0 = { (-1) u + s λ ia } D − a =0 ; i ∈ Z L . (13)Let us discuss the implications of the above equation forthe three following cases (i-iii):(i) For s (cid:54) = u and even L (which necessarily holds if u =1),the minimal Landau-level degeneracy is L/
2. As illus-tration, the two valley-centered Fermi surfaces ( L =2) inthe transition-metal dichalchogenide WSe are mutuallyconstrained by T symmetry [class II-B, u =0 , s =1]; theLandau levels are nondegenerate but exhibit a symmetricsplitting, as elaborated in Sec. V A.(ii) For s (cid:54) = u and odd L , the minimal Landau-leveldegeneracy is L ; if D is also odd, it is necessary that { λ } contains either 0 or π . This is exemplified by threeof four disconnected Fermi pockets ( L =3) per valley ofbilayer graphene (with trigonal warping); we refer tothe three Fermi pockets that are mutually constrained bythe composition of T and six-fold rotational symmetry c z [class II-B, u =0 , s =1]. For pedagogy, it is instructiveto consider a model of bilayer graphene with spinlesselectrons, hence the orbit associated to each Fermipocket is nondegenerate ( D =1). The three-fold degen-erate λ = π reflects that each spinless pocket encircles aDirac point. (iii) We would demonstrate that the same set of Fermipockets, if equipped with multiple point-group symme-tries, can belong to multiple symmetry classes in Tab. I.For example, the above-mentioned three pockets are also invariant under the three-fold rotational symmetry c z ( u = s =0). We may therefore apply Eq. (13) with u = s ,which generally implies that the the minimal Landau-level degeneracy is L . Again, L =3 in the present exam-ple. We remind the reader that this three-fold degeneracyconstraint was consistently implied by T c z ; in addition, T c z implies a stronger constraint that { λ } contains ei-ther 0 or π . V. CASE STUDIES
The utility of Tab. I is illustrated in the following casestudies of existing conventional and topological metals,which were introduced in the previous section [Sec. IV C].
A. Orbits mutually constrained by time-reversalsymmetry: Application to graphene and transitionmetal dichalcogenides
Our first study encompasses materials with two time-reversal-related valleys in their band dispersion, as ex-emplified by graphene and monolayer WSe . We willdemonstrate: (i) how orbits in time-reversal-invariant( T ) solids can nevertheless develop a nonzero orbitalmagnetic moment, and (ii) the role of point-group sym-metry in discretizing the Berry phase of valley-centeredorbits.To explain (i), we point out that T symmetry relatesthe magnetic moment of wavepackets at k and − k [cf. Eq.(12)]; this mapping in k -space distinguishes the symme-try transformation of magnetic moments in solids fromthat in atoms. This allows for valley-centered orbits thatare separated in k -space to individually develop a mag-netic moment – since time reversal relates one valley tothe other, the net magnetic moment must vanish.In more detail, let us consider a finite chemical po-tential (as measured from the Dirac point) where valley-centered orbits ( o i ) are disconnected; we introduce herea valley index i ∈{ , } . The two orbits are mutually con-strained as T ◦ o = o [class II-B, u =0 , s =1]; each of o i isalso self-constrained as T i ◦ o i = o i [class I,s=1]. T i im-poses reality for the orbital component of wavefunctionsat each k , leading to M =0 [cf. Eq. (12) with D = s =1and ˇ g = − I ]. In analyzing graphene (which has negligi- ble spin-orbit coupling), it is useful to first neglect thespinor structure of its wavefunction and then account forthe Zeeman effect after – this was implicit in our pre-vious assumption of D =1. Such ‘spinless’ wavefunctionstransform in an integer-spin representation ( F =0 in Tab.I); the corresponding A in Tab. I should be interpretedas Eq. (2) without the Zeeman term.The Berry phase of graphene is π , and therefore˜ λ i := φ i R + φ i B =0+ π for each valley labelled by i =1 ,
2, asis consistent with the reality constraint in [class I, s =1]of Tab. I. We have added an accent to ˜ λ i to remindourselves that it is the purely-orbital contribution to λ i .Further accounting for the Zeeman contribution, λ i ± = ˜ λ i ± π g m c m , (14)with ˜ λ i = π , ± distinguishing two spin species, and m c the cyclotron mass. The symmetric splitting of λ ± about π implies an invariance under inversion: { λ i + , λ i − } = {− λ i + , − λ i − } mod 2 π . T symmetry imposesthe mutual constraint: { λ , λ − } = {− λ , − λ − } mod 2 π, (15)which follows from [class II-B, u =0 , s =1] of Tab. I. Eqs.(14-15) together imply { λ , λ − } = { λ , λ − } . To reca-pitulate, we have reproduced the well-known fact thatgraphene’s Landau levels are valley-degenerate but spin-split by the Zeeman effect.The valley degeneracy of the Landau levels may besplit by breaking spatial-inversion ( i ) symmetry, e.g.,with a hexagonal BN or SiC substrate. In zero field,the i asymmetry is predicted to produce a band gap of 53and 260 meV respectively – this may be interpreted asa nonzero Semenoff mass for the Dirac fermion. Sinceeach valley-centered orbit is no longer self-constrainedby T i , it develops a nonzero orbital moment (as was firstnoted in Ref. 61), as well as a non-integer φ i B /π . Conse-quently, Eq. (14) remains valid with ˜ λ i deviating from π ,as we confirm with a first-principles calculation in Fig.1(a). While { λ i + , λ i − } (that is associated to one valley)is no longer invariant under inversion, we remark thatthe invariance persists for the full set: { λ , λ − , λ , λ − } ,owing to the mutual constraint by T symmetry [cf. Eq.(15)].The predicted Semenoff masses induced by substratestend to be small – it is instructive to compare grapheneto WSe , a transition metal dichalcogenide with a largeSemenoff mass due to the natural absence of inversionsymmetry in its space group. WSe is similar to graphenein that its low-energy bands (at zero-field) are also cen-tered at two valleys, but differs from graphene in thatits bands are spin-split by spin-orbit coupling. Fig. 1(b)illustrates our calculated λ as a function of energy, with i =1 , λ = − λ owing to T sym-metry; at 0 . λ − λ ≈ π/ FIG. 1. (a) Characterization of a single valley of non-centrosymmetric graphene on hBN substrate, assuming atranslation-invariant, commensurate phase. Left inset: plotof φ B (blue), φ R (orange), ˜ λ := φ B + φ R (green) and λ ± [cf. Eq.(14)] (red) against energy. Right inset: plot of the orbit area( S ) vs E . (b) Characterization of a single valley of mono-layer WSe . Left inset: plot of φ B (blue), φ R (orange), and λ := φ B + φ R + φ Z (red) vs E . Right inset: plot of S ( E ). (a-b)are obtained by ab initio calculations detailed in App. D. B. Orbits self-constrained by mirror/glidesymmetry: Application to topological crystallineinsulators λ/π ∈ Z may originate from crystalline symmetry alone.This occurs for orbits that are self-constrained by a mir-ror/glide symmetry [IIA, u =1 , s =0], which may arise intime-reversal-invariant and magnetic solids. Any circu-lar orbit in [IIA,1 ,
0] intersects a mirror/glide-invariantline at two points denoted by k a and k b ; at each point,the assumed-nondegenerate band transforms in either theeven or odd representation of mirror/glide. The real-ity condition for e iλ in [IIA,1 ,
0] of Tab. I has a sim-ple interpretation: λ =0 if the representations at k a and k b are identical, and π if the two representations aredistinct. The former is exemplified by a band thatis nondegenerate at all k ⊥ bounded by o , which impliesthat o is contractible to a point – due to continuity ofthe representation along the high-symmetry line, the rep-resentations at k a and k b must be identical. λ = π oc-curs iff there is an odd number of linear band touchingsalong the segment of the mirror line within o – at eachtouching (a Dirac point), the mirror/glide representationflips discontinuously. λ = π is exemplified by the mirror-symmetric surface states of the SnTe-class of topolog-ical crystalline insulators, as well as by glide-symmetricDirac cones in band-inverted, nonsymmorphic metals. C. Effect of field orientation on the crystallinesymmetry of extremal orbits: Application to 3DWeyl fermion
The previous examples demonstrate that deformingthe crystal structure offers a way to tune λ . For 3Dsolids, we may continuously tune between integer andnon-integer λ/π without explicitly breaking any symme-try – by modifying the orientation of the field with re-spect to the crystal structure, we effectively modify thesymmetry of the extremal orbit. For specific orientations, the extremal orbit may be invariant under a point-groupsymmetry that stabilizes integer-valued λ/π . We shallillustrate this with a 3D topological metal whose type-IWeyl points lie on a high-symmetry plane that is invari-ant under T s z (the composition of time reversal witha two-fold screw rotation), as exemplified by strainedWTe . When the field is aligned parallel to the screwaxis, the maximal orbit of a Weyl point lies within the T s z -invariant plane [class-I, s =1], and is characterizedrobustly by λ = π . As illustrated in Fig. 2, λ deviatesfrom π when the field is tilted away, owing to the ab-sence of any symmetry for the tilted maximal orbit. FIG. 2. (a) Fermi surface of a Weyl fermion centered ona generic wavevector in a high-symmetry plane; red circleindicates an extremal orbit that has no symmetry. (b) For amodel of the Weyl fermion, we plot λ as a function of the fieldorientation, which we parametrize by spherical coordinatesillustrated in (a). D. Orbits self-constrained by time-reversalsymmetry
Let us consider orbits ( o ) which are self-constrainedas T ◦ o = o [II-A, u =0 , s =1]; these are orbits that lie in T -invariant cross-sections of the Brillouin torus, and en-circle Kramers-degenerate points. The contexts in whichself-constrained 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 λ ± , where ± distinguishesthe two spin species; recall that ˜ λ and λ ± differ onlyby the Zeeman splitting, as shown in Eq. (14). Theorbit-averaged orbital moment vanishes, and ˜ λ =0 reflectsthe trivial Berry phase of band extrema at T -invariantwavevectors. Combining ˜ λ =0 with Eq. (14), we obtainthe following zero-sum rule: λ + + λ − =0 mod 2 π .(b) This zero-sum rule is also satisfied for spin-degenerateorbits in spin-orbit-coupled solids: λ + λ =0 mod 2 π .(c) Both the orbital moment and Zeeman effect aver-age out, and λ = π reflects the nontrivial Berry phaseassociated to the Kramers-degenerate Dirac point at T -invariant wavevectors. Demonstration
In all cases, the propagator A [ o ] satisfies¯ T A ∗ ¯ T -1 = A , ¯ T ¯ T ∗ = ( − F A -1 , ¯ T -1 = ¯ T † , (16)with ¯ T K an antiunitary representation of T [cf. row 4of Tab. I with N = L =2]. The second equation may becontrasted with the usual antiunitary representation (de-noted ˜ T K ) of T , which satisfies ˜ T ˜ T ∗ =( − F ; the addi-tional factor of A -1 in Eq. (16) indicates that Eq. (16)represents an extension of the magnetic point group bycontractible loop translations (represented by A ) in k -space, which generalizes a previous work on noncon-tractible loop translations. The first equation in Eq. (16) impliesdet A = e i (cid:80) Da =1 λ a = ±
1, with D =2 and 1 for cases(b) and (c) respectively; in case (a), we should interpretdet A = e i ˜ λ = ±
1. For cases (a) and (b), we may furtherrestrict det A to +1 by the following argument: whilepreserving det A , contract o T -symmetrically to aninfinitesimal loop that encircles the T -invariant point.In the classes of (a-b), the band dispersion is extremizedat T -invariant points, hence the band velocity v ( k )vanishes along the infinitesimal loop. det A is thussolely contributed by the non-geometric one forms whichdepend inversely on the velocity. Further applyingthe identity: det exp (cid:82) V ( k ) | d k | =exp (cid:82) Tr V ( k ) | d k | forany matrix V ( k ), and the time-reversal constraints:Tr M ( k )= − Tr M ( − k ) [cf. Eq. (12) with g = T ] andTr σ ( k )= − Tr σ ( − k ), we derive the desired result. Theabove proof required path-ordering and matrix traces incase (b), where A , M and σ are two-by-two matrices;such matrix operations are not needed in case (a).In case (c), we contract o T -symmetrically to an in-finitesimal loop encircling the Kramers-degenerate Diracpoint – since the velocity is finite in this limit, the non-geometric contribution to λ vanishes; what remains is the π -Berry-phase contribution. This completes the demon-stration.One implication of (a-b) for spin-degenerate orbits thatare self-constrained by T symmetry: the net phase offsetΘ [cf. Eq. (10)], of two time-reversal-related fundamentalharmonics, can only be 0 or π . The former occurs if thetwo values of { λ } are closer to 0 than to π , and vice versa. E. Global constraint on { λ } fortime-reversal-symmetric solids Let us impose a global constraint on the set [denoted { λ } ] of all λ ia that are contributed by closed orbits ofbulk states in a T -symmetric solid. To clarify, for a d -dimensional solid, a ‘bulk state’ is spatially extended in d directions. By combining our results for orbits that areself-constrained [cf. Sec. V D] and mutually constrained[cf. Sec. V A] by time-reversal ( T ) symmetry, we wouldshow that { λ } comprise only of inversion-invariant pairs,i.e., pairs that are symmetric about zero. In our countingof ‘pairs’, each closed orbit with energy degeneracy D contributes D values of λ , independent of whether anyof these D values are mutually degenerate, or degeneratewith λ from a distinct orbit. A corollary of this result isthe global sum rule: (cid:80) λ =0 mod 2 π .Let us first prove this claim for spin-degenerate orbits.All orbits in a T -symmetric solid are either self- or mu-tually constrained by T symmetry. As proven in Sec.V D, each self-constrained orbit contributes an inversion-invariant pair { λ, − λ } . Utilizing [class II-B, u =0 , s =1]of Tab. I, the net contribution of any T -related pairof spin-degenerate orbits is an inversion-invariant quar-tet { λ , λ , − λ , − λ } ; this was exemplified by our casestudy of graphene in Sec. V A.We would next prove the global constraint for spin-orbit-coupled solids lacking spatial-inversion symmetry[case (c) in Sec. V D]. By similar reasoning as inthe previous paragraph, each T -related pair of spin- non degenerate orbits contributes { λ, − λ } . A self-constrained orbit necessarily contributes λ = π , as provenin Sec. V D. Furthermore, self-constrained orbits alwayscome in pairs; this follows because each self-constrainedorbit encircles a T -invariant wavevector, and there arealways 2 d such wavevectors in a d -dimensional lattice.This completes the proof.It is instructive to draw analogies between our resultand the Nielsen-Ninomiya fermion-doubling theorem (asapplied to T -invariant solids). In full generality, thetheorem states that there are always an equal numberof left- and right-handed Weyl fermions on a lattice; asimple implication is that Weyl fermions come in pairs.If Weyl fermions exist at wavevectors which are not T -invariant, they necessarily come in pairs due to T sym-metry mapping k to − k . In the absence of point-groupsymmetry in the space group of the lattice, bands alwayscome in pairs that touch at isolated T -invariant wavevec-tors – such Kramers-degenerate points are also Weylpoints. That there are even number of Weyl fermions ona lattice thus complements our previous observation thatthere are an even number of Kramers-degenerate points.If the Fermi surface encircles Kramers-degenerate Weylpoints instead of intersecting them, we recover our pre-vious claim that self-constrained orbits come in pairs.To recapitulate, the global zero-sum rule for { λ } is ex-emplified by all case studies of bulk orbits in this work;more generally it constrains the bulk oscillatory phenom-ena of all T -symmetric solids.Let us extend our discussion to closed orbits con-tributed by surface states of a 3D T -symmetric solid;these orbits lie in a 2D surface Brillouin zone. By a ‘sur-face state’, we mean a state localized spatially to a surfacethat is translation-invariant in two directions; this surfacelies at the interface between vacuum and a bulk that issemi-infinite in the direction orthogonal to the surface. Ifthe solid is spin-orbit-coupled, surface bands are generi-cally nondegenerate, and we would use nearly the sameargument (given above) for spin-nondegenerate, bulk or-bits; the sole difference is that the fermion-doubling the-orem does not apply, and it is possible to have an oddnumber of self-constrained orbits (associated to an oddnumber of surface Dirac fermions). One implication isthat χ s := (cid:80) λ (summed over surface orbits) may equal0 or π . If the bulk of the solid is insulating, χ s maybe viewed as a Z index that classifies insulators fromthe perspective of its surface magnetotransport; this isequivalent to the Z strong classification of 3D insu-lators in Wigner-Dyson class AII. If the solid has negli-gible spin-orbit coupling, then the surface orbits satisfythe same type of global constraint as for bulk orbits –to demonstrate this, one may utilize the above argumentfor spin-degenerate bulk orbits. F. Spin-degenerate orbits in spin-orbit-coupled,centrosymmetric metals: Application to the 3DDirac metal Na Bi A zero-sum rule also applies individually to each T i -invariant orbit ( o in class [I, s =1]); this rule applieswhether or not the orbit is self-constrained by T sym-metry. We have previously discussed such orbits in thecontext of graphene; here, we extend our discussion to T i -invariant, spin-degenerate orbits in SO-coupled solids.The reality of det A and contractibility of o together im-ply λ + λ =0 mod 2 π .Let us apply this result to a 3D massless Diracfermion – a four-band touching point between conically-dispersing bands which are spin-degenerate at genericwavevectors (owing to T i symmetry). We shall assumethat the 3D Dirac point is not centered at a T -invariantpoint, hence an orbit at finite chemical potential (as mea-sured from the Dirac point) is T i -invariant but not T -invariant, e.g., such Dirac points in the topological metalNa Bi are stabilized by three-fold rotational symmetry. 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 . λ , ≈± π/ k z , as indicated bya dashed line on the plot). In Fig. 3(b), we further plot λ a ( µ, ¯ k z ) for a range of µ below the Dirac point. Theab-initio calculation of λ a is detailed in App. D.Let us discuss the experimental implications for Na Bi.A single Dirac point contributes to the oscillatory mag-netization a term of the form Eq. (4), with D =2 and λ , ( µ, ¯ k z ) plotted in Fig. 3(b). Na Bi actually has two T -related Dirac points, which both lie on the rotation-invariant line through Γ. Assuming the Fermi surfacesof the two Dirac points are disconnected, the effect of thesecond Dirac point is simply to double the magnitude ofthe oscillations. This follows because the maximal orbits
FIG. 3. Characterization of the 3D Dirac metal Na Bi. (a)Left: surface of fixed energy ( ¯ E = − .
08 eV measured fromthe 3D Dirac point); contours of fixed energy and k z are il-lustrated by red circles. Right: plot of { λ a ( ¯ E, k z ) } a =1 , fora field parallel to the trigonal axis ( k z ). k z is measured inunits where the reciprocal period (0 . -1 ) is unity; the min-imum (resp. maximum) k z on the plot corresponds to thesouth (resp. north) poles of the constant-energy surface. Theextremal orbit occurs at k z = ¯ k z ( ¯ E ) indicated by a dashedline. (b) λ , ( E, ¯ k z ( E )) vs E , for a range of E below the Diracpoint. ( o , o ) on disconnected Fermi surfaces are mutually con-strained by T symmetry, hence from class [II-B, u =0 , s =1]of Tab. I, { λ , λ } = {− λ , − λ } mod 2 π . We remind thereader that λ i = − λ i owing to T i symmetry, hence in com-bination { λ , λ } = { λ , λ } .To reiterate, the two Dirac points contribute identi-cally to the oscillatory magnetization, which we plot inFig. 4 for a temperature of 0.37 K and three representa-tive Fermi energies. These plots are obtained from Eq.(4), which receive as inputs: (i) the ab-initio calculatedband parameters { S, S zz , m c , λ , } for the maximal orbit,as well as the Dingle lifetime τ =10 − s. For the rangeof B plotted, kT / (cid:126) ω c ∼ − and ω c τ ∼
1, which impliesthat the Dingle factor (rather than temperature smear-ing) is the limiting factor in observing higher harmonics;for the parameters considered, the third harmonic doesnot appreciably modify the plots. The second harmonic isespecially transparent at µ c := − λ = π/ − λ [cf. Fig. 3(b)] leads to the complete de-structive interference of all time-reversed pairs of oddharmonics (including the fundamental harmonic) [cf. Eq.(9)], i.e., the dHvA period is effectively halved. The pointof destructive interference intermediates two energy in-tervals where the phase offset Θ [cf. Eq. (10)] in thefundamental harmonic differs by π : Θ= π for µ>µ c andcloser to the Dirac point [cf. green shaded region in Fig.3(b)]; Θ=0 for µ ≤ µ c [blue shaded region]. Our case studydemonstrates that the experimental tunability of µ in Na Bi should expose a wide range of T -symmetricinterference patterns.For the parameter ranges in Fig. 3 and Fig. 4, theFermi surfaces enclosing each Dirac point are indeed dis-connected. As µ is increased beyond ∼
10 meV, the twoFermi surfaces intersect at a saddlepoint and merge intoa dumbbell Fermi surface. For a field parallel to the trig-onal axis, this merged Fermi surface has a single minimalorbit o (encircling the saddlepoint), in addition to theaforementioned two maximal orbits. At slightly higher0 FIG. 4. Magnetization oscillations calculated for Na Bi us-ing Eq. (4) at T=0.37K, τ =10 − s and Fermi energies (a) -72meV, (b) -69meV, and (c) -66meV. The corresponding λ , are indicated by dashed lines in Fig. 3(b). The sum of thetwo lowest harmonics [ r =1 , r =1, dashedblue line) is inadequate. The unit of magnetization is chosenas 10 − µ B / ˚A , with µ B the Bohr magneton. energy ( ∼
20 meV), the dumbbell transforms into an el-lipsoid, and correspondingly the three extremal orbitsmerge into a single maximal orbit (denoted o ). o (andalso o ) is T -invariant [under class II-A, u =0 , s =1] and T i -invariant [class I, s =1]; this exemplifies how a singleorbit having multiple symmetries may fall under differ-ent classes in Tab. I. Both classes impose (consistently)a zero-sum-rule constraint on { λ } .
1. Comment on magnetotransport experiments of 3D Diracmetals
We conclude this section by commenting on the inter-pretation of several magnetotransport experiments on 3DDirac metals.
The experimental data in these workshave only been fitted to the fundamental harmonic ( r =1)in the traditional Lifshitz-Kosevich formula. The phaseoffset in this fitted fundamental harmonic is interpretedas φ B − φ M ± π/ φ B = π that is inferred from such a fittingis viewed as a smoking gun for a 3D Dirac metal. Theseinterpretations have been justified by appealing tothe one-band Onsager-Lifshitz quantization rule [Eq. (1)with D =1].We propose an alternative interpretation that never-theless makes sense of the robustness of their measure-ments, across different experimental groups and slightly-varying samples – this robustness is solely a consequenceof T symmetry, rather than anything intrinsic to a 3D-Dirac Fermi surface. In other words, Θ= π should not beviewed as a smoking gun for a 3D Dirac metal.From our perspective, spin-degenerate orbits in spin-orbit-coupled metals are semiclassically described by themulti-band quantization rule [cf. Eq. (1)] and a gener-alized Lifshitz-Kosevich formula [cf. Eqs. (4-8)]. Thephase corrections in these expressions encode the non-abelian Berry, orbital-moment and Zeeman one-forms [cf.Eq. (2)]. Owing to the generic non-commutivity of thethree one-forms in the path-ordered integral [cf. Eq. (2)], λ ± cannot be decomposed into a sum of Berry, orbital- moment and Zeeman phases. As described in Eqs. (4-10), the phase offset (Θ) of thesummed fundamental harmonic, after substracting theMaslov and Lifshitz-Kosevich corrections, is a suitably-defined average of λ , . As proven in Sec. V D, Θ is re-stricted to 0 or π for orbits that are self-constrained by T symmetry. We believe that previous measurements of Θ have been misinterpreted as a measurement of theBerry phase. The convergence in measurements of sev-eral groups indicate the quantization of Θ due to T sym-metry. However, Θ= π is not an intrinsic property of 3DDirac metals, as we have exemplified by our case studyof Na Bi. Moreover, we point out that Θ= π has beenmeasured in conventional (i.e., non-topological) metals,even in the first metal (3D, single-crystal Bismuth) forwhich quantum oscillations were reported. As far backas 1954, it was remarked by Dhillon and Shoenberg that:‘For Bismuth good agreement was found by the theoreti-cal formula except that the signs of the fundamental andodd harmonics had to be reversed.’ Here, the ‘theoret-ical formula’ refers to the traditional Lifshitz-Kosevichformula without the Θ correction to the phase. In alater work, the quantization of Θ is implicit in the formof Shoenberg’s spin ‘reduction factor’, which involvesa phenomenological g-factor in the presence of spin-orbitcoupling.One contribution of this work is to present an analytic formula for the effective g-factor through the spectrumof Eq. (2). Our proposed formula is calculable from ab initio wavefunctions [as elaborated in App. D], whichallows for a quantitative comparison between theory andexperiments. Another contribution of this work is therecognition that the quantization of Θ originates from T symmetry; hence in magnetic materials, Shoenberg’sspin ‘reduction factor’ is not generally applicable, but ourgeneralized Lifshitz-Kosevich formulae in Sec. III remainvalid.We further emphasize that Θ is an incomplete char-acterization of the Fermi-surface orbit; a more completecharacterization involves measuring the individual valuesof λ , , which carry the complete non-abelian informationabout the Fermi-surface wavefunction. Two methods areavailable to measure λ , : (i) where either of ω c τ (cid:28) kT (cid:29) (cid:126) ω c applies, one needs to measure both the ampli-tude and the phase offset of the fundamental harmonic; λ , is then extracted from the r =1 terms in Eq. (4) andEqs. (9-10). (ii) Where neither of the above inequali-ties apply, the experimental data should be fitted to thefull harmonic expansion of Eq. (4), and λ , extractedaccordingly. Note that the zero-sum rule for λ , impliesthat there is only one parameter to be fitted.Going beyond general symmetry constraints, one mayask if Θ and λ , may be analytically calculated in a case-by-case basis. Given material-specific information aboutthe magnetic point group of the 3D Dirac point, as wellas the relevant symmetry representations that span thelow-energy Hilbert space at the Dirac point, one mayindeed determine Θ from a k · p analysis, as exemplified1in the next section. G. Calculating λ , from k · p analysis: Applicationto 3D massive Dirac fermions in Z topologicalinsulator Bi Se In certain symmetry classes, and for sufficiently small(but nonzero) chemical potentials ( µ ) relative to theDirac point, 3D massless Dirac fermions may exhibit λ ≈ λ ≈ π ; this would imply from Eq. (10) that Θ= π .Precisely, we restrict | µ | to be small enough that a lin-earized k · p Hamiltonian is a good approximation – howsmall | µ | needs to be to fulfil the above criterion dependson material-specific band parameters. As far as the lin-ear approximation is valid, λ , that follows from thesubsequent k · p analysis is approximately π , hence ouruse of ≈ throughout this section. λ , should manifestin the dHvA oscillations for sufficiently weak field, suchthat the semiclassical approximation remains valid, i.e., l S ( µ ) (cid:29) Let us demonstrate how λ , ≈ π arises for 3D mass-less Dirac fermions centered at a wavevector (Γ) witha magnetic point group ( P ) that combines T symmetrywith the D d point group. λ , ≈ π relies not just on P , but also on symmetries that are absent in P ; theseadditional symmetries emerge only at long wavelength( k → ), where the averaged-out crystal field has higher symmetry than the magnetic point group.Our study of D d is motivated by the identicalsymmetry class of Bi Se , a well-known Z topologi-cal insulator. For the symmetry representations ofBi Se , the k · p Hamiltonian around Γ (of the 3D Bril-louin torus) assumes the form: H = (cid:126) (cid:0) v ( k x σ + k y σ )+ wk z σ (cid:1) ⊗ τ + M ( k ) I × ⊗ τ . (17)to quadratic accuracy in k . Here, σ j and τ j are Pauli matrices spanning spin and orbital (la-belled P + z and P − z ) subspaces respectively, and M ( k )=∆+ b k z + b ( k x + k y ) with ∆, the Dirac mass, van-ishing at the critical point between trivial and topologicalinsulator.Assuming ∆=0 for now, [ H ( k ) , I × ⊗ τ ]= O ( k ) is anemergent conservation law of the linearized H . Block-diagonalizing H with respect to I × ⊗ τ = ±
1, we derivetwo decoupled 3D Weyl Hamiltonians H ± with oppo-site chirality, each satisfying the time-reversal constraint T H ± ( k ) T -1 = H ± ( − k ) with T = iσ K squaring to − I ina half-integer-spin representation ( F =1). Independentof the field orientation, the spin-degenerate extremal or-bit effectively decouples to two nondegenerate ( D =1) ex-tremal orbits which are individually T -invariant – thisimplies from a previous demonstration [case (c) of Sec.V D] that λ , ≈ π , which may be viewed as the Berryphase ( φ B ) for each decoupled orbit. This would effec-tively imply a single set of dHvA harmonics (indexedby r in Eq. (4)) with twice the usual amplitude, and a Lifshitz-Kosevich phase correction of − π/ λ , ≈ π ± π ∆ /mv with m the free-electron mass, for a field parallel to the trigo-nal axis ( k z ). It is worth remarking that φ B reduces from π , and the phase ( φ R ) associated to the orbital momentincreases from 0, however their sum ( φ B + φ R ) remainsfixed to π . In fact, naturally-occurring Bi Se does not lie at thecritical point, i.e., the corresponding 3D Dirac fermionis massive: we estimate ∆ /mv ≈ .
13 utilizing fitted pa-rameters from an ab-initio calculation in Ref. 83; utiliz-ing Eq. (10) with | ∆ /mv | < .
5, we further derive Θ= π .This splitting of λ manifests as two sets of harmonics inthe dHvA oscillation, which should be measurable uti-lizing techniques outlined in Sec. V F. The tunabilityof ∆ by doping or strain suggests that the dHvAoscillations are correspondingly tunable – in particular, | ∆ /mv | =0 . π .There exist galvano-magnetic evidence that supportsthe quantization of Θ. The fundamental Shubnikov-de Haas harmonic in the transverse resistivity ofBi Se has been fitted as∆ ρ xx ∝ cos( l S +Θ+ φ M − π/ ≈ cos[ l S +2 π (0 . , (18)with the higher harmonics suppressed by temperature.The above proportionality may be derived from: (a) ρ xx ≈ σ xx /σ xy , which is valid for a large Hall angle σ xy (cid:29) σ xx , and (b) the proportionality between the oscil-latory components of the transverse conductivity and thelongitudinal magnetic susceptibility: ∆ σ xx ∝ m c ∆ χ ; the cyclotron mass ( m c ) is positive for the lone electron-like orbit in Bi Se , and the oscillatory susceptibility(∆ χ ) is obtained from differentiating Eq. (4) (with D =2)with respect to B and keeping only the fastest oscilla-tory terms. The proportionality stated in (b) is valid formetals having only a single extremal orbit, and for suf-ficiently weak fields such that µ (cid:29) (cid:126) ω c . The fit in Eq.(18) implies Θ=0 . − / / ≈ ≈ π/ π/ π/ Se , produced the same value for Θ.Our explanation for this robustness: independent of thefield orientation, the extremal orbit is invariant undertime-reversal symmetry, which quantizes Θ to 0 or π . Incontrast, crystalline symmetries of an extremal orbit de-pend sensitively on the field orientation, as explained inSec. V C.2 VI. DISCUSSION
In fermiology, the shape of the Fermi surface is de-ducible from the period of dHvA oscillations; here, wepropose that the topology of the Fermi-surface wavefunc-tions are deducible from the phase offset ( λ a ) of dHvA os-cillations in 3D solids, as well as in fixed-bias oscillationsof the differential conductance in scanning-tunneling mi-croscopy. λ a manifests as the complete , subleading [ O (1)]corrections to the Bohr-Sommerfeld quantization rulesfor closed orbits without breakdown [cf. Eq. (1)]; we haveformulated λ a as eigen-phases of a propagator A definedin Eq. (2).In certain solids and for certain field orientations withrespect to a crystal axis, (cid:80) Da =1 λ a (with D the degeneracyof the low-energy band subspace) are topologically invari-ant under deformations of the zero-field Hamiltonian thatrespect the symmetry of the orbit, as well as preserves theglobal shape of the orbit. Precisely, globally-equivalentorbits correspond to a graph with a homotopy equiva-lence defined in Ref. 29. To identify orbits with robustlyinteger-valued (cid:80) Da =1 λ a /π , as well as identify the degen-eracy of Landau levels, we classify symmetric orbits intoten (and only ten) classes summarized in the first threecolumns of Tab. I; the rest of the table describes thecorresponding 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 unitarygenerated by the Berry connection, and may be viewedas the purely-geometric component of A ; φ B,a are definedas the eigen-phases of W , and may be viewed as the non-abelian generalization of the Berry phase. W providesa purely-geometric characterization of bands that is inti-mately related to the topology of wavefunctions over theBrillouin torus. While W (resp. φ B,a ) genericallydiffers from A (resp. λ a ) due to the orbital moment andthe Zeeman effect, W and A transform identically undersymmetry, and therefore satisfy the same constraints. Inparticular, if we define { φ B } as the set of φ B contributedby all bulk orbits in a T -symmetric solid, then { φ B } com-prises only pairs of φ B which individually satisfy a zero-sum rule, in close analogy with the global constraint on { λ } that is described in Sec. V E.Higher-order (in | B | ) corrections to the Bohr-Sommerfeld quantization rule become relevant in higher-field experiments that intermediate the semiclassicaland quantum limits; these corrections may be in-terpreted as zero-field magnetic response functions. These corrections are accounted for in the general-ized Lifshitz-Kosevich formulae [cf. Eqs. (4-8)] by thesimple replacement l S ( E )+ λ a ( E ) → l S ( E )+ λ a ( E, | B | ) , with λ a expandable asymptotically in powers of | B | ; λ a ( E, | B | =0)= λ a ( E ) is obtained from diagonalizing thepropagator A [cf. Eq. (2)]. It would be interesting to ex-tend our symmetry analysis to the field-dependent com-ponent of the phase offset. ACKNOWLEDGMENTS
The authors thank Zhijun Wang and Ilya K. Drozdovfor their expert opinions on 3D Dirac metals and tunnel-ing spectroscopy. We also thank Yang Gao and Qian Niufor communicating their work on higher-order correctionsin the quantization rule. We acknowledge support bythe Yale Postdoctoral Prize Fellowship (AA), NSF DMRGrant No. 1603243 (LG), the Ministry of Science andTechnology of China, Grant No. 2016YFA0301001 (CWand WD), and the National Natural Science Foundationof China, Grants No. 11674188 and 11334006 (CW andWD).The appendix is organized as follows.(A) For bands of any degeneracy and symmetry, wederive the oscillatory magnetization and density of statesin the de Haas-van Alphen and related effects.(B) We derive the eigen-phases ( λ a ) of the propaga-tor for the 3D massive Dirac fermion with D d symmetry.(C) We prove the reality of the orbital component of e iλ for a class of two-band Hamiltonians, which generalizesa previous result for a certain class of massive Diracfermions in Ref. 98.(D) We detail the ab initio methods used in the calcu-lation of Fig. 1, and provide further details about thegraphene-on-boron-nitride case study.(E) We provide the k · p model Hamiltonian (of a Weylfermion) that was used in the calculation of Fig. 2. Appendix A: de Haas-van Alphen-type oscillationsfor bands of any degeneracy and symmetry
The following derivation of Eqs. (4-8) is a simple gen-eralization of Roth’s calculation for spin-degeneratebands in 3D solids. Consider a free-fermion system im-mersed in a magnetic field; we assume that near theFermi energy, bands are D -fold degenerate, and the Lan-dau levels { E a,j } a ∈{ ,...,D } are obtained semiclassicallyfrom Eq. (1). The grand-canonical potential is con-tributed by D sets of sub-Landau levels:Ω = − kT D D (cid:80) Da =1 (cid:82) ∞−∞ dk z (cid:80) ∞ j =0 φ ( E a,j ) , − kT D D (cid:80) Da =1 (cid:80) ∞ j =0 φ ( E a,j ) , (A1)for 3D and 2D metals respectively. The above expressioninvolves φ ( ε ) := log (cid:16) e − β ( ε − µ ) (cid:17) , (A2)3with β =1 /kT , and the Landau-level degeneracy factor(per unit volume, or area): D D := 14 π l , D D := 12 πl . (A3)Utilizing the standard Poisson summation and followingessentially Roth’s calculation in Ref. 16, the oscillatorycomponent of Ω is derived as δ Ω = (2 π ) / kT D D l | S zz | / × D (cid:88) a =1 ∞ (cid:88) r =1 cos (cid:2) r ( l S + λ a − φ M ) ± π/ (cid:3) r / sinh (2 π kT / (cid:126) ω c ) (cid:12)(cid:12)(cid:12)(cid:12) µ, ¯ k z , (A4)for 3D metals, and δ Ω = 2 kT D D D (cid:88) a =1 ∞ (cid:88) r =1 cos (cid:2) r ( l S + λ a − φ M ) (cid:3) r sinh (2 π kT / (cid:126) ω c ) (cid:12)(cid:12)(cid:12)(cid:12) µ , (A5)for 2D metals. Employing the thermodynamic defini-tion of magnetization as M = − ∂ Ω /∂B , and keeping onlythe fastest oscillatory term in the semiclassical limit( l S (cid:29) Employing the following relation: g ( µ )= − ∂ Ω /∂µ | T =0 between the density of statesand the zero-temperature grand-canonical potential, wedifferentiate Eq. (A4) [at zero temperature, and inclu-sive of the Dingle factor] twice to obtain the oscillatorycomponent of the 3D density of states: δg ( E ) = 1 √ π / (cid:126) ω c l | S zz | / × D (cid:88) a =1 ∞ (cid:88) r =1 e − rπω c τ cos (cid:2) r ( l S + λ a − φ M ) ± π/ (cid:3) r / (cid:12)(cid:12)(cid:12)(cid:12) E, ¯ k z . (A6)By convolving this quantity with the derivative of theFermi-Dirac distribution (as in Eq. (7)), we derive thatthe oscillatory component of G is Eq. (8). This derivationis aided by the identity − (cid:90) ∞−∞ dε f (cid:48) T ( ε − µ − E ) Φ( ε ) = πD sin πD Φ( µ + E ) , (A7)with D := kT ( ∂/∂µ ); the above differential op-erator is defined through the Taylor expansion: x/ sin x =1+ x / x / . . . , which is well-known inits application to the Sommerfeld expansion. Applyingthe above identity with Φ= δg , and retaining only thequickest oscillatory terms [i.e., only the terms derivedfrom D acting on the cosine function: πD sin πD cos(Γ) = πkT d Γ /dµ sinh( πkT d Γ /dµ ) cos(Γ) (cid:21) , (A8)we arrive at the desired result. Appendix B: 3D massive Dirac fermion with D d symmetry The linearized k · p Hamiltonian from Eq. (17) H ( k ) = (cid:126) (cid:0) v ( k x σ + k y σ )+ wk z σ (cid:1) ⊗ τ + ∆ I × ⊗ τ (B1)with a nonzero mass ∆. This Hamiltonian has an emer-gent reflection symmetry:˘ r z H ( k ⊥ , k z ) ˘ r -1 z = H ( k ⊥ , − k z ) , ˘ r z = − iσ τ , (B2)which is broken by terms which are cubic in k . If thefield is oriented parallel to the trigonal axis ( k z ), themaximal orbit lies at k z = 0 owing to the just-mentionedsymmetry. It is here that [ H ( k ⊥ , , ˘ r z ] = 0, and wemay block-diagonalize H with respect to even and oddrepresentations of ˘ r z : H ± = (cid:126) v ( k x γ + k y γ ) ± ∆ γ . (B3) H ± describe two 2D massive Dirac fermions with oppo-site chirality. It is known for each of H ± that φ R + φ B = π is independent of the Semenoff mass (∆); we rederivethis result in our language in in Sec. C, where we furtherclarify the applicability of their result to more generaltwo-band Hamiltonians. What remains is to calculatethe Zeeman contribution to λ . Employing that the spinoperator S z is represented (cid:126) γ z /
2, and v ⊥ = v (cid:126) | k | (cid:112) ( (cid:126) v k ) + ∆ , (B4)the Zeeman phase is φ Z = g (cid:126) m (cid:90) o | d k | v ⊥ ∆ (cid:112) ∆ + ( (cid:126) v k ) ≈ π ∆ mv . (B5)Employing the following parameters from Ref. 83: ∆ =0 . v = 6 . × ms − , we obtain φ Z ≈ . π . Appendix C: Reality of the orbital component of e iλ For a general dispersion (cid:15) ( k ), the area of constant en-ergy contour is a function of (cid:15) : S ( (cid:15) ). For a general scalarvalued function f ( k ), we define a F ( (cid:15) ) := (cid:82) S ( (cid:15) ) f ( k ) d k .Then the average of f on the contour is defined by: f := dFd(cid:15) / dSd(cid:15) = (cid:90) (cid:34) fv x dk y / dSd(cid:15) . (C1)Thus φ R = (cid:90) (cid:33) − l BMv x dk y = l BM dSd(cid:15) , (C2)where M is orbital magnetization. The sign changecomes from a change in orientation in path. For twoband Hamiltonian in the form of H ( k ) = (cid:32) ∆ r ( k ) e iθ ( k ) r ( k ) e − iθ ( k ) − ∆ (cid:33) , (C3)4energy dispersion are symmetric with respect to zero.Here, r is real and ∆ is k independent. In this case,orbital magnetization are simply related to Berry curva-ture by M = e (cid:126) c Ω (cid:15) , where Berry curvature is defined by Ω = ∇ × X ( k ). Then φ R = l eB(cid:15) (cid:126) c Ω dSd(cid:15) = Ω (cid:15) dSd(cid:15) = (cid:15) dφ B d(cid:15) , (C4)and φ R + φ B = (cid:15) dφ B d(cid:15) + φ B = d(cid:15)φ B d(cid:15) . (C5) d(cid:15)φ B d(cid:15) is evaluated to be ± π (cid:82) (cid:33) d k · ∇ k θ for conductionband and valence band respectively, which is fixed to 0or π (mod 2 π ).Any two-band Hamiltonian that can be transformedinto Eq. (C3) (i.e., with at least one k -independentmultiplicative coefficient of a Pauli matrix) is alsocharacterized by real e i ( φ R + φ B ) . We propose a sufficientcondition for the existence of this transformation, whichwe denote by U in what follows. Sufficient condition
In general, for any three-by-threespecial orthogonal matrix ( R ), there exists a k -independent unitary U such that U † d · σ U = ( R ( U ) d ) · σ , R ( U ) ∈ SO(3) , (C6)where d ( k ) is a real three-vector, and σ = ( σ , σ , σ )are Pauli matrices. Let us separate each component of d into its k -independent and -dependent parts: d j ( k ) = d j + d j ( k ). If d x , d y , d z are linearly-dependent functionsof k , such a transformation is always possible. Indeed,suppose d y and d z are linearly independent, and by as-sumption we might express d x = αd y + βd z for α, β ∈ R .Then the desired basis transformation maps d → d (cid:48) = α + β ( − d + αd + βd ) (C7)such that d (cid:48) is k -independent. The above coefficientsof d j may be viewed as a real three-vector with unitnorm. Therefore, the transformation from d → d (cid:48) might be viewed as a special-orthogonal rotation ( R )of a Euclidean coordinate system. Since R existsthat is k -independent, so thus U . This completes thedemonstration.As an example, if d ( k ) is at most linear in k , with k = ( k x , k y ) restricted to a Brillouin two-torus, theabove condition is satisfied owing to the existence ofonly two linearly-independent terms: k x and k y , andtherefore φ R + φ B is fixed to either 0 or π . This class oftwo-band, linearized Hamiltonians includes, as a specialcase, the massive-Dirac Hamiltonian [cf. Eq. (C3)] first discussed in Ref. 98. Appendix D: ab initio calculations leading to Fig. 1and further details about the graphene-hBN casestudy
Our ab initio calculations are carried out withinthe framework of density functional theory, as imple-mented in quantum espresso software package . Norm-conserving potentials and Perdew-Burke-Ernzerhof(PBE) exchange-correlation functional are used to de-scribe electron-ion, electron-electron interactions respec-tively. Most of the physical quantities are obtainedby Wannier interpolation, where maximally-localizedWannier functions are obtained by wannier90 code .For the graphene-boron-nitride heterostructure, vdWcorrection is further included. Calculation of the prop-agator A has been coded into python package owl, which has been developed for general-purpose Wannierinterpolation.When graphene is placed on a boron-nitride substrate,the Dirac point splits with an energy gap of approx-imately 60 meV. The cyclotron mass m c ( E ) [definedas ( (cid:126) / π ) dS/dE , with S ( E ) the k-space area of theconstant-energy band contour], when evaluated at theband edge, is approximately 1000 eV /c with c the speedof light; this mass is approximately 0.002 times the free-electron mass. Appendix E: k · p model of Weyl fermions with C T symmetry In spinor basis, T is represented by iσ y K , C z (two-foldrotation about (cid:126)z ) by iσ z ; we have picked a basis wherethe orbital component of the basis functions transformstrivially under C z . We consider a k · p Hamiltonianexpanded around a generic point on the k z = 0 plane; ithas the symmetry C T H ( k ) C − T = H ( k x , k y , − k z ) (E1)with C T := T C z represented by iσ x K . A symmetry-allowed Hamiltonian assumes the form H = ( Ak x + Bk y ) I + ( ak x + ck y + e ( k x + k y )) σ x + bk x σ y + dk z σ z . (E2)We do not exhaust all possible quadratic terms – oneis sufficient to demonstrate the anisotropy of λ . One mayverify that the quadratic term causes φ R + φ B to deviatefrom π (see App. C). Furthermore, k y σ y is not present inthe Hamiltonian since it can be eliminated by a rotationin the k x − k y plane. We choose the parameters to be A= 0.2, B = 0.2, a = 1, c = 0.3, b = 1, d = 0.2, e = 0.5(eV · ˚A) to produce Fig. 2 in the main text. In Fig. 2,constant energy contour is evaluated at -0.02eV.5 R. Peierls, Zeitschrift f¨ur Physik , 763 (1933). L. Onsager, The London, Edinburgh, and Dublin Philo-sophical Magazine and Journal of Science , 1006 (1952). L. M. Lifshitz and A. Kosevich, Dokl. Akad. Nauk SSSR , 963 (1954). W. J. de Haas and P. M. van Alphen, Proc. NetherlandsRoy. Acad. Sci. , 1106 (1930). L. W. Shubnikov and W. J. de Haas, Proc. NetherlandsRoy. Acad. Sci. , 130 (1930). D. Shoenberg,
Magnetic oscillations in metals (Cam-bridge University Press, 1984). N. W. Ashcroft and N. D. Mermin,
Solid state physics (Thomson Learning, 1976). T. Champel and V. P. Mineev, Philosophical MagazinePart B , 55 (2001). W. Kohn, Phys. Rev. , 1460 (1959). L. Roth, Journal of Physics and Chemistry of Solids ,433 (1962). E. I. Blount, Phys. Rev. , 1636 (1962). M. Wilkinson and R. J. Kay, Phys. Rev. Lett. , 1896(1996). M. V. Berry, Proc. R. Soc. Lond A , 45 (1984). G. P. Mikitik and Y. V. Sharlai, Phys. Rev. Lett. , 2147(1999). M.-C. Chang and Q. Niu, Phys. Rev. B , 7010 (1996). L. M. Roth, Phys. Rev. , 434 (1966). H. J. Fischbeck, Physica Status Solidi (b) , 11 (1970). Y. Gao and Q. Niu, Proceedings of the Na-tional Academy of Sciences G. Zil’berman, JETP , 208 (1957). S. Jeon, B. B. Zhou, A. Gyenis, B. E. Feldman, I. Kimchi,A. C. Potter, Q. D. Gibson, R. J. Cava, A. Vishwanath,and A. Yazdani, Nature Materials , 851 (2014). I. Zeljkovic, Y. Okada, M. Serbyn, R. Sankar, D. Walkup,W. Zhou, J. Liu, G. Chang, Y. J. Wang, M. Z. Hasan,F. Chou, H. Lin, A. Bansil, L. Fu, and V. Madhavan,Nature Materials , 318 (2015). I. A. Luk’yanchuk and Y. Kopelevich, Phys. Rev. Lett. , 166402 (2004). A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W.Ludwig, Phys. Rev. B , 195125 (2008). A. Kitaev, AIP Conf. Proc. , 22 (2009). B. Simon, Phys. Rev. Lett. , 2167 (1983). F. Wilczek and A. Zee, Phys. Rev. Lett. , 2111 (1984). G. H. Wannier and D. R. Fredkin, Phys. Rev. , 1910(1962). G. Nenciu, Rev. Mod. Phys. , 91 (1991). A. Alexandradinata and L. Glazman, ArXiv e-prints(2017), arXiv:1710.04215 [cond-mat.other]. M. H. Cohen and L. M. Falicov, Phys. Rev. Lett. , 231(1961). M. Y. Azbel, JETP , 891 (1961). J. B. Keller, Annals of Physics , 180 (1958). A. Alexandradinata and L. Glazman, ArXiv e-prints(2017), arXiv:1708.09387 [cond-mat.mes-hall]. Proceedings of the Royal Society of London A: Mathemat-ical, Physical and Engineering Sciences , 500 (1952),http://rspa.royalsocietypublishing.org/content/211/1107/500.full.pdf. R. B. Dingle,
Proceedings of the Royal Society of London.Series A. Mathematical and Physical Sciences , , 517 (1952). T. Champel, Phys. Rev. B , 054407 (2001). S. K. Kushwaha, J. W. Krizan, B. E. Feldman, A. Gye-nis, M. T. Randeria, J. Xiong, S.-Y. Xu, N. Alidoust,I. Belopolski, T. Liang, M. Z. Hasan, N. P. Ong, A. Yaz-dani, and R. J. Cava, APL Materials , 041504 (2015),https://doi.org/10.1063/1.4908158. I. Giaever and K. Megerle, Phys. Rev. , 1101 (1961). L. Sun, D.-J. Kim, Z. Fisk, and W. K. Park, Phys. Rev.B , 195129 (2017). C. J. Chen,
Introduction to scanning tunneling microscopy (Oxford University Press, 2008). J. Dhillon and D. Shoenberg, Philosophical Transactionsof the Royal Society of London A: Mathematical, Physicaland Engineering Sciences , 1 (1955). M. H. Cohen and E. I. Blount, The Philosoph-ical Magazine: A Journal of Theoretical Ex-perimental and Applied Physics , 115 (1960),http://dx.doi.org/10.1080/14786436008243294. Where L>
1, we consider L symmetry-related propagators A i with an additional index i ∈ Z L ; each of A i is a matrixof dimension D , and its eigen-phases are denoted by λ ia ,with a ∈ Z D . A. Alexandradinata, Z. Wang, and B. A. Bernevig, Phys.Rev. X , 021008 (2016). K. Shiozaki, M. Sato, and K. Gomi, ArXiv e-prints(2017), arXiv:1701.08725 [cond-mat.mes-hall]. M. Lax,
Symmetry principles in solid state and molecularphysics (Wiley-Interscience, 1974). A. Alexandradinata, X. Dai, and B. A. Bernevig, Phys.Rev. B , 155114 (2014). While it is known that φ B = π for a time-reversal-invariant2D Dirac fermion, it is not generally appreciated that φ R = φ Z =0. L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. ,106803 (2007). L. Fu and C. L. Kane, Phys. Rev. B , 045302 (2007). J. E. Moore and L. Balents, Phys. Rev. B , 121306(2007). R. Roy, Phys. Rev. B , 195322 (2009). T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, andL. Fu, Nat. Comm. , 982 (2012). H. B. Nielsen and M. Ninomiya, Nucl. Phys. B , 20(1981). E. McCann and V. I. Fal’ko, Phys. Rev. Lett. , 086805(2006). A. Varlet, M. Mucha-Kruczyski, D. Bischoff, P. Simonet,T. Taniguchi, K. Watanabe, V. Falko, T. Ihn, and K. En-sslin, Synthetic Metals , 19 (2015), reviewsof Current Advances in Graphene Science and Technology. X. Xu, W. Yao, D. Xiao, and T. F. Heinz, Nature Physics , 343 (2014). G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly,and J. van den Brink, Phys. Rev. B , 073103 (2007). S. Y. Zhou, G.-H. Gweon, A. V. Fedorov, P. N. First,W. A. de Heer, D.-H. Lee, F. Guinea, A. H. Castro Neto,and A. Lanzara, Nat Mater , 916 (2007). G. W. Semenoff, Phys. Rev. Lett. , 2449 (1984). D. Xiao, W. Yao, and Q. Niu, Phys. Rev. Lett. , 236809(2007). The approximate equality of φ B + φ R to π may be under- stood from a linearized, two-band model without next-nearest-neighbor hoppings . C. R. Woods, L. Britnell, A. Eckmann, R. S. Ma, J. C.Lu, H. M. Guo, X. Lin, G. L. Yu, Y. Cao, R. V. Gor-bachev, A. V. Kretinin, J. Park, L. A. Ponomarenko,M. I. Katsnelson, Y. N. Gornostyrev, K. Watanabe,T. Taniguchi, C. Casiraghi, H.-J. Gao, A. K. Geim, andK. S. Novoselov, Nat Phys , 451 (2014), article. L. Muechler, A. Alexandradinata, T. Neupert, andR. Car, Phys. Rev. X , 041069 (2016). J. H¨oller and A. Alexandradinata, ArXiv e-prints (2017),arXiv:1708.02943 [cond-mat.other]. A. A. Soluyanov, D. Gresch, Z. Wang, Q. Wu, M. Troyer,X. Dai, and B. A. Bernevig, Nature , 495 (2015). In principle, there may be a four-fold band touching (e.g.,3D Dirac point) at a T -invariant point that is acciden-tal or imposed by crystalline symmetries. Note that suchband touchings are not stable under a T -symmetric per-turbation that breaks every other symmetry except, re-spectively, spin SU(2) symmetry in class (a), and i sym-metry in class (b). The value of det A cannot change dueto this perturbation, since such value is quantized only by T symmetry. The result of this argument is that we mayalways simplify the dispersion at the T -invariant point toa band extremum, and apply the proof in the main text.This perturbation argument is developed more carefullyin Ref. 29. By a T -symmetric perturbation, one may always simplifythe dispersion to a Dirac point. This argument is analo-gous to that presented in the previous footnote. Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu,H. Weng, X. Dai, and Z. Fang, Phys. Rev. B , 195320(2012). J. Xiong, S. Kushwaha, J. Krizan, T. Liang, R. J. Cava,and N. P. Ong, EPL (Europhysics Letters) , 27002(2016). J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan,M. Hirschberger, W. Wang, R. Cava, and N. Ong, Science , 413 (2015). H. Lan-Po and L. Shi-Yan, Chinese Physics B , 117105(2016). Z. J. Xiang, D. Zhao, Z. Jin, C. Shang, L. K. Ma, G. J.Ye, B. Lei, T. Wu, Z. C. Xia, and X. H. Chen, Phys. Rev.Lett. , 226401 (2015). A. Pariari, P. Dutta, and P. Mandal, Phys. Rev. B ,155139 (2015). D. Kumar and A. Lakhani, Physica Status Solidi (RRL) , 636 (2015). A. Narayanan, M. D. Watson, S. F. Blake, N. Bruyant,L. Drigo, Y. L. Chen, D. Prabhakaran, B. Yan, C. Felser,T. Kong, P. C. Canfield, and A. I. Coldea, Phys. Rev.Lett. , 117201 (2015). Y. Zhao, H. Liu, C. Zhang, H. Wang, J. Wang, Z. Lin,Y. Xing, H. Lu, J. Liu, Y. Wang, S. M. Brombosz, Z. Xiao,S. Jia, X. C. Xie, and J. Wang, Phys. Rev. X , 031037(2015). If D =1, the three one-forms commute, and therefore λ may by decomposed. We remark that the π Berry phaseof a Weyl point originates from an integral of the abelian
Berry connection in a single-band subspace. For small orbits encircling a high-symmetry wavevector,this effective g-factor may alternatively be calculated inthe effective-mass approximation. For larger orbits, theeffective-mass approximation is inappropriate, but Eq. (2) remains valid. Our method is essentially identical to pre-vious works by Roth , Mikitik and Sharlai; the onlydifference is that a basis choice (gauge) has been chosen inprevious works such that their one-forms are traceless. In comparison, no gauge-fixing is assumed in Eq. (2), andtherefore Eq. (2) is the most natural object to calculatenumerically. As mentioned in Sec. III, this procedure requires an inde-pendent determination of | S zz | . Note this implies a double constraint on | µ | . M. Tinkham,
Group Theory and Quantum Mechanics (Dover, 2003). H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C.Zhang, Nat. Phys. , 438 (2009). C.-X. Liu, X.-L. Qi, H. Zhang, X. Dai, Z. Fang, and S.-C.Zhang, Phys. Rev. B , 045122 (2010). This robustness is an artifact of a restricted class of k · p models described in App. C. The robustness is not gen-erally valid at larger chemical potentials where the effectof other bands become significant. A similar observationwas first found in Ref. 98 in the context of graphene witha Semenoff mass. M. Brahlek, N. Bansal, N. Koirala, S.-Y. Xu, M. Neupane,C. Liu, M. Z. Hasan, and S. Oh, Phys. Rev. Lett. ,186403 (2012). W. Liu, X. Peng, C. Tang, L. Sun, K. Zhang, andJ. Zhong, Phys. Rev. B , 245105 (2011). S. M. Young, S. Chowdhury, E. J. Walter, E. J. Mele,C. L. Kane, and A. M. Rappe, Phys. Rev. B , 085106(2011). Y. Liu, Y. Y. Li, S. Rajput, D. Gilks, L. Lari, P. L.Galindo, M. Weinert, V. K. Lazarov, and L. Li, NaturePhysics , 294 (2014). E. Adams and T. Holstein, Journal of Physics and Chem-istry of Solids , 254 (1959). L. M. Roth and P. N. Argyres, in
Semiconductors andSemimetals , edited by R. K. Willardson and A. C. Beer(Academic Press, New York and London, 1966), Vol. 1. K. Eto, Z. Ren, A. A. Taskin, K. Segawa, and Y. Ando,Phys. Rev. B , 195309 (2010). A. M. Kosevich and V. V. Andreev, JETP , 637 (1960). E. Lifshitz and L. P. Pitaevskii,
Physical Kinetics (But-terworth Heinemann, 1979). In addition, there are auxiliary conditions such as (cid:126) ω c (cid:28) ( kT µ ) / if kT and (cid:126) /τ (cid:28) (cid:126) ω c . It has been assumedin Ref. 94 that impurity-induced scattering between spin-split Landau levels is absent; we are not certain that thisassumption is generally valid in spin-orbit-coupled metals. J. Zak, Phys. Rev. Lett. , 2747 (1989). D. Gresch, G. Aut`es, O. V. Yazyev, M. Troyer, D. Van-derbilt, B. A. Bernevig, and A. A. Soluyanov, Phys. Rev.B , 075146 (2017). J. N. Fuchs, F. Pi´echon, M. O. Goerbig, and G. Montam-baux, The European Physical Journal B , 351 (2010). A derivation can be found in the subsection titled ‘TheSommerfeld expansion’ in Ref. 108.
P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, Davide Ceresoli, G. L. Chiarotti, M. Co-coccioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris,G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis,Anton Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, Stefano Paolini, A. Pasquarello,L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentz- covitch, J. Phys.: Condens. Matter , 395502 (2009). M. Schlipf and F. Gygi, Computer Physics Communica-tions , 36 (2015).
J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza,D. Vanderbilt, and N. Marzari, Computer Physics Com-munications , 2309 (2014).
W. Reckien, F. Janetzko, M. F. Peintinger, and T. Bre-dow, J. Comput. Chem. , 2023 (2012). C. Wang, X. Liu, L. Kang, B.-L. Gu, Y. Xu, andW. Duan, Phys. Rev. B , 115147 (2017). Website of owl: https://mistguy.gitlab.io/owl/.
G. P. Mikitik and Y. V. Sharlai, Journal of Experimentaland Theoretical Physics , 747 (1998). D. P. Arovas,