Magnetizations and de Haas-van Alphen oscillations in massive Dirac fermions
aa r X i v : . [ c ond - m a t . m e s - h a ll ] J a n Magnetizations and de Haas-van Alphen oscillations in massive Dirac fermions
F. R. Pratama, ∗ M. Shoufie Ukhtary, and Riichiro Saito
Department of Physics, Tohoku University, Sendai 980-8578, Japan (Dated: January 6, 2021)We theoretically study magnetic field, temperature, and energy band-gap dependences of magne-tizations in the Dirac fermions. We use the zeta function regularization to obtain analytical expres-sions of thermodynamic potential, from which the magnetization of graphene for strong field/lowtemperature and weak field/high temperature limits are calculated. Further, we generalize the resultby considering the effects of impurity on orbital susceptibility of graphene. In particular, we showthat in the presence of impurity, the susceptibility follows a scaling law which can be approximatedby the Faddeeva function. In the case of the massive Dirac fermions, we show that a large band-gapgives a robust magnetization with respect to temperature and impurity. In the doped Dirac fermion,we discuss the dependence of the band-gap on the period and amplitude of the de Haas-van Alpheneffect.
I. INTRODUCTION
Historically, theoretical study on the magnetic prop-erties of graphene [1–3] can be traced back to a paperby McClure [4] in 1956. He showed that the diamag-netism in undoped graphene is largely contributed bycoalescence of the massless Dirac electrons at the va-lence bands to the zeroth Landau levels (LLs) at the K and K ′ valleys in the hexagonal Brillouin zone in thepresence of an external magnetic field. At the zerothLLs, the free energy increases with the increasing mag-netic field, thus graphene shows orbital diamagnetism.Saito and Kamimura [5] showed that the orbital param-agnetism appears in graphene intercalation compounds.Furthermore, Raoux et al. [6] demonstrated that the dia-magnetism in a two-dimensional (2D) honeycomb latticecan be tuned into paramagnetism by introducing an ad-ditional hopping parameter, where the Berry phase isvarying continuously between π and 0.A method to derive analytical formula for orbital sus-ceptibility of a massive Dirac system is developed byKoshino and Ando [7, 8]. In their method, the Euler-Maclaurin expansion formula is applied to calculate thethermodynamic potential in the presence of the mag-netic field. They showed that the pseudospin param-agnetism is responsible for a singular orbital suscepti-bility inside the band-gap region and disappears whenthe chemical potential enters the valence or conductionband [7–9]. The expansion formula was first used byLandau [10] to demonstrate the orbital diamagnetism inmetals. In the context of graphene-related materials, theEuler-Maclaurin expansion is employed to calculate or-bital susceptibilities of transition-metal dichalcogenides(TMDs) [11] and the Weyl semimetals [12]. However,the magnetization as a function of magnetic field B andtemperature T can not be obtained by using the Euler-Maclaurin expansion formula. It is because that the mag-netization diverges due to an infinite number of the LLs ∗ pratama@flex.phys.tohoku.ac.jp formed in the valence bands that areincluded in the cal-culation of the thermodynamic potential, unless a cut-offof the LLs is introduced [13]. The calculated magnetiza-tion M as a function of B has a form M ∝ C − C B ( C > C becomes infinite withincreasing the number of the LLs, while in the case ofa conventional metal, only LLs in the conduction bandsare considered. Moreover, the Euler-Maclaurin expan-sion method is valid only when the spacing of the LLs ismuch smaller than the thermal energy k B T , where k B isthe Boltzmann constant. [10].In the recent study by Li et al. [14], the magnetiza-tion of undoped graphene is measured for a wide rangeof B and temperature T . In the strong B /low T limit,it is shown that the magnetization is proportional to thesquare-root of the magnetic field and diminishes linearlywith increasing temperature ( M ∝ −√ B +constant × T ),while in the weak B /high T limit, it is observed that themagnetization is proportional to B and inversely propor-tional to T ( M ∝ − B/T ). The experimental data andnumerical calculations are fitted into a Langevin func-tion, from which the properties of magnetization for thetwo limiting cases can be deduced [14].To avoid the divergence in the magnetization, we de-rive analytical formula for thermodynamic potentials ofthe Dirac systems by using the zeta function regular-ization. In the context of quantum field theory, thismethod was used by Cangemi and Dunne [15] to calcu-late the energy of relativistic fermions in magnetic field.In graphene-related topics, the zeta function regulariza-tion was employed by Ghosal et al. [16] to explain theanomalous orbital diamagnetism of graphene at T = 0 Kand by Slizovskiy and Betouras [17] to show the non-linear magnetization of graphene in a strong B at high T . In this paper, we derive the analytical formula formagnetization of graphene for both strong B /low T andweak B /high T limits by the zeta function regularization,which reproduces the Langevin fitting to the experimen-tal observation. Further, we discuss the effect of impurityon the orbital susceptibility of graphene. We show thatin the presence of impurity, the orbital susceptibility asa function of temperature follows a scaling law which isapproximately given by the so-called Faddeeva function.The effects of energy band-gap on the magnetization inundoped and doped Dirac systems are also discussed. Inthe undoped case, large band-gaps in TMDs yield rela-tively small but robust magnetizations with respect totemperature and impurity. In the doped case, we showthat the opening of the band-gap is observed from the di-minishing amplitude of the de Haas-van Alphen (dHvA)oscillation at T = 0 K. This phenomenon can not beobtained by the Euler-Maclaurin expansion method [4].The paper is organized as follows. In Sec. II, wepresent analytical methods for calculating the LLs, ther-modynamic potential, and magnetization. In Sec. III,calculated results are discussed. In Sec. IV, conclusionis given. II. CALCULATION METHODSA. The Landau levels of massive Dirac fermions
As a starting point, we consider a massive Dirac systemwith a band gap of ∆ >
0. We employ a 2 × λ in the Hamilto-nian [18]. The energy dispersions are approximated bya linear function of momentum p = ( p x , p y ), and theZeeman term is neglected because we only consider twocases: (1) at low temperature, the Zeeman splitting ismuch smaller than the Landau levels (LLs) separationand (2) at high temperature for TMDs. In the presenceof an external magnetic field B = B ˆ z , the momentum ac-quires an additional term by the Peierl substitution, i.e. p → p + e A . The vector potential A is related to B by B = ∇× A . By choosing the Landau gauge A = (0 , Bx ),the Hamiltonian is given by [2, 3, 19]:ˆ H τs = (cid:20) ∆ / v F { τ p x − i ( p y + eBx ) } v F { τ p x + i ( p y + eBx ) } − ∆ / λτ s (cid:21) . (1)Here, v F is the Fermi velocity whose typical value for theDirac fermions is ∼ m / s, τ = +1 ( −
1) is the indexfor the K ( K ′ ) valley and s = +1 ( −
1) is the index forspin-up (spin-down).To solve Eq. (1), we define annihilation and creationoperators ˆ a ≡ [ ℓ B / ( √ ~ )][ ip x + ( p y + eBx )] and ˆ a † ≡ [ ℓ B / ( √ ~ )][ − ip x + ( p y + eBx )], where ℓ B = p ~ / ( eB )is the magnetic length. In term of the annihilation andcreation operators, the Hamiltonian for a given τ and s reduces to ˆ H ξ = (cid:20) ∆ / − i ~ ω c ˆ O τ i ~ ω c ˆ O † τ − ∆ / λξ (cid:21) , (2)where ξ ≡ τ s , ˆ O + ≡ ˆ a , ˆ O − ≡ ˆ a † , and ω c = √ v F /ℓ B = p v F eB/ ~ is the cyclotron frequency of the Dirac (cid:1) n [ e V ] (cid:2) (cid:3)(cid:4) (c) (cid:0) n [ e V ] (cid:5) K n=0 (a) (cid:6) (cid:7) (cid:8) (cid:9) (cid:10) K' n=0 (b)K (cid:11) n [ e V ] (cid:12) (e) K n=0n=0 (cid:13)(cid:14) n=0 (cid:15)(cid:16) (d) K'(f) K' (cid:17)(cid:18) n= (cid:19) (cid:20) (cid:21) (cid:22) (cid:23) (cid:24) FIG. 1. (Color online) The LLs ( n = − n = 5) of themassive Dirac systems. The figures in the left (right) sidecorrespond to the LLs at the K ( K ′ ) valley. (a), (b): LLs ofa gapped graphene (∆ = 100 meV, λ = 0, and v F = 10 m / s)for B = 0 − (∆ = 1 .
66 eV, λ = 75 meV and v F = 5 . × m / s) for B = 0 −
20 T. TheLLs for the spin-up and spin-down electrons are shown withred solid and blue dashed lines, respectively. fermions. ξ = +1( −
1) represents the spin-up (spin-down) electron at the K valley or spin-down (spin-up)electron at K ′ valley. The n -th LLs ǫ ξn and wave func-tion | Ψ ξn i are given by the eigenvalues and eigenvectorsof Eq. (2), respectively, as follows: ǫ ξn = ξλ τ ( n ) r ( ~ ω c ) | n | + (cid:16) ∆ ξ (cid:17) (3)and | Ψ ξn i = 1 q | ǫ ξn − ξλ | − i q | ǫ ξn + ∆2 − λξ || α τn i q | ǫ ξn − ∆2 || β τn i , (4)where we define ∆ ξ ≡ ∆ − ξλ as a shorthand notation.sgn τ ( n ) is the sign function defined by sgn τ ( n ) = − n < τ ( n ) = +1 for n >
0. For n =0, a non-trivial wave function is satisfied by choosingsgn + (0) = − − (0) = +1. The eigenvectors | α τn i and | β τn i are opposite for the K and K ′ valleys,i.e. | α + n i = | β − n i ≡ || n | − i and | α − n i = | β + n i ≡ || n |i ,where || n |i is a normalized eigenvector of the opera-tors ˆ a † and ˆ a † , such that ˆ a || n |i = p | n ||| n | − i andˆ a † || n |i = p | n | + 1 || n | + 1 i . The zeroth LLs at the K and K ′ valleys exist at the valence and the conduction bands,respectively [7, 8]. The presence of only one zeroth LLin each valley is confirmed by first-principle calculationsfor hexagonal boron nitride (h-BN) and MoS [20]. For∆ < λ = 0, a non-trivial wave function for n = 0is satisfied by sgn + (0) = +1 and sgn − (0) = −
1, as inthe case of topological silicene [21]. Nevertheless, in thisstudy we only consider ∆ > K and K ′ valleys is same as used in references [7, 8, 22]which is opposite of those references [11, 23].In Fig. 1, we plot the LLs ( n = − n = 5) of agapped graphene [(a) and (b)] and MoS [(c)-(f)] at the K and the K ′ valleys as a function of the magnetic field B . In (a) and (b), the LLs of the gapped graphene with∆ = 100 meV, λ = 0, and v F = 10 m / s show √ B dependences, because ∆ / ~ ω c (72 . B = 4 T). In Fig. 1 (c)-(f), theLLs of MoS at the conduction bands [(c) and (d)] andthe valence bands [(e) and (f)] are shown, where we adopt∆ = 1 .
66 eV, λ = 75 meV and v F = 5 . × m / s [18, 22–24]. We can see that the SOC generates spin-splittingbetween the spin-up (red solid lines) and the spin-down(blue dashed lines) electrons except for the zeroth LLsat the K ′ valley [Fig. 1(d)]. For the valence band, aspin-splitting 2 λ = 150 meV occurs for the zeroth LLs atthe K valley [Fig. 1(e)]. The LLs are linearly dependentfor B = 0 −
20 T because ∆ / ~ ω c = 0 .
086 eV for B = 20 T. B. Thermodynamic potential and magnetization
The thermodynamic potential per unit area at tem-perature T in the presence of the magnetic field is givenby Ω = − β eBh X ξ = ± ∞ X n = −∞ ln[1 + e − β ( ǫ ξn − µ ) ] ≡ X ξ = ± (Ω − + Ω + ) , (5)where β = 1 / ( k B T ) and µ is the chemical potential. Thepre-factor eB/h represents the Landau degeneracy perunit area. For expository purposes, we define Ω − andΩ + as the potentials for the LLs at the valence and con-duction bands, respectively. More specifically, Ω − is thethermodynamic potential for the n < n = 0 LL at the K valley, and Ω + is the thermodynamic potentialfor the n > n = 0 LL at the K ′ valley. Theupper limit of n depends on the magnitude of the LLsspacing compared with the thermal energy k B T . Whenthe LLs spacing is much larger than k B T , the summationonly includes the occupied LLs. When the LLs spacingis much smaller than k B T , the summation can be takenfor all LLs. In order to avoid divergence in Eq. (5), theinfinite summations of the LLs are expressed by the zetafunction, from which a finite value from an infinite sum-mation can be obtained through the method of analyticalcontinuation.In the absence of SOC ( λ = 0), we drop the summationby the index ξ because all the LLs are degenerate for thespin-up and spin-down electrons, and the thermodynamicpotential is multiplied by g s = 2 to account the spindegeneracy. After obtaining an analytical expression ofΩ, the magnetization is calculated by M = − ∂ Ω ∂B . (6)In the presence of impurity, the magnetization for agiven scattering rate γ is calculated by convolution of M in Eq. (6) with a Lorentzian profile as follows [25–29]: M ( µ, γ ) = γπ Z ∞−∞ dεM ( ε ) 1( ε − µ ) + γ . (7)The parameter γ is related to the self-energy due im-purity scattering, and γ is inversely proportional to therelaxation time of the quasiparticle. For simplicity, weassume that γ is independent of B and T , and thereforethe susceptibility as a function of temperature is givenby χ ( µ, γ, T ) = [ ∂M ( µ, γ, T ) /∂B ] B =0 .
1. Thermodynamic potential for ~ ω c ≫ k B T, ∆ ξ ≫ k B T First, let us derive the thermodynamic potential for ~ ω c ≫ k B T . Since we consider electron-doped system,we get µ > ∆ /
2. The logarithmic function in Eq. (5) isapproximated by ln[1 + exp {− β ( ǫ ξn − µ ) } ] ≈ − β ( ǫ ξn − µ )which is valid in the case of − β ( ǫ ξn − µ ) ≫ k B T or T → n
0, Ω − is expressed by (see Appendix A 1 forderivation),Ω − = − eBh X ξ = ± " ~ ω c ζ (cid:16) − , Γ ξ (cid:17) − ∆ ξ , (8)where we define Γ ξ ≡ ∆ ξ / (2 ~ ω c ), and the infinite sum-mation of the LLs is given by the Hurwitz zeta functionwhich is defined by (see for example reference [30]), ζ ( p, q ) ≡ ∞ X k =0 k + q ) p . (9)It is noted that the chemical potential µ does not ap-pear in the expression of Ω − for the electron-doped sys-tem. Eq. (8) gives the intrinsic diamagnetism of TMDsat T = 0 K. A similar result was derived by Sharapovet al. [25] for a gapped graphene, which is obtained byintroducing an ultraviolet cut-off in the calculation of thethermodynamic potential. By using the asymptotic formof ζ ( p, q ) for q → ∞ and p = 1, ζ ( p, q ) ∼ q − p p − q − p pq − (1+ p ) , (10)the Ω − of a heavy Dirac fermion (∆ ξ / ≫ ~ ω c ) is givenby Ω − ≈ X ξ = ± π ∆ ξ ( ~ v F ) + eBh ( ~ ω c ) ξ ! . (11)In the calculation of Ω ( e )+ , we define an integer ν ξ asthe index of the highest occupied LLs in the conductionbands since the summations of the LLs in the thermody-namic potential are carried out up to the ν ξ -th levels, asfollows: ν ξ ≡ $ ( µ − ξλ/ − (∆ ξ / ( ~ ω c ) % ≡ ⌊ ˜ ν ξ ⌋ , (12)where the floor function ⌊ x ⌋ is defined by the greatestinteger smaller than or equal to x . It is noted that when λ = 0, we drop subscript ξ in Eq. (12). When we in-troduce the step function Θ( µ − ∆ /
2) as a threshold, weconfirm that µ > ∆ / ( e )+ , as follows (see Appendix A 2 for derivation):Ω ( e )+ = − eBh X ξ = ± " ∆4 + µ (cid:16) ν ξ + 12 (cid:17) − ξλ ν ξ + 1)+ ~ ω c ( ζ (cid:16) − , Γ ξ + ν ξ + 1 (cid:17) − ζ (cid:16) − , Γ ξ (cid:17)) × Θ( µ − ∆ / . (13)In Eq. (13), we need the fact that the finite summationof the LLs is expressed by a subtraction of two zeta func-tions as follows [30]: N X n =0 ( n + q ) − p = ζ ( p, q ) − ζ ( p, q + N + 1) . (14)In Appendix B, we show that the numerical calculationof the left-hand side of Eq. (14) reproduces the analyti-cal expression on the right-hand side. Therefore, for theelectron-doped case, the total thermodynamic potentialis given by the summation Ω − + Ω ( e )+ [Eqs. (8) and (13)].In the case of hole-doped system, on the other hand ν ξ is now defined as the highest occupied LLs in the valence bands. By using the same procedure as Ω ( e )+ , Ω ( h ) − is givenas follows:Ω ( h ) − = 2 eBh X ξ = ± " ∆4 + | µ | (cid:16) ν ξ + 12 (cid:17) + ξλ ν ξ + ~ ω c ( ζ (cid:16) − , Γ ξ + ν ξ + 1 (cid:17) − ζ (cid:16) − , Γ ξ (cid:17)) × Θ( | µ | − ∆ ξ / . (15)The total thermodynamic potential is given by the sub-traction Ω − − Ω ( h ) − [Eqs.(8) and (15)] because Ω ( h ) − repre-sents the ”potential” of the unoccupied LLs. It is notedthat the electron- and hole-doped Dirac systems giveidentical thermodynamic potentials in the case of λ = 0,because of the electron-hole symmetry.
2. Thermodynamic potential for ~ ω c ∼ ∆ ≪ k B T Now, let us derive the thermodynamic potential of agapped graphene ( ~ ω c ∼ ∆, λ = 0) at high temper-ature. By considering µ > ∆ /
2, the total thermody-namic potential is given by Ω = Ω − + Ω ( e )+ . Here, weexpress Ω ( e )+ = Ω + − Ω ′ + , where Ω + and Ω ′ + are thethermodynamic potentials from the entire LLs and fromLLs n = ν to n = ∞ at the conduction bands, respec-tively. In Appendix A 3, we show that Ω ′ + is negligible for µ ≪ k B T . We expand the logarithmic and exponentialterms in the first line of Eq. (5) to calculate Ω ∓ which isvalid for [ − ǫ ξ ∓ n + µ ] ≪ k B T . The thermodynamic poten-tial Ω = Ω − + Ω − is given by a power series of ( ~ ω c ) asfollows:Ω = 4 g s β eBh ∞ X ℓ =0 Li − ℓ ( − e βµ ) ( β ~ ω c ) ℓ (2 ℓ )! " ζ ( − ℓ, Γ ) − Γ ℓ ≡ ∞ X ℓ =0 Ω ℓ , (16)where Γ = ∆ / (2 ~ ω c ) for clarity. It is noted that the de-pendence of Ω on µ is given by the polylogarithm func-tion Li s ( z ) ≡ P ∞ k =1 z k /k s , which converges for | x | > ζ ( − ℓ, x ) = −B ℓ +1 ( x ) / ( ℓ + 1) [30], where B ℓ ( x ) isthe Bernoulli polynomial, B ( x ) = x − / ( z ) = − ln(1 − z ) [30], Ω is given byΩ = k B T g s π ∆ ( ~ v F ) ln[1 + e βµ ] . (17)Ω is proportional to the square of band-gap ∆ as wellas temperature T and does not depend on the magneticfield. Thus the Ω can be interpreted as the amount M [ - A ] -5 -4-3-2 M [ - A ] M S [ - A ] T=0 K25 K 50 K75 K100 K10 T8 T6 T4 TB=2 TB=2 T4 T6 T 8 T10 T (a)(b)(c)
FIG. 2. (Color online) Magnetization of graphene at ~ ω c ≫ k B T limit as a function of magnetic field and temperature. (a) M as a function of B = 0 −
10 T at fixed temperatures, (b) M and (c) M S as a function of T ( T = 0 −
100 K) for severalvalues of B . The calculations from the analytical formula[Eq. (23)] and the Langevin function [Eq. (24)] are depictedby symbols and lines, respectively. of thermal energy required to excite electrons from thevalence bands across the band-gap. Ω ℓ with ℓ = 0 gives alinear response of the magnetization as a function of B ,as shall be discussed in the next section. III. RESULTSA. Magnetization of graphene ( ∆ = 0 , λ = 0 ) In the case of graphene, the LL is given by ǫ ξn =sgn τ ( n ) = ~ ω c √ n and the zeroth LLs ǫ ξ = 0 are sharedbetween the valence and conduction bands at the K and K ′ valleys, respectively. By separating the zeroth LLs M [ - A ] T=200 K225 K275 K300 K250 K -0.25-0.2-0.15-0.1-0.05 M [ - A ] -0.5-0.45-0.4-0.35 (a)(b)(c) -1 [K -1 ] M / B [ - A T - ] FIG. 3. (Color online) Magnetization of graphene at ~ ω c ≪ k B T limit as a function of magnetic field and temperature.(a) M as a function of B = 0 − . M and (c) M/B as a function 1 /T ( T = 200 −
300 K)for several values of B . The calculations from the analyticalformula [Eq. (23)] and the Langevin function [Eq. (24)] aredepicted by symbols and lines, respectively. from the n = 0 LLs, the thermodynamic potentials Ω − and Ω + for ~ ω c ≫ k B T are given byΩ − = − g s β eBh " ln { e βµ } + 2 ∞ X n =1 ln { e β ( ~ ω √ n + µ ) } = − g s eBh " β ln { e βµ } + ~ ω c ζ (cid:16) − , (cid:17) − µ , (18)andΩ + = − g s β eBh " ln { e βµ } + 2 ν X n =1 ln { e − β ( ~ ω √ n + µ ) } = − g s eBh " β ln { e βµ } + µν + ~ ω c ( ζ (cid:16) − , ν (cid:17) − ζ (cid:16) − , (cid:17)) , (19)respectively. We can confirm that by putting ∆ = λ = 0,and T → µ = 0), thetotal thermodynamic potential is therefore given byΩ( µ = 0) = − g s eBh " ~ ω c ζ − , ! + 1 β ln(2) ≡ Ω B + Ω S . (20)Since the thermodynamic potential can be equivalentlyexpressed by Ω = E − T S , where E and S are, respec-tively, the internal energy and entropy, we define Ω B andΩ S in Eq. (20) to denote the thermodynamic potentialassociated with the n < S of thezeroth LLs, respectively.In the case of ~ ω c ≪ k B T , the magnetization is givenby Ω in Eq. (16) by putting ∆ = 0 as follows:Ω = − g s β eBh ∞ X ℓ =1 Li − ℓ ( − e βµ ) ( β ~ ω c ) ℓ (2 ℓ )! B ℓ +1 ℓ + 1 . (21)Since Li − ( z ) = z/ (1 − z ) , we reproduce the formula forsusceptibility of graphene derived by McClure [4] χ = − e v F πk B T sech (cid:16) µ k B T (cid:17) . (22)It is noted that Eq. (22) is valid for any temperature T > B = 0 to calculate χ , thus thecondition ~ ω c ≪ k B T is always satisfied.From Eqs. (20) and (21), the magnetization of undopedgraphene ( µ = 0) is given by M = − . π e / v F ~ / √ B + 2ln(2) π e ~ k B T, ( ~ ω c ≫ k B T ) , − . e v F π Bk B T + O ( B ) , ( ~ ω c ≪ k B T ) . (23)It is noted that in Eq. (23) for ~ ω c ≪ k B T , only oddpowers of B survive because the Bernoulli number of B ℓ +1 is zero for even ℓ > B and T are fitted into a Langevin function L ( x ) = coth( x ) − /x as follows: M = − . π e / v F √ B ~ / L p ~ v F eBα ( T ) √ k B T ! . (24)The temperature dependence of M is approximated bythe function α ( T ) ≡ C/ ( C + √ T ), where C = 45 K / .Since L ( x ) ∼ x/ x → x → ∞ ,the magnetization is given by [14]: M ≈ − . π e / v F ~ / √ B + 0 . √ π e ~ k B T, ( ~ ω c ≫ k B T ) , − . √ e v F π Bk B T , ( ~ ω c ≪ k B T ) . (25)By comparing Eq. (23) with Eq. (25), the analytical for-mula reproduces experimentally observed B and T de-pendences of the magnetization of graphene both for ~ ω c ≫ k B T and ~ ω c ≪ k B T . It suggest that the zetaregularization works reasonably.In Fig. 2(a), we plot M ( B, T ) for ~ ω c ≫ k B T as afunction of B for several values of T by the analyticalexpression Eq. (23) (symbols) and the Langevin functionEq. (24) (lines). We can see that Eq. (23) works wellat temperature as high as T = 100 K for B > B < ~ ω c ≫ k B T does not satisfy for B < M as a function of T for several values of B .The linear dependence of M on T originates from theentropy S of electrons which coalesce to the zeroth LLs[see Eq. (20)]. In Fig. 2(c), we plot M S ≡ M + C √ B ,with C ≡ . e / v F / ( π ~ / ) [see Eq. (23)] as a functionof T for several values of B , which is a ’paramagnetic’contribution from the entropy of the zeroth LLs. Here,the functions M S are aligned into a straight line, whichshows that M S is independent of B .In Fig. 3(a), we plot M ( B, T ) for ~ ω c ≪ k B T as afunction of B for several values of T by the analyticalexpression Eq. (23) (symbols) and the Langevin functionEq. (24) (lines), where the linear B dependences of M for B ≤ . T = 200 K, especially for weak B ∼ . B , Eqs. (23) and (24) begin to show some discrepancies.In Fig. 3(b), M is plotted as a function of 1 /T for severalvalues of B . In Fig. 3(c), the function M/B is plotted asa function of 1 /T and is aligned into a straight line whichillustrates M ∝ − B/T dependence. From Eq. (16), it isinferred that the linear response of M with increasing B is originated from the contribution of the entire LLs atthe valence and conduction bands. -4.5-3-1.50 M [ - A ] (cid:25) =5 meV (cid:26) =40 meV (cid:27)=(cid:28)(cid:29)(cid:30) (cid:31) ! =1 eV T= 0 K FIG. 4. (Color online) Magnetization of massive Diracfermions (∆ = 5 meV ,
40 meV , . B = 0 −
10 T at T = 0 K. -15 Ω [ - J m -" ] -40 =100 meV, , =0 meV T -1 , < =0 meV > =100 meV, ?@ABCDE FGHIJK LMN OPQ RST UVW B -1 [T -1 ]-4 XY Z M [ - A ] [\] ^ _ FIG. 5. (Color online) The dHvA oscillations of the Diracfermions for µ = 100 meV , ∆ = 0 meV (solid lines), µ =70 . , ∆ = 0 meV (dashed lines), and µ = 100 meV , ∆ =141 . B. Magnetization for massive Dirac fermions T = 0 K , λ = 0 In Fig. 4, we plot the magnetization of the massiveDirac fermions as a function of B for several values of∆ at T = 0 K. We can see that the magnetization un- dergoes a gradual change from M ∝ −√ B to M ∝ − B dependences with increasing ∆ for B = 1 −
10 T, whichindicates that the anomalous orbital diamagnetism for∆ = 0 disappears with opening the gap. In this case, thespacing of the LLs which is initially p | n | dependence be-comes constant ( ~ ω c ) / ∆ with increasing ∆. This pro-cess can be observed by the transition from the topo-logical to the trivial phases of undoped silicene in whichthe band-gap can be controlled by applying an externalelectric field [31] perpendicular to the silicene plane. Asimilar transition is predicted in graphene with increas-ing temperature for the same reason, in which the p | n | dependence of the LLs in the valence bands is responsiblefor the M ∝ −√ B dependence. When the thermal en-ergy becomes larger than the cyclotron energy, the effect p | n | spacing of the LLs on the magnetization becomesno more important and thus the Dirac system shows lin-ear response M ∝ − B at a high T .The oscillation of magnetization or susceptibility in auniform magnetic field, which is known as the de Haas-van Alphen (dHvA) effect, has been observed experimen-tally in a quasi-2D system of graphite [32]. Much of thetheoretical studies on the magnetic oscillations in the 2Dsystems [33–36] have been carried out within the frame-work of a generalized Lifshits-Kosevich (LK) [37] theory,which was proposed to account the magnetic oscillationsin metals. In the LK theory, the dHvA effect is expressedby adding an oscillatory term to the Euler-Maclaurin for-mula (also known as the Poisson summation formula) forcalculating the thermodynamic potential [37].By assuming a fixed µ and λ = 0 in Eqs. (8) and (13),let us discuss the effect of band-gap on the period andamplitude of the dHvA oscillation at T = 0 K. In thebottom panel of Fig. 5(a), we plot Ω for µ = 100 meV and∆ = 0 meV as a function of inverse magnetic field 1 /B .At several values of 1 /B (labelled as 1 /B ν , ν = 1 , , ... ),we observe peaks of Ω which indicate the local maximaof potential and the peaks are separated by a period of0 .
13 T − . At 1 /B ν , the ν -th LLs at the K and K ′ valleysexactly match the chemical potential µ and thus we get ν = µ − (∆ / ~ v F eB ν = ~ A F πeB ν , (26)where A F = πk F = π [ µ − (∆ / ] / ( ~ v F ) is the areaof the Fermi surface of the Dirac system. The rightmostside of Eq. (26) is the relation derived by Onsager [38]to demonstrate that the dHvA oscillation can be utilizedto reconstruct the Fermi surface in metals. The periodof the dHvA oscillation in the massive Dirac system isgiven as follows [4, 25]: P = 1 B ν +1 − B ν = 2 ~ v F eµ − (∆ / . (27)In the middle and upper panels of Fig. 5(a), we plot Ω byadopting µ = 100 / √ ≈ . , ∆ = 0 meV and µ = 100 meV , ∆ = 100 √ ≈ . µ relative to the band gap ∆. This method is originallyproposed by Sharapov et al. [25] to detect the openingof band-gap in graphene with keeping µ constant. Ex-perimentally, the band-gap opening was observed [39]in epitaxially grown graphene on SiC substrate, where∆ ≈ .
26 eV is observed by breaking of sublattice sym-metry due to the graphene-substrate interaction.In Fig. 5(b) we plot the magnetization for the corre-sponding values of µ and ∆ provided in Fig. 5(a), wherethe oscillations exhibits a sawtooth-like feature. It isknown from the LK theory that the sawtooth-like os-cillations is a characteristic of the dHvA effect in 2Dsystems [25, 33, 40, 41]. We show that the sawtoothoscillation of M can be derived from the zeta func-tions in Eqs. Eq. (8) and (13). By using the formula ∂ζ ( p, q ) /∂q = − pζ ( p + 1 , q ), the magnetization is ex-pressed analytically as follows: M = 4 eh " ~ ωζ − , φ ! + µ ν + 12 ! − µ e ν X ν δ ( e ν − ν ) − ~ ωζ , φ !( Γ + e ν X ν δ ( e ν − ν ) ) , (28)where we define φ ≡ Γ + ν + 1. Thus, the sawtooth-like oscillation in the magnetization originates from thedelta function at ˜ ν = ν in Eq. (28), which is the resultof differentiation of the floor function in the expressionof ν [ ∂ ⌊ x ⌋ /∂x = P n ∈ Z δ ( x − n )]. Physically, the deltafunction indicates the occupations of electrons occupyingdiscrete LLs. With the increase of temperature, impurityscattering, electron-electron interactions, and electron-phonon interactions [25, 42–44], the LLs become broadin which the delta function is replaced by a Lorentzianfunction to account the broadening by the interactions.As a result, the oscillation of magnetization becomes lesssharp. The effects of the broadening on the dHvA os-cillation can be incorporated by the convolution of thethermodynamic potential at T = 0 K with the distribu-tion functions for temperature and impurities, as givenin references [25, 45, 46].We observe that in the cases of ∆ = 0 (solid and dashedlines in Fig. 5(b)), the smaller µ not only yields a decreas-ing frequency but also a weaker amplitude in the oscil-lation. When we consider the cases of p µ − (∆ / =70 . / (2 ~ ω )], and therefore the opening ofthe band-gap decreases the amplitude of the magnetiza-tion as the functions ζ ( − / , φ ) and ζ (1 / , φ ) possess thesame signs for a given φ (see Appendix B). The effect ofSOC on the magnetic oscillation in TMDs with the Zee- M [ - A ] T= 200 K300 K ‘ = 0 eV ab cde f g = 1 eV hijk FIG. 6. (Color online) Magnetization of undoped and doped( µ = 1 eV) MoS as a function B = 0 −
10 T at T = 200 and300 K. man splitting will be the subject of the next work. C. Magnetization of TMDs ( ∆ = 0 , λ = 0 ) For TMDs in a magnetic field B up to ∼
10 T, theLLs that are given by Eq. (3) are approximated by ǫ ξn ≈ ξλ − ∆ / − ( ~ ω c ) | n | / ∆ ξ and ǫ ξn ≈ ∆ / ~ ω c ) | n | / ∆ ξ for sgn τ ( n ) = − τ ( n ) = +1, respectively. Thus,the LLs separation is inversely proportional to the band-gap, i.e. ( ~ ω c ) / ∆ ξ . This approximation is also valid forheavy Dirac fermions such as h-BN where ∆ ≈ λ = 0. The thermodynamic potentials forTMDs in the case of ( ~ ω c ) / ∆ ξ ≪ k B T are given byΩ − = − β eBh X ξ = ± ∞ X l =1 Li − l [ − e β ( µ − ξλ +∆ / ] (cid:16) β ∆ ξ (cid:17) l × ( ~ ω c ) l l ! B l +1 l + 1 , (29)andΩ + = − β eBh X ξ = ± ∞ X l =1 Li − l [ − e β ( µ − ∆ / ] (cid:16) − β ∆ ξ (cid:17) l × ( ~ ω c ) l l ! B l +1 l + 1 . (30)Note that in Eqs. (29) and (30), the summations be-gin from l = 1 because the terms of l = 0 cancel withthe thermodynamic potentials originated from the zerothLLs n = 0 (see Appendix A 4 for derivation). Thus, theentropy of electrons at the zeroth LLs is not manifestedin a linear T dependence of as in the case of graphene.By keeping only the first term of Ω, the thermodynamicpotential of TMDs are given byΩ ≈ e v F B π X ξ = ± ξ sinh h β ∆ ξ i cosh h β ∆ ξ i + cosh h β (cid:16) µ − ξλ (cid:17)i . (31)From Eq. (31), we can see that Ω ∝ B and thereforethe magnetization is linearly proportional to B , whichprevails only for heavy Dirac fermions.In Fig. 6, we plot M of undoped MoS as a function of B for B = 0 −
10 T at T = 200 K and 300 K, where themagnetization does not change with the increasing tem-perature from T = 200 K to 300 K. Thus, even thoughthe magnitude of magnetization in a heavy Dirac fermiondecreases with the increasing band gap, the magnetiza-tion is robust for temperature. In fact, by comparingEq. (31) with (11), we infer that the magnetizations ofthe heavy Dirac fermions with a given ∆ at T = 0 K andfinite temperatures are equal, provided that ∆ / ≫ k B T .In Fig. 6, we also plot M for a doped case ( µ = 1 eV)at T = 200 K where the magnetization becomes zero,which demonstrates the effect of pseudospin paramag-netism [7, 8]. By neglecting the spin-orbit coupling, thesusceptibility of the massive Dirac fermion as a functionof ∆ and T is given by χ = − e v F π ∆ sinh h β ∆2 i cosh h β ∆2 i + cosh[ βµ ] , (32)which reproduce the result by Koshino and Ando [7, 8]with the Euler-Maclaurin formula. D. Susceptibility of grapehene and TMDs withimpurity
Finally, we analyse the effect of impurity scattering onthe orbital susceptibility of the Dirac fermions by usingEq. (7) for susceptibility. In the case of graphene, weapproximate the function sech ( βε/
2) in Eq. (22) by aGaussian function exp[ − ( Cβε ) ], where C is a constantdefined by C ≡ √ ln2 / [ √ √ ≈ .
447 (see Ap-pendix C for detail). The solution of the Voigt profile(convolution of the Gaussian with the Lorentzian func-tions) is given by the real part of the Faddeeva function w ( z ) as follows [30]: V ( x, y, σ ) ≡ yπ Z ∞−∞ dt exp[ − t / (2 σ )]( t − x ) + y = Re[ w ( z )] , (33)where σ is the standard deviation of the Gaussian func-tion, z = ( x + iy ) / ( √ σ ), and w ( z ) is the Faddeeva func-tion defined by w ( z ) ≡ e − z i √ π Z z dte t ! . (34) χ [ - A T - ] l mnopq r stuvwxyz
10 5 10 5 0 5 {|
10 0 10 }~(cid:127) (cid:128)(cid:129)(cid:130) (cid:131)(cid:132)(cid:133) (cid:134) B T/ γ (cid:135)(cid:136)(cid:137) (cid:138) / (cid:139) =0 (cid:140) (cid:141)(cid:142)(cid:143) χ / χ FIG. 7. (Color online) The susceptibility of graphene withimpurity as a function of temperature. (a) The calculation of χ as a function of T = 0 −
300 K for several values of γ and µ by the Faddeeva function (solid lines) numerical calculations(symbols). (b) The scaled susceptibility χ/χ as a functionof k B T /γ for µ/γ = 0 , , and 4 with the Faddeeva function. χ [ (cid:144)(cid:145) A T - ] (cid:146) (cid:147) (cid:148)(cid:149) (cid:150)(cid:151)(cid:152)
15 meV10 meV5 meV0 meV
FIG. 8. (Color online) The susceptibility of undoped MoS with impurity as a function of temperature T = 0 −
300 K for γ = 0 , ,
10, and 20 meV.
Therefore, in the presence impurity, the orbital suscepti-bility of graphene is approximately given by χ ( µ, γ ) ≈ − e v F πk B T Re[ w ( z ′ )] , (35)where we define z ′ ≡ Cβ ( µ + iγ ).In Fig. 7(a), we plot the susceptibility graphene asa function of temperature for several values of γ and µ by using Eq. (35) (lines) as well as by numerical calcu-lation of the convolution by using the sech ( βε/
2) func-tion (symbols). We can see that the approximation withthe Fadddeeva function is in a good agreement with thenumerical calculation. For comparison, we show the sus-ceptibility of undoped graphene without impurity χ byputting µ = 0 in Eq. (22) [ χ = − ( ev F ) / (6 πk B T )],which is inversely proportional to the temperature. From0Fig. 7(a), the susceptibility for non-zero γ is finite as T → µ = 0, we observe minimum values of χ at finite temperatures. For γ = 5 meV, the minimumvalue becomes smaller and shift to the higher temper-ature as we increase µ from 10 meV to 20 meV. Thepresent method reproduces the calculation by Nakamuraand Hirasawa [28], where the susceptibility of graphenewith impurity is approximated by the Sommerfeld expan-sion and also shows the minimum values in the suscepti-bility as a function of temperature. In Fig. 7(b), we plot χ/χ as a function of k B T /γ . For a given ratio µ/γ , thecurves shown in Fig. 7(a) follow the scaling law shownin Fig. 7(b). Therefore, the advantage of using the Fad-deeva function is that the susceptibility of graphene inthe presence of the impurity scattering is approximatelyscaled by the function Re[ w ( z ′ )].In Fig. 8, we numerically calculate the susceptibilityof undoped MoS as a function of T for several valuesof γ . Here, χ does not change with increasing T . Aswe increase γ , the magnitude of χ decreases with thesame rate, which means that for a given temperature,the susceptibility decreases linearly as a function of γ . IV. CONCLUSION
In this study, the analytical expressions of thermo-dynamic potentials and magnetizations of the Diracfermions are derived by using the technique of zeta func-tion regularization. There are four main results obtainedin the study: (1) the analytical formula reproduce theLangevin fitting for the magnetization of graphene fortwo limits of ~ ω c ≫ k B T and ~ ω c ≪ k B T . (2) We de-rive the formula for the magnetization of heavy Diracfermions and show that the magnetization is robust withrespect to temperature and impurity scattering. (3) Thescaling law for the susceptibility of graphene with im-purity scattering can be approximated by the real partof the Faddeeva function, and (4) the gap effect on thedHvA oscillation at T = 0 K are discussed from the prop-erty of the zeta function in the thermodynamic potential.All results by taking zeta function regularization repro-duces the previous work by taking some limits. Thus,the zeta function regularization is justiied without anyexceptions. V. ACKNOWLEDGEMENT
FRP acknowledges MEXT scholarship. MSU acknowl-edges JSPS KAKENHI Grant Number JP18J10199.RS acknowledges JSPS KAKENHI Grant NumberJP18H01810.
Appendix A: Derivation of Eqs. (8), (13), (16), and(29)1. Derivation of Eq. (8)
By substituting ǫ ξn = ξλ/ − ~ ω c q | n | + Γ ξ for n − , we getΩ − = eBh X ξ = ± " ∞ X n =0 + ∞ X n =1 ξλ − ~ ω c q n + Γ ξ − µ , (A1)where the summation operator which begins from n = 0( n = 1) indicates the sum of the LLs at the K ( K ′ ) valley.Now, let us shift the index of the summation from n = 1to n = 0 for the term − ~ ω c q n + Γ ξ as follows:Ω ( e ) − = eBh X ξ = ± " ξλ ∞ X n =1 ξλ − µ − ∞ X n =1 µ + ∆ ξ − ~ ω c ∞ X n =0 q n + Γ ξ . (A2)The first and second summations are expressed by theRiemann zeta function ζ ( p ) = P ∞ k =1 /k p , while the thirdsummation is expressed by the Hurwitz zeta function[Eq. (9)]. By using P ∞ n =1 = ζ (0) = − / − as given byEq. (8).
2. Derivation of Eq. (13)
By substituting ǫ ξn = ξλ/ ~ ω c q | n | + Γ ξ for n > ( e )+ is given byΩ ( e )+ = eBh X ξ = ± " ν ξ X n =1 + ν ξ X n =0 ξλ ~ ω c q n + Γ ξ − µ × Θ( µ − ∆ / eBh X ξ = ± " ξλ − µ ! (2 ν ξ + 1) − ∆ ξ ~ ω c ν ξ X n =0 q n + Γ ξ Θ( µ − ∆ / . (A3)Here, the summation operator which begins from n =1 ( n = 0) indicates the summations of the LLs at the K ( K ′ ) valley. By expressing the summation of n asa difference of two zeta functions [see Eq. (14)], we getEq. (13) in the main text. With the same method, byusing ǫ ξn = ξλ/ − ~ ω c q | n | + Γ ξ for n
0, we canderive Eq. (15) for Ω ( h ) − .1
3. Derivation of Eq. (16)
The thermodynamic potential Ω ′ + for T > ′ +0 ≡ Ω ′ + ( T = 0 K) with( − ∂f /∂ε ) = β sech [ β ( ε − µ ) / /
4, where f ( ε ) is theFermi distribution function [4, 25]. By using ǫ ξn = ~ ω c p | n | + Γ , Ω ′ +0 is given byΩ ′ +0 = 2 g s eBh ∞ X n = ν [ ǫ ξn − µ ]= 2 g s eBh " ~ ω c ζ (cid:16) − , Γ + ν + 1 (cid:17) + (cid:16) ν − (cid:17) µ . (A4)Here, in the second line of Eq. (A4) we use P ∞ n = ν = ζ (0) − P ν − n =1 = − ν +1 /
2. Since Γ + ν ≈ µ / ( ~ ω c ) and byconsidering µ ≫ ~ ω c , the zeta function ζ ( − / , Γ + ν +1)can be approximated by using Eq. (10). Ω ′ +0 and Ω ′ + aregiven byΩ ′ +0 ( µ ) ≈ π g s ( ~ v F ) " µ − ∆ µ − µ ( ~ ω c ) − ( ~ ω c ) µ . (A5)andΩ ′ + ( µ ) = β Z ∞−∞ dε Ω ′ +0 ( ε ) sech h β ε − µ ) i , (A6)respectively, where Ω ′ +0 ( ε ) is given by substituting µ inEq. (A5) to variable ε . Ω ′ +0 ( ε ) is an odd of ε . In thecase of µ ≪ k B T , the function sech [ β ( ε − µ ) /
2] canbe approximated as an even function. Therefore, we getΩ ′ + ( µ ) ≈ T > µ ≪ k B T .Now, by expanding the logarithmic and exponen-tial functions in the expression of thermodynamic po-tential for ~ ω c ∼ ∆ ≪ k B T , and by using ǫ ξ − n = − ~ ω c p | n | + Γ , Ω − is given byΩ − = − g s β eBh " ∞ X n =0 + ∞ X n =1 ∞ X k =1 ( − k − k e βµk exp( − βǫ ξ − n k )= g s β eBh " ∞ X n =0 + ∞ X n =1 ∞ X k =1 ( − e βµ ) k k ∞ X l =0 ( − βkǫ ξ − n ) l l != g s β eBh ∞ X l =0 ∞ X k =1 ( − e βµ ) k k − l ( β ~ ω c ) l l ! × " ∞ X n =0 ( n + Γ ) l/ − Γ l = 2 g s eBβh ∞ X l =0 Li − l ( − e βµ ) ( β ~ ω c ) l l ! " ζ (cid:16) − l , Γ (cid:17) − Γ l . (A7)In the third line of Eq. (A7), we switch the order of sum-mations with indices k and l in order to express the µ dependence of Ω − in term of the polylogarithm function P ∞ k =1 ( − e βµ ) k /k − l = Li − l ( − e βµ ), and we shift the sumoperator which begins from n = 1 to n = 0 in order to ad-just the Hurwitz zeta function. Similarly, the expressionfor Ω + is given byΩ + = 2 g s eBβh ∞ X l =0 Li − l ( − e βµ ) ( − β ~ ω c ) l l ! " ζ (cid:16) − l , Γ (cid:17) − Γ l . (A8)When we add Eqs. (A7) and (A8), the terms odd l dis-appear, while the terms for l = 0 , , , ... are doubled.By expressing l = 2 ℓ , we get Eq. (16).
4. Derivation of Eq. (29)
The Ω − for TMDs is given as follows:Ω − = − β eBh X ξ = ± " ∞ X n =0 + ∞ X n =1 ln[1 + exp {− β ( ǫ ξ − n − µ ) } ] . (A9)By using ǫ ξ − n = ξλ − ∆ / − ( ~ ω c ) | n | / ∆ ξ and separatingthe term n = 0 from n = 1 to ∞ in the sum of n = 0 to ∞ , we haveΩ − = − β eBh X ξ = ± ln[1 + e β ( µ − ξλ +∆ / ] − β eBh X ξ = ± ∞ X n =1 ln[1 + e β { µ − ξλ +∆ / ~ ω c ) n/ ∆ ξ } ] ≡ Ω ′− + Ω ′′− . (A10)Here, we define Ω ′− and Ω ′′− as the thermodynamic po-tentials for the zeroth LL at the K valley and n − ′′− , we getΩ ′′− = 2 β eBh X ξ = ± ∞ X k =1 [ − e ( µ − ξλ +∆ / ] k k × ∞ X l =0 (cid:16) βk ∆ ξ (cid:17) l ( ~ ω c ) l l ! ∞ X n =1 n l = 2 β eBh X ξ = ± ∞ X l =0 Li − l [ − e β ( µ − ξλ +∆ / ] (cid:16) β ∆ ξ (cid:17) l × ζ ( − l )= ∞ X l =0 Ω ′′ l (A11)Because Li ( z ) = − ln(1 − z ) and ζ (0) = − /
2, we getΩ ′′ = eB/ ( βh ) P ξ = ± ln[1 + e ( µ − ξλ +∆ / ] = − Ω ′− . Asa result, only the terms for l > − [Eq. (29)]. Similarly, by using ǫ ξn =∆ / ~ ω c ) | n | / ∆ ξ for n >
0, we can derive Eq. (30).2 -20-100 ζ (cid:153)(cid:154)(cid:155)(cid:156)(cid:157) p= (cid:158)(cid:159)(cid:160)¡¢£⁄ ζ ¥ƒ§¤'“«‹› fifl(cid:176)–†‡·(cid:181)¶•‚„ p= ”»…‰(cid:190)¿(cid:192)`´ˆ˜¯ ˘˙¨(cid:201)˚¸(cid:204)˝˛ˇ—(cid:209)(cid:210) (cid:211)(cid:212)(cid:213)(cid:214)(cid:215)(cid:216)(cid:217)(cid:218)(cid:219)(cid:220)(cid:221)(cid:222) (cid:223)(cid:224)Æ (cid:226)ª(cid:228) (cid:229)(cid:230)(cid:231) ŁØŒº(cid:236)(cid:237)(cid:238)(cid:239)(cid:240)æ(cid:242)(cid:243)(cid:244) p= ı(cid:246)(cid:247)łøœß (cid:252) Z (cid:253)(cid:254)(cid:255) q Z -505 (cid:11)(cid:12)(cid:13)(cid:14)(cid:15)(cid:16)(cid:17)(cid:18) (cid:24)(cid:25)(cid:26)(cid:27)(cid:28)(cid:29)(cid:30)(cid:31) FIG. B.1. (Color online) Plot of ζ ( p, q ) with q = 1 − p = − /
2, (b) p = 1 /
2. In the insets of (a) and (b)we show that the values of ζ ( p, q ) for p = − / p = 1 / q > .
66 and q > .
3, respectively. Plot of ζ ( p, q, N +1) with N = 100 for (c) p = − /
2, (d) p = 1 /
2. Thecomparison between the functions S ( p, q, N ) and Z ( p, q, N )for given p , q , and N are shown in (e) and (d). ! " -100 -6-4 +, ./1 ζ ( p , q + N + ) SZ SZp=1/2(a)p=−1/2 p=1/2p=−1/2(c) (b)(d)
FIG. B.2. (Color online) Plot of ζ ( p, q + N + 1) with N =10 ⌊ q ⌋ for (a) p = − /
2, (b) p = 1 /
2. The comparison betweenthe functions S ( p, q, N ) and Z ( p, q, N ) for given p , q , and N are shown in (e) and (f). Appendix B: Numerical calculations of the zetafunction
In Fig. (B.1) (a) and (b), we plot the zeta func-tion ζ ( p, q ) for p = 1 / p = − /
2, respectively.The value of ζ (1 / , q ) is negative and deceases mono-tonically for q > .
66 as shown in the inset in (a).The value of ζ ( − / , q ) diverges at q = 0 and changesigns at q ≈ .
3. In (c) and (d), we substitute q to q + N + 1 and take N = 100 for explaining the changeof ζ ( p, q ) − ζ ( p, q + N + 1). In (e) and (f) we com- pare the functions S ( p, q, N ) ≡ P Nn =0 ( n + q ) − p and Z ( p, q, N ) ≡ ζ ( p, q ) − ζ ( p, q, N + 1), similar to the left-hand and right-hand sides of Eq. (14), respectively. Itis observed that the two functions exactly identical for q = 0 to 10. It is noted that both the functions S and Z are continuous and do not explain the oscilla-tory behaviour of the thermodynamic potential in thedHvA effect. By changing the constant N to 10 ⌊ q ⌋ foran example, the function ζ ( p, q, N + 1) shows step-likebehaviour as shown in Fig. B.2(a) and (b) because of thenature of the function ⌊ q ⌋ . In (c) and (d), we comparethe functions S ( p, q, N ) and Z ( p, q, N ) for p = 1 / p = − /
2, respectively. As in the previous case, the twofunctions match each other. Therefore, the analyticalexpressions of the thermodynamic potentials for doped aDirac fermion is numerically verified.
Appendix C: Approximation of sech ( βǫ/ to aGaussian function
56 89 :;< ε [meV]0 =>?@ABCDEFGH fI ε JKLM ε ) TNOP QRU VWXY Z[\] ^_‘a b
FIG. C.1. (Color online) Comparison between the functions f ( ε ) and g ( ε ) (thin solid lines) for T = 25 , , ,
200 and300 K.
For given secant-hyperbolic and the Gaussian distribu-tions, F ( ε ) ≡ sech( ε/W ) and G ( ε ) ≡ exp[ − ε / (2 σ )], re-spectively, the half-width of the distributions are given byHW F = ln(2+ √ W and HW G = p σ . By solvingHW F = HW G and choosing W = 2 /β , the Gaussian ap-proximation for the function f ( ε ) ≡ sech ( βε/
2) is givenby g ( ε ) ≡ exp[ − ( Cβε ) ], where C = √ ln2 / [ √ √ ≈ .
447 as defined in the main text. In Fig. C.1, wecompare f ( ε ) and g ( ε ) (thin solid lines) for several val-ues of temperature. The distribution g ( ε ) has a smallertail compared with f ( ε ), which is the origin of discrepan-cies between the numerical calculation and the Faddeevaapproximation in the calculation of χ ( µ, γ ). [1] M. Sepioni, R. R. Nair, S. Rablen, J. Narayanan, F.Tuna, R. Winpenny, A. K. Geim, and I. V. Grigorieva, Phys. Rev. Lett., , 207205 (2010). [2] M. O. Goerbig, Rev. Mod. Phys. , 1193 (2011).[3] Y. Fuseya, M. Ogato, H. Fukuyama, J. Phys. Soc. Jpn. , 012001 (2015).[4] J. W. McClure, Phys. Rev. , 666 (1956).[5] R. Saito and H. Kamimura, Phys. Rev. B , 7218(1986).[6] A. Raoux, M. Morigi, J.-N. Fuchs, F. Pi´echon, and G.Montambaux, Phys. Rev. Lett. , 026402 (2014).[7] M. Koshino and T. Ando, Phys. Rev. B , 195431(2010).[8] M. Koshino and T. Ando, Solid State Commun. ,1054 (2011).[9] H. Fukuyama, Y. Fuseya, M. Ogato, A. Kobayashi, andY. Suzumura, Physica B , 1943 (2012).[10] L. Landau. Z. Phys. , 629 (1930).[11] T. Cai, S. A. Yang, X. Li, F. Zhang, J. Shi, W. Yao, andQ. Niu, Phys. Rev. B , 115140 (2013).[12] M. Koshino and I. F. Hizbullah, Phys. Rev. B , 045201(2016).[13] X.-Y. Yan and C. S. Ting, Phys. Rev. B , 104403(2017).[14] Z. Li, L. Chen, S. Meng, L. Guo, J. Huang, Y. Liu, W.Wang, and X. Chen, Phys. Rev. B , 094429 (2015).[15] D. Cangemi and G. Dunne, Ann. Phys. (N.Y.) , 582(1996).[16] A. Ghosal, P. Goswami, and S. Chakravarty. Phys. Rev.B , 115123 (2007).[17] S. Slizovskiy and J. J. Betouras. Phys. Rev. B , 125440(2012).[18] G.-B. Liu, W.-Y. Shan, Y. Yao, W. Yao, and D. Xiao,Phys. Rev. B. , 085433 (2013).[19] W.-K. Tse and A. H. MacDonald, Phys. Rev. B ,205327 (2011).[20] J. L. Lado and J. Fern´andez-Rossier, 2D Mater. , 035023(2016).[21] C. J. Tabert and E. J. Nicol, Phys. Rev. Lett. ,197402 (2013).[22] F. Qu, A. C. Dias, J. Fu, L. Villegas-Lelovsky, and D. L.Azevedo, Sci. Rep. , 41044 (2017).[23] X. Li, F. Zhang and C. Niu, Phys. Rev. Lett. , 066803(2013).[24] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Phys.Rev. Lett. , 196802 (2012).[25] S. G. Sharapov, V. P. Gusynin, and H. Beck, Phys. Rev.B. , 075104 (2004). [26] M. Koshino and T. Ando, Phys. Rev. B. , 235333(2007).[27] M. Nakamura, Phys. Rev. B. , 113301 (2007).[28] M. Nakamura and L. Hirasawa, Phys. Rev. B. , 045429(2008).[29] C. J. Tabert, J. P. Carbotte, and E. J. Nicol, Phys. Rev.B. , 035423 (2015).[30] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C.W. Clark (eds.), NIST Handbook of Mathematical Func-tions, Cambridge Univ. Press (2010).[31] M. Ezawa, Eur. Phys. J. B , 363 (2012).[32] I. A. Luk’yanchuk and Y. Kopelevich, Phys. Rev. Lett. , 166402 (2004).[33] T. Champel and V. P. Mineev, Philos. Mag. B , 55(2001).[34] I. A. Luk’yanchuk, Low Temp. Phys. , 45 (2011).[35] A. R. Wright and R. H. McKenzie, Phys. Rev. B ,085411 (2013).[36] K. Kishigi and Y. Hasegawa, Phys. Rev. B , 085427(2014).[37] I. M. Lifshits and A. M. Kosevich, Zh. Eksp. Teor. Fiz. , 730 (1955) [Sov. Phys. JETP , 636 (1956)].[38] L. Onsager, Phil. Mag. , 1006 (1952).[39] 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.
770 (2007).[40] F. Escudero, J. S. Ardenghi, and P. Jasen, J. Phys. Con-dens. Matter , 285804 (2019).[41] F. Escudero, J. S. Ardenghi, and P. Jasen, Eur. Phys. J.B , 93 (2020).[42] C. H. Yang, F. M. Peeters, and W. Xu, Phys. Rev. B , 146401(2015).[46] S. Becker and M. Zworski, Commun. Math. Phys.
941 (2019).[47] Y. Kubota, K. Watanabe, O. Tsuda, and T. Taniguchi,Science , 932 (2007).[48] K. K. Kim, A. Hsu, X. Jia, S. M. Kim, Y. Shi, M. Hof-mann, D. Nezich, J. F. Rodriguez-Nieva, M. Dresselhaus,T. Palacios, and J. Kong, Nano Lett. , 161 (2011).[49] G. Cassabois, P. Valvin, and B. Gil, Nat. Photonics10