Landau levels, response functions and magnetic oscillations from a generalized Onsager relation

A generalized semiclassical quantization condition for cyclotron orbits was recently proposed by Gao and Niu \cite{Gao}, that goes beyond the Onsager relation \cite{Onsager}. In addition to the integrated density of states, it formally involves magnetic response functions of all orders in the magnetic field. In particular, up to second order, it requires the knowledge of the spontaneous magnetization and the magnetic susceptibility, as was early anticipated by Roth \cite{Roth}. We study three applications of this relation focusing on two-dimensional electrons. First, we obtain magnetic response functions from Landau levels. Second we obtain Landau levels from response functions. Third we study magnetic oscillations in metals and propose a proper way to analyze Landau plots (i.e. the oscillation index $n$ as a function of the inverse magnetic field $1/B$) in order to extract quantities such as a zero-field phase-shift. Whereas the frequency of $1/B$-oscillations depends on the zero-field energy spectrum, the zero-field phase-shift depends on the geometry of the cell-periodic Bloch states via two contributions: the Berry phase and the average orbital magnetic moment on the Fermi surface. We also quantify deviations from linearity in Landau plots (i.e. aperiodic magnetic oscillations), as recently measured in surface states of three-dimensional topological insulators and emphasized by Wright and McKenzie \cite{Wright}.


Introduction
The quantization of closed cyclotron orbits for Bloch electrons in the presence of a magnetic field leads to the formation of Landau levels (LLs) [5]. When the Fermi level crosses the LLs, physical quantities such as the resistance, the magnetization or the density of states feature quantum magnetic oscillations [6]. The analysis of the latter is considered as a powerful way of extracting informations on the band structure of metals [7]. Band structure here refers to both the zero-field band energy spectrum n (k) and to the geometry of the cell-periodic Bloch states |u n (k) . The standard way to analyze magnetic oscillations is to use Onsager's quantization condition [2] to obtain Landau levels for Bloch electrons and the Lifshitz-Kosevich formula [8] that describes the temperature and disorder dependence of the amplitude of oscillations.
As early as 1966, Roth generalized Onsager's condition by including inter-band effects in the semiclassical quantization condition of closed cyclotron orbits up to second order in the magnetic field [3] (see also Appendix B for other important contributions). In a recent insightful paper, Gao and Niu [1] proposed a further extension by systematically including higher-order corrections in the magnetic field in a compact and thermodynamic manner. Their equation generalizes Onsager's relation [2] (n + 1 which relates N 0 ( ), the zero-field integrated density of states (IDoS) at the Fermi energy , to the degeneracy eB/h of a LL, where B > 0 is the modulus of the magnetic field, and the Landau index n, which is an integer [here we assumed a sample area A = 1]. The relation obtained by Gao and Niu reads where N ( , B) is the smoothed (i.e. without the magnetic oscillations [9]) IDoS in the presence of a magnetic field B, M 0 ( ) is the spontaneous magnetization, χ 0 ( ) is the magnetic susceptibility, R p are higher order magnetic response functions and the prime denotes the derivative with respect to the energy M 0 = ∂ M 0 . These quantities N 0 , M 0 , χ 0 , R p are taken at zero temperature and in the limit of zero magnetic field, and depend on the Fermi energy . The conceptual important difference between (1) and (2) is that whereas the first relies only on the zero field energy spectrum (or more precisely on the zero-field IDoS), the supplementary terms in the second explicitly involve information on the zero field cell-periodic Bloch states and on interband coupling. This extra information is carried by quantities such as the magnetization and the susceptibility. Indeed, the magnetization involves the orbital magnetic moment and the Berry curvature [10], and the susceptibility involves not only the curvature of the energy spectrum but also the Berry curvature and the quantum metric [11].
Interestingly equation (2) may be rewritten in a form similar to (1) [n + γ( , B)] eB h = N 0 ( ) , but involving an energy and magnetic-field dependent phase-shift instead of a mere constant 1 2 . Roth wrote a formally similar equation (see (42) in [3]), however she did not relate the first (resp. second) order correction to the derivative of the magnetization (resp. magnetic susceptibility) and did not obtain the higher-order corrections in the field. In the following, we will call eq. (2) the Roth-Gao-Niu relation.
There are at least three important consequences of this generalized Onsager quantization condition: 1) When LLs are known analytically, their expression can be inverted to obtain response functions analytically (actually their derivatives ∂M 0 ∂ and ∂χ 0 ∂ ). 2) Response functions can be used to obtain LLs in cases where the latter are hard to obtain exactly and for which response functions are known by other means (i.e. linear response theory).
3) The phase of magnetic oscillations (related to γ( , B)), such as Shubnikov-de Haas oscillations in the longitudinal resistance or de Haas-van Alphen oscillations in the magnetization, can be derived from the Roth-Gao-Niu relation. This helps to analyze Landau plots (i.e. index n of oscillations as a function of the inverse magnetic field 1/B). In particular, we can obtain a simple formula for the zero-field phase-shift γ( , 0) in terms of the Maslov index, the Berry phase and the orbital magnetic moment (averaged over the Fermi surface). Also we can see the general structure of the γ( , B) and study deviations from linear Landau plots, i.e. aperiodic magnetic oscillations, as recently measured in topological insulator surface states and discussed by Wright and McKenzie [4].
The aim of the present paper is to elaborate on these three consequences. The article is organized as follows. We first review the Roth-Gao-Niu quantization condition and its validity (section 2), then present the three type of consequences: from LLs to response functions (section 3) on the examples of the graphene monolayer and bilayer, a generic 2D semiconductor and the Rashba model; from response functions to LLs (section 4) on the examples of the Hofstadter, the semi-Dirac and the bilayer models; and an analysis of magnetic oscillations in section 5. Eventually, we give a conclusion in section 6. Some material is also presented in appendices.

Roth-Gao-Niu quantization condition
We start by shortly reviewing the Roth-Gao-Niu quantization condition and then discuss its validity. It is argued [1] that for a half-filled LL (i.e. when the Fermi energy is exactly at the energy n of the nth LL) the number of states below the Fermi level at zero temperature is where Ω( , B, T ) is the non-oscillatory part (the smooth part) of the grand potential and the prime denotes the partial derivative with respect to the chemical potential Ω = ∂ Ω. Here, we assumed that the sample area A = 1 so that the degeneracy of a LL is N φ = eBA h = B with = 1 and e = 2π such that the flux quantum φ 0 = h/e = 1. The temperature is assumed to be larger than the typical splitting between LLs: T ω c ∝ B with T → 0 and B → 0. This requirement ensures that one is working with the low field expansion of the smooth (non-oscillatory) part of the grand potential (see the discussion in Appendix A). Note that the above relation is far from obvious as the left hand side denotes the quantum mechanical IDoS at zero temperature, while the right hand side is for the smoothed IDoS, typically obtained at temperature higher than the separation between LLs.
It is assumed that the smooth grand potential Ω may be written as a power series in B (see below for a discussion of the validity of this assumption). Its expansion is written as where M 0 is the zero-field (spontaneous) magnetization, χ 0 is the zero-field magnetic susceptibility and R p ( ) = − ∂ p Ω ∂B p ( , B = 0) with p ≥ 3 are higher order magnetic response functions (all at zero temperature). From the expansion, it follows that where N 0 ( ) = −Ω 0 ( ) is the zero-field IDoS. This is the Roth-Gao-Niu relation [1]. When keeping only N 0 ( ) in the r.h.s., it reduces to the Onsager quantization condition [2]. When keeping first order correction to the Onsager relation, N 0 ( ) = (n + 1 2 − M 0 ( ))B, it shows the appearance of Berry phase type of corrections -hidden in M 0 -and recovers various results scattered in the literature. When keeping second order corrections, N 0 ( ) = (n + 1 2 − M 0 ( ))B − χ 0 ( ) B 2 2 , it is formally equivalent to the Roth quantization condition [3], albeit written in a much more compact and transparent form. This will be discussed in detail in section V.
We now discuss the validity of the Roth-Gao-Niu relation. There are several issues: (i) A first issue concerns its application to a system with degeneracies such as several valleys or spin projections. In fact, it is only meaningful for each species separately [1]. This implies that on the right hand side of Eq.(7), the effective zero field quantities N 0 ( ), M 0 ( ) and χ 0 ( ) are species dependent and thus do not correspond to directly measurable equilibrium thermodynamic responses. In particular, for time reversal invariant systems, there is no finite thermodynamic spontaneous magnetization but the species dependent effective spontaneous magnetization contribution M 0 ( ) might nevertheless appear finite.
(ii) Equation (7) is valid for electron-like but not for hole-like cyclotron orbits. In the latter case, it becomes (8) where N tot is the total number of electrons when the band is full and n ∈ N. Actually is the DoS. Generally speaking N tot is a constant that is either determined by the occupation of a full band in the case of a model defined on a lattice or should be determined from self-consistency in the case of an effective low-energy model. See also Appendix D where the Onsager quantization conditions for hole orbits is discussed.
(iii) Semiclassical quantization condition for closed cyclotron orbits (either Roth-Gao-Niu or Onsager) predicts perfectly degenerate Landau levels and does not account for lattice broadening of Landau levels into bands [12]. Indeed, it does not include tunneling between cyclotron orbits belonging either to different valleys or to another Brillouin zone, i.e. magnetic breakdown [13]. All these phenomena are beyond the present description. Neglecting magnetic breakdown is consistent with the assumption that the grand potential admits a series expansion in (positive) integer powers of B and does not contain terms such as e −#/B .
(iv) The Roth-Gao-Niu relation fails to capture singular behaviors in response functions, that may appear at specific energies corresponding to band contacts or edges. More precisely, it misses step functions or delta functions in N 0 , M 0 , χ 0 , etc. As a first example, the Roth-Gao-Niu relation does not account for the McClure diamagnetic delta-peak in the susceptibility at the Dirac point of graphene [14]. Similarly for massive Dirac fermions (gapped graphene or boron nitride), the step functions at the gap edges of the diamagnetic susceptibility plateau [15] are not accounted for by the present formalism. In the first example, at the energy of the band contact, the grand potential actually behaves as B 3/2 and therefore does not admit an expansion in integer powers of B, as is assumed in the Roth-Gao-Niu relation.
(v) The right-hand side of eq. (2) is an asymptotic series in powers of B. Here we stress that the small parameter of this expansion is B at fixed (n + 1 2 )B. The fact that the lefthand side of eq. (2) is fixed comes from the quantization occurring at constant energy and from Onsager's relation (n + 1 2 )B ≈ N 0 ( ) = constant. Because (n + 1 2 )B is fixed, the small parameter B can also be thought as 1 n+1/2 , where one recognizes the usual semi-classical criterion of large n.

From Landau levels to magnetic response functions
When available, one may use the knowledge of the exact LLs to obtain magnetic response functions such as the magnetization and the susceptibility. In some cases, it may be easier to proceed that way rather than to try to compute directly these response function using linear response theory. In the following, we treat four examples in detail: (1) gapped graphene monolayer, (2) gapped graphene bilayer, (3) gapped Dirac electrons for a semiconductor and (4) the Rashba model. The three first examples are gapped and have a valence and a conduction band. The fourth example has a Fermi surface that can be electron or hole-like. It is therefore essential to be able to treat both electrons and holes so that the Roth-Gao-Niu relation should be adapted to account both for electrons and holes contributions. It is thus convenient to regroup equations (7) and (8) into a single relation where s = sign( ) is the sign of the energy and sN 0 ( )+N tot Θ(−s) is simply the IDoS of either electrons N el ( ) = N 0 ( ) or holes N h ( ) = N tot − N 0 ( ) depending on s. It gives the Landau index n as a Laurent series in the magnetic field B starting with a 1/B term: n = ∞ p=−1 c p B p . The main result of this section is that the Roth-Gao-Niu relation indeed allows one to recover many results concerning the energy-derivative of magnetic response functions. However, it fails to recover singular behaviors such as step functions, Dirac delta functions, etc.

Gapped graphene monolayer with Zeeman effect
We consider the low energy description of a gapped graphene monolayer (e.g. boron nitride) in the presence of a Zeeman effect. The Hamiltonian (with ξ = ±1 the valley index) is [22] H ξ = (vξΠ x τ x + vΠ y τ y + ∆τ z )σ 0 + ∆ Z σ z τ 0 , where Π = p + 2πA is the gauge-invariant momentum, τ are sublattice pseudo-spin Pauli matrices, σ are real spin Pauli matrices, ∆ Z = g 2 µ B B is the Zeeman energy (µ B = e 2me = π me is the Bohr magneton), the magnetic field B > 0 is assumed to be along the z direction perpendicular to the conduction plane and ∆ is a staggered on-site energy. In the following we take units such that v = 1. The LLs are where λ = ±1 is the band index, σ = ±1 is the spin index,n = n + 1−λξ 2 (n ∈ N such that n ∈ N if λξ = +,n ∈ N * if λξ = −) and ω v ≡ √ 4πB. It is important to treat correctly the n = 0 LL in order to account for the parity anomaly. Therefore: Here the comparison with equation (9) gives where Θ ≡ Θ[| | − ∆] is defined in order to avoid cluttered notation and here λ = sign( ). These quantities are plotted in Figure 1, together with the zero-field energy spectrum and the LLs. The DoS agrees with that found by Koshino and Ando [15]. TRS implies that when summing over spin and valley indices, the spontaneous magnetization should vanish and therefore M 0,ξ, 2 comes from the peculiar behavior of the winding number W = λξ of gapped graphene, which leads to a phase shift γ 0 (see section 5 below) which is energy and magnetic field independent, γ 0 = 1 2 − M 0,ξ = 1 2 − W 2 [16]. The quantization γ 0 = 0 mod. 1 is actually only exact for the simplified two-band model of gapped graphene that we consider here, see [4,17,18] for a discussion. After integrating χ 0 , one recognizes two contributions, the Pauli spin paramagnetism χ spin ( ) = ( g 2 µ B ) 2 ρ 0,ξ,σ ( ) and a flat orbital susceptibility χ orb ( ) = const, which cannot be obtained from this method. Actually, Koshino and Ando find that the orbital susceptibility is piecewise constant χ orb ( ) = − π 3∆ Θ(∆ − | |) [15] rather than constant. This means that the derivative is singular at the gap edge: χ orb ( ) = π 3∆ λδ(| |−∆). The singularity is not captured by the present approach based on the Roth-Gao-Niu quantization condition. In the gap closing limit ∆ → 0, the Koshino-Ando orbital susceptibility recovers the singular diamagnetic delta peak χ orb ( ) = − 2π 3 δ( ), first obtained by McClure [14]. Whenever the energy spectrum in the presence of a magnetic field is the sum of an orbital part and of a spin part (i.e. without cross-terms) as in eq. (11), the magnetic response functions are just the sum of an orbital and of a spin response functions. Generally, for systems with TRS, all odd derivatives (R 1 = M 0 , R 3 , etc.) are proportional to the valley index ξ, such that upon summing over both valleys, they vanish as expected.

Gapped graphene bilayer
The low energy gapped graphene bilayer Hamiltonian is [19] where ξ = ± is the valley index, m is an effective mass and 2∆ is a gap. The exact LLs are [19] (n ≥ 0) where ω 0 = 2πB m ,n = n + 1 + λξ with λ = ± the band index (the n = 0 and n = 1 LLs given in [1], which are supposed to be those for λ = + and ξ = +, are not correct). Inverting this relation we obtain: Substituting this relation in the l.h.s. of eq.(9) and expanding as a series in powers of B, we obtain Order by order comparison with the r.h.s. implies that where as before Θ ≡ Θ[| | − ∆] and λ = sign( ). The above results are plotted in Figure 2. They are compatible with those obtained in [1] except for the sign of M 0,ξ [20]. The derivative of the spontaneous magnetization vanishes when summed over both valleys, as it should for a system with TRS. Also M 0,ξ = −ξ corresponds to a winding number W = 2λξ as shown in [16] through the relation γ 0 = 1 2 − W 2 for the phase shift (see section 5 below). The result for the derivative of the susceptibility agrees with the response function found by Safran [21] where t ⊥ is a high-energy cutoff related to interlayer hoping (see also [23,24]).
We conclude this section by discussing a related but different case, namely that of a quadratic band crossing point that occurs at a time-reversal invariant momentum, such as the M point of the checkerboard (Mielke) lattice [25]. The low-energy Hamiltonian reads The main differences with the gapless graphene bilayer Hamiltonian, eq. (14) with ∆ = 0, are in the involved Pauli matrices, which reflect time-reversal invariance and the fact that there is a single valley. In this situation, the Landau levels are n,λ (B) = λω 0 n(n + 1), where n = 0, 1, 2... and λ = ±. Note that the change in Pauli matrices does not affect the energy spectrum but it affects the labelling (i.e. the quantum numbers n and λ), which matters to deduce the response functions. Inverting this equation, we find which, in contrast to eq. (16) does not contain the valley term λξ. Then we deduce : We then correctly find that the magnetization (actually its derivative) is zero in this case, while it was finite and of opposite signs in the previous case with two valleys. In both cases the total magnetization is zero due to TRS.

Two-dimensional semiconductors with gapped Dirac electrons
The low energy Hamiltonian of gapped Dirac electrons in two-dimensional semiconductors such as transition metal dichalcogenides (see e.g. [17]) is where ξ = ±1 is the valley index and (τ x , τ y , τ z ) are pseudo-spin Pauli matrices. In the following, we take units such that v = 1. The exact LLs are given by (for n ≥ 0): where and with λ = ±1 the band index. Inverting this relation, we obtain where we have defined 1 The assumption that 1 mz > 1 m insures that the Fermi surface is made of a single piece. Using relation eq. (9), and expanding to order B 2 we obtain successively: and where Θ ≡ Θ[| | − ∆] and λ = sign( ). These results are plotted in Figure 3. Despite a complicated expression, they correspond to simple functions. The derivative of the magnetization vanishes when summed over both valleys as expected for a system with TRS. This is actually the case of all odd response functions.
Equation (28) shows that, for this 2D semiconductor model, the valley magnetization contribution is energy dependent and not quantized, in contrast to gapped Dirac electrons (see section 3.1). As a consequence, it would contribute an energy-dependent phase shift in magnetic oscillations (see section 5 and [4,18]).

Rashba model
As a last example, we consider the case of a two-dimensional electron gas with Rashba spinorbit coupling and a Zeeman effect, known as the Rashba model [26]. In contrast to the previous examples, the spectrum of this model is gapless. The Hamiltonian reads [26] where m is the mass, v the velocity (i.e. the Rashba parameter quantifying the spin-orbit coupling), ∆ Z the Zeeman energy and the magnetic field B > 0 is assumed to be along the z direction perpendicular to the conduction plane. The Pauli matrices (σ x , σ y , σ z ) describe the true spin and there is a single valley (e.g. at the A point of the hexagonal Brillouin zone in the case of BiTeI, see [27]). This is also an effective description of the surface states of threedimensional topological insulators [4,28,29]. However, there is a subtle distinction between the case of the Rashba model and that of the topological insulator: in the first case, the Fermi surface is made of two disconnected pieces because the parabolic term in the Hamiltonian is not a small correction to the linear term, in contrast to the second case in which there is a single Fermi surface. The correct treatment of the Rashba model therefore requires great care. In the following we take units such that v = 1 in addition to = 1 and e = 2π. The dispersion relation is made of two bands λ (k) = k 2 2m +λk where λ = ± is a band index (note that here λ is not the same as the sign of the energy) and k = |k|. The Fermi surface is made of two pieces: an inner circle with radius k i and an outer circle with radius k o (see Figure   4). At positive energy, the inner circle belongs to the λ = + band (k i = −m + m 1 + 2 m ) and the outer circle to the λ = − band (k o = m + m 1 + 2 m ). Both circles are electronlike. At negative energy (− m 2 < < 0), the Fermi surface is again made of two circles but both the inner and the outer circles belong to the λ = − band (k i = m − m 1 − 2| | m and k o = m + m 1 − 2| | m ) and the λ = + band no longer intersects the Fermi energy. The inner circle is hole-like and the outer circle is electron-like. For future need, we introduce a branch index b = o/i = ± to designate the two Fermi circles: o for outer (+) and i for inner (−). Therefore which is valid for b = ± and for positive or negative . In the case of the Rashba model, it is essential to consider the two Fermi circles. However, when describing surface states of topological insulators, the Rashba model is only valid for small wave-vectors and therefore only the inner circle should be considered (in Figure 4(a) or (b), it means that only the full line should be considered and not the dashed one). The latter model was already discussed by Gao and Niu [1]. It leads to inconsistencies such as a non-vanishing magnetization despite TRS.
The LLs are [26,30] n,λ = ω 0n + λ ω 2 vn + where and we assume thatg < 1. In particular, n = n + 1−λ 2 with n = 0, 1, 2, ..., which means thatn ∈ N when λ = + andn ∈ N * when λ = − so that the parity anomaly at n = 0 is correctly taken into account. Indeed atn = 0 there is a single LL (λ = + with our choice), whereas for anyn > 1, there are always two LLs (one with λ = + and one with λ = −). It will be important to keep in mind that some of these LLs are electron-like (inner and outer circles at positive energy and outer circle at negative energy) and some are hole-like (inner circle at negative energy), see Figure 4(c) and 5(a).
Inverting this equation with care, we find . Upon expansion in powers of B, this gives Equation (9) was obtained for electron-like LLs and is therefore valid only for the outer circle (for all energies) and for the inner circle at positive energy. However, we need to modify it for the hole-like LLs corresponding to the inner circle at negative energy. In that case, see equation (8), it becomes where N tot is a constant. By comparing eq. (34) with eq. (9) and (35), we obtain the IDoS, the magnetization and the susceptibility. We start with the IDoS. For > 0, the IDoS is featuring a 1D-like van Hove singularity when → − m 2 corresponding to the minimum of the mexican hat dispersion relation. In (35), we took N tot = 0 in order that the total IDoS be continuous at = 0. The IDoS agrees with The change of sign in front of k 2 i is due to the fact that the inner Fermi circle is hole-like. In the end, restoring the units, the IDoS is which tends to m π 2 Θ( ) in the v → 0 limit as expected (see Figure 5). The derivative of the magnetization is M 0,b = − b 2 so that M 0 = b M 0,b = 0 and the total magnetization M 0 ( ) = const as a consequence of the compensation between the inner and outer Fermi circles. TRS actually implies that M 0 ( ) = 0. Compare with [30].
We obtain the derivative of the susceptibility as: For positive energy and the inner circle b = i = −, this result agrees with [1] (see in particular the supporting information of this article). The total (derivative of the) susceptibility is therefore and is displayed in Figure 5. Except for a delta-like singularity at = 0, which is the McClure result mentioned in the section on graphene [14], this compares well with the susceptibility found by Suzuura and Ando [30] χ 0 ( ) = π m (1 −g) 2 the derivative of which is: However it slightly disagrees with [27]. In the vicinity of = 0, these authors find a diamagnetic peak ∝ −

From magnetic response functions to Landau levels
One may also use the Roth-Gao-Niu quantization condition in order to obtain an approximate analytical expression of LLs. As an input, one needs zero-field quantities such as the IDoS N 0 ( ), the spontaneous magnetization M 0 ( ) and the magnetic susceptibility χ 0 ( ).
The general structure of LLs obtained by including successive powers of B in the Roth-Gao-Niu relation (7) can be shown to be where x ≡ (n + 1 2 )B and the functions g p (x) depend on the magnetic response functions. At zeroth order in B, g 0 (x) = N −1 0 (x) = (0) n , which is the Onsager result. At first order, and at second order: At pth order, the deviation with respect to the exact LLs is expected to be of order B p+1 (at fixed x).

Square lattice tight-binding model
We consider the square lattice tight-binding model in a perpendicular magnetic field (the famous Hofstadter model [32]) and choose units such that a = 1 (and not A = 1 here), t = 1 and φ 0 = 1 ( = 1 and e = 2π) so that the flux per plaquette is f = Ba 2 /φ 0 = B. The idea is to use the Roth-Gao-Niu relation to reproduce features of the Hofstadter spectrum. This example is particularly interesting since the quantities entering the Roth-Gao-Niu relation are known analytically. In particular, we recover a semi-classical expression of the LLs up to third order in the magnetic flux, originally obtained by Rammal and Bellissard [33].
As the magnetization vanishes, the Roth-Gao-Niu relation reads: Time-reversal symmetry implies that the expansion in powers of f of the grand potential only contains even powers. The above expression is therefore valid up to order f 4 . The zero-field DoS ρ 0 ( ) = N 0 ( ) is [34] ρ 0 ( ) = 1 where K(x) is the complete elliptic integral of the first kind [35,36]. The susceptibility χ 0 has been also calculated and can be written in the form [24,37]: E( ) = 0 νρ 0 (ν)dν is the total energy and Q ν (z) is a Legendre function of the second kind [35]. The derivative χ 0 ( ) is Inserting in Eq.(44) the expressions of the IDoS N 0 ( ) obtained from integrating Eq.(45) and of the derivative χ 0 ( ), we obtain the position n (f ) of the LLs in the inverted form: One recovers the Onsager quantization rule in the limit of a constant susceptibility (χ 0 → 0). By solving this equation, one obtains the field dependence of the LLs n (f ). They are plotted in Fig. 6 and are compared with the exact levels of the Hofstadter spectrum. It provides a very good description of the global evolution of the Hofstadter bands. More quantitatively, As shown in Fig. 6, this excellent fit works as long as a semi-classical quantization relation is applicable, that is as long as the number of states in a broadened Landau level equals the Landau degeneracy f A. When f max (n) = 1 2(n+1) , the broadened level n overlaps with its symmetric in energy (coming from the top of the band) and the broadened level begins to empty instead of following the Landau degeneracy.
From the implicit Eq. (48), one can easily obtain a series expansion of the field dependence n (f ). From Eqs. (45,46), we obtain the following expansions, defining µ = + 4 Remarkably, we recover a result obtained by Rammal and Bellissard (equation (3.28) in [33]) after a lengthy calculation from a totally different method. As we are using a Roth-Gao-Niu relation valid up to order f 3 included (but not at order f 4 as we do not know the response functions R p with p ≥ 4), there is no point in getting the semi-classical LLs beyond that order. Note that using the Onsager relation instead already leads to problems at order f 2 for the LLs. The above LLs can be rewritten in terms of the combinaison 2n + 1 as: GN The Onsager quantization gives the following LLs On (2n + 1) 3 + 0 4 + . . . , (53) which are a function of (2n + 1)f only. In order to compare these two results and their order in the semi-classical expansion, we take (2n + 1)f as a constant and count the powers of f (see the corresponding discussion in point (v) of 2). The comparison of (52) and (53) gives (where n = 0, 1, 2, · · · is the Landau index) and start to deviate beyond that.
The approximate quantization condition (44) works remarkably well when compared to the numerics up to f max (n), see Fig. 6 and 7. We can therefore assume that this equation is qualitatively correct beyond its derived range of validity and up to fluxes of order f ∼ 1/2. In terms of magnetic oscillations (we anticipate on section 5), it gives the following Landau plot (oscillation index n as a function of f −1 ) and clearly shows a non-linear behavior ∝ 1/f −1 on top of the familiar linear slope f −1 + constant, see Fig. 8.

Two-dimensional semi-Dirac Hamiltonian
We now consider the 2D semi-Dirac Hamiltonian describing the merging point of two Dirac cones with opposite chirality [38,39]. We use the Roth-Gao-Niu relation to obtain LLs from response functions. The semi-Dirac Hamiltonian [38] is We use units such that = 1, v y = 1 and m = 1 (the energy unit is therefore mv 2 y , and that of length is /(mv y )). The zero field spectrum is ± (k) = ± k 4 x /4 + k 2 y , the density of states is ρ 0 ( ) = C √ (when ≥ 0) with C = Γ(1/4) 2 4π 5/2 ≈ 0.187857. The integrated DoS is This is the number of states below energy > 0 in the conduction band. Onsager's quantization condition gives the following semiclassical LLs [38]: We now improve on these LLs by using the Roth-Gao-Niu relation. We first truncate eq. (7) at second order in the field and wish to solve for the LLs n using our knowledge of the magnetic response functions. As there is a single valley and time-reversal symmetry, we know that M 0 = 0. In addition, we computed the susceptibility from linear response theory [40] and obtained (for > 0) [41] so that Therefore we have to solve the following equation to obtain the LLs (at next order in the semiclassical/small field expansion): This is a quadratic equation in 3/2 . We find that n = SC n g(n) with g(n) = This gives g(0) ≈ 0.918 (instead of 0.808 [38]), g(1) ≈ 0.992 (instead of 0.994 [38]) and g(n ≥ 2) ≥ 0.997 (instead of ∼ 1 [38]). At fixed (n + 1 2 )B, SC n is a constant and g(n) ∝ B 0 + B 2 + O(B 4 ). The above equation is equivalent to writing with γ n = n + 1

Gapped graphene bilayer
From the knowledge of the magnetic response functions of the gapped graphene bilayer, we now seek to obtain approximate LLs. This section can be seen as the reversed path (LLs from response functions) compared to what we did in section 3.2, in which we obtained response functions from exact LLs. This case was already studied by Gao and Niu [1], who considered a single valley, resulting in a small mistake in the M 0 term [20] so that we feel it is worth treating here. At zero order, the LLs (0) n,λ,ξ are solution of : The large n behavior of (0) n is the same as that of the exact LLs n,λ,ξ given in eq. (15). The major difference is that (0) n is valley independent and therefore for λξ = −1 the n = 0, 1 levels are not degenerate, moreover they appear field dependent. More quantitatively, to this order the deviation from exact LLs is at fixed nB. This gives , which up to the overall sign agrees with [1]. At first order, the LLs (1) n,λ,ξ are solution of with M 0 = −λξ [1]. To this order, the LLs acquire a valley index dependency such that for λξ = −1 the n = 0 and n = 1 levels are now degenerate, however they still keep an unphysical field dependence. To this order the deviation from exact LLs is: n,λ,ξ ) 2 = − This gives (1) n,λ,ξ − n,λ,ξ = λ , which disagrees with [1]. It is remarkable that the structure (n + 1 2 + λξ) 2 = n 2 + n(2λξ + 1) + O(n 0 ) present in eq. (67) agrees with the exact LLs (15) of structure (n + 1 + λξ)(n + λξ) = n 2 + n(2λξ + 1) + O(n 0 ) in the semi-classical limit, i.e. when n → ∞. In contrast, the LLs at zeroth order, see eq. (65), have a structure (n + 1 2 ) = n 2 + n + O(n 0 ), which does not feature the valley index ξ dependence.
Summarizing, this example shows that the deviation δ (p) is expected to be of order B p+1 when B → 0 at fixed nB. The fact that δ (2) turns out to be of order B 4 rather than B 3 is due to the vanishing of R 3 (see section III.B). The main difference between the exact and the approximate LLs concerns the n = 0 and 1 levels. The exact n = 0, 1 levels are degenerate and field independent, whereas the approximate ones become degenerate at high enough order but remain field dependent at any order. The fact that (p) n − n is generally of order B p+1 is also discussed in [1].

Generalities
Magnetic oscillations occur in different physical quantities in a metal when the magnetic field or the electron density is varied. They are due to the quantization of closed cyclotron orbits into LLs and their traversal by the Fermi energy [6]. Most notably they occur in the longitudinal magnetoresistance (as found by Shubnikov and de Haas) [43], in the magnetization (de Haas-van Alphen oscillations) [44], in the density of states, which can be measured by scanning tunneling microscopy (see e.g. [45]) or in a quantum capacitance (see e.g. [46]), and in other quantities [7].
Magnetic oscillations are typically analyzed assuming either that the chemical potential is fixed or that the number of charge carriers is fixed [47], in which case the chemical potential oscillates with the magnetic field [48]. Experimental systems are closer to one or the other limit depending on the precise setup (isolated sample, contacted sample, substrate, etc). In the following, for simplicity, we assume that the chemical potential is fixed and analyze magnetic oscillations under this assumption, having mainly DoS oscillations in mind.
Magnetic oscillations can be decomposed in harmonics [7] (see also below). In the low field limit, when the cyclotron frequency is smaller than the temperature or the disorder broadening of LLs, the oscillations in the DoS are dominated by a fundamental harmonic. The latter has a sinusoidal shape modulated by an amplitude that decays exponentially with the temperature and the disorder broadening as described by the Lifshitz-Kosevich [8] and the Dingle reduction factors [7]. The usual way to analyze the oscillations is to index the maxima (or minima, see Appendix C) of the sinusoid by an integer n and to plot this integer as a function of the inverse magnetic field 1/B. This is known as a Landau plot or as a Landau fan diagram, see Fig. 8. In simple cases, the oscillations are periodic in 1/B and this plot is well fitted by a straight line: The slope gives the oscillation frequency F ( ) in teslas and the intercept at 1/B → 0 gives the zero-field phase-shift γ 0 , which in simple cases is either a half-integer ("normal electrons") or an integer ("Dirac electrons"). In more complicated cases, it appears that the plot is not a straight line meaning that the oscillations are no-longer periodic in 1/B and that the value found for the zero-field phase-shift is neither an integer nor half an integer and is therefore hard to interpret [4].
The theoretical approach to a Landau plot starts either from the exact LLs, which are rarely known, or more often from a semiclassical quantization condition as pioneered by Onsager and Lifshitz [2]. For a detailed account of magnetic oscillations in two-dimensional electron gases, see Champel and Mineev [49]. In particular, this reference discusses consequences on magnetic oscillations in thermodynamic quantities due to fixing the number of electrons rather than the chemical potential. For a general reference on magnetic oscillations, see the classic book by Shoenberg [7].

Aperiodic magnetic oscillations and non-linear Landau plot
Here, we will analyze the Landau plot starting from the Roth-Gao-Niu quantization condition, which appears to be more general than other semiclassical quantization conditions, and resolves certain issues especially concerning the phase of magnetic oscillations and the deviation from linearity of the Landau plot. For simplicity, we assume that there is a single species (valley, spin, etc.) here.
where S( ) = (2π) 2 N 0 ( ) is the reciprocal space area of a closed cyclotron orbit at energy and γ is a phase-shift given by: This shows that the phase-shift γ is in general a function of both the Fermi energy and the magnetic field B (see also Wright and McKenzie [4]) and that it has an expression as a power series in the field B. gives: This means that the magnetic oscillations are no longer 1/B-periodic because of the B dependence of the phase-shift γ( , B). If γ( , B) = γ 0 ( ) does not depend on the magnetic field, then the oscillations are strictly 1/B-periodic with a period 1/N 0 ( ) but harmonics p = 1, .., ∞ are phase shifted by 2πpγ 0 ( ), see also [18]. It is also possible to integrate the above equation to obtain the IDoS including magnetic oscillations: For example, in the case of electrons in the continuum, we have N 0 ( ) = m 2πB and γ( , B) = 1 2 and we recover equation (92). Inverting equation (72), we obtain [4]: The first term in equation (76) (the slope of n as a function of B −1 ) gives the frequency F of magnetic oscillations as: The frequency of oscillations is directly related to the density of carriers n c as F ( ) = φ 0 n c where φ 0 = h/e is the flux quantum (if there were several species, the frequency would be F ( ) = φ 0 n c /n s where n s is the number of species). The second term gives the phase-shift in the zero field limit [1]: The third term gives the deviation from linearity (or curvature) in the Landau plot: In the end, the expansion reads which is a convenient expression to fit experiments on magnetic oscillations [4,50] (in ref. [4], |C( )| . The deviations from linearity are actually even richer in principle, see equation (9). The general structure being that of a series in power of B in addition to the usual 1/B term: n = ∞ p=−1 c p B p . We also note that, apart from the frequency of oscillations that only depends on the zerofield energy spectrum, all other terms involve the geometry of Bloch states (or equivalently the energy spectrum in a finite field).

Phase of magnetic oscillations: Berry phase and average orbital magnetic moment
The second term of the above equation (76) for n as a function of and B can be easily related to a usual expression in terms of the Berry phase and of the orbital magnetic moment (see also [1]). We use the expression of the magnetization as a function of the Berry curvature and of the orbital magnetic moment (this is the expression of the so-called "modern theory of magnetization", see [10] for review and also [51]) where here α is a band index, k is a Bloch wavevector in the first Brillouin zone, M α (k) is the orbital magnetic moment and Ω α (k) is the Berry curvature. The Heaviside step function comes from the zero temperature limit of the Fermi function n F ( α ) = 1 e ( α− )/T +1 . Differentiating with respect to the chemical potential where is the Berry phase [52] (for an iso-energy contour of energy ) and M ≡ ρ 0 ( ) −1 α,k M α (k)δ( − α (k)) is the average of the orbital magnetic moment over the Fermi surface. Therefore: Equation (83) is an important result, which was already obtained in several works, see Refs. [1,16,18,51,53,54]. The three terms in γ 0 ( ) are: the Maslov index 1/2 due to two caustics in the cyclotron orbit, the Berry phase Γ( ) and the Wilkinson-Rammal (WR) phase M ρ 0 ( ) due to the orbital magnetic moment. The last two terms occur at order in the semiclassical quantization. When inversion symmetry is present, one recognizes the often quoted result γ 0 = 1 2 − Γ 2π between the phase of magnetic oscillations and the Berry phase [55]. It is then often assumed that γ 0 can only take two values, either 1/2 for normal electrons (with zero Berry phase Γ = 0) or 0 for massless Dirac electrons (with a quantized Berry phase of Γ = π). However, in general, the Berry phase Γ( ) depends on the energy and is not quantized, see e.g. [16]. In addition, there is a second term in γ 0 , the WR phase, involving the orbital magnetic moment. Note that these two terms (Berry phase and WR phase) both involve the cell-periodic Bloch wavefunctions and can not be obtained from the zero-field energy spectrum alone. The phase of magnetic oscillations therefore crucially depends on the geometry of the Bloch states. We now give an alternative derivation of the above equation for γ 0 ( ) following [16]. The quantization condition can be written as including only the Berry phase correction at the price of shifting the energy of the cyclotron orbit by a Zeeman-like term −MB due to the orbital magnetic moment: Expanding the l.h.s. at first order in B, we obtain as ρ 0 ( ) = N 0 ( ) and S( ) = 4π 2 N 0 ( ). Therefore which agrees with the above expression (83) for γ 0 ( ). In order to illustrate the above discussion of the phase of magnetic oscillations, we present two simple examples. First, in the case of doped graphene ( > 0), Γ( ) = π does not depend on , M = 0 and γ( , B) = 0. Indeed, the Berry curvature is delta-like localized at K and K points of the first Brillouin zone (it is like a very thin flux tube carrying a π flux, i.e. half a flux quantum), and so is the orbital magnetic moment as M α (k) = 2π α (k)Ω α (k) in a twoband model. Therefore, at finite doping > 0, the average of the orbital magnetic moment over the Fermi surface indeed vanishes M = 0 but not the Berry phase Γ( ) = π [16].
Second, in the case of boron nitride or gapped graphene, there is a staggered on-site energy ±∆ that breaks the inversion symmetry of the honeycomb lattice so that the Berry curvature and the orbital magnetic moment are no longer delta-like. We consider a finite doping > ∆ with the Fermi level in the conduction band. One has an energy dependent Berry phase Γ( ) = π(1− ∆ ), an average orbital magnetic moment M = − π /v 2 ∆ which is non-vanishing (as the Berry curvature is Ω + = − ∆v 2 2 3 + ) and the density of states is = 0 despite the fact that both the Berry phase and the orbital magnetic moment depend on [16].
In some simple cases, the zero-field phase-shift γ 0 is indeed quantized to either 0 or 1/2 (modulo 1) independently of the energy. This is the case in the above examples of graphene, of boron nitride or of a graphene bilayer, for which it turns out that γ 0 = 1 2 − W 2 is independent of the energy and given by a winding number W [16]. However, this is not generally the case even in the presence of particle-hole symmetry [4], as explained in [17].

From a measurement of magnetic oscillations to response functions
The measurement of magnetic oscillations as a function of the Fermi energy could in principle be used to obtain the (per species) magnetization and the susceptibility by integration [1]. One starts from a measurement of the position of maxima or minima in a physical quantity featuring magnetic oscillations to realize a Landau plot. This plot is then fitted by a relation such as and compared with equation (9). From the frequency c −1 ( ), one obtains the density of states by integration: ρ 0 ( ) = d c −1 ( ) + const. From the constant c 0 ( ), one obtains the magnetization by integration: M 0 ( ) = d 1 2 + c 0 ( ) + const. From the coefficient c 1 ( ), one obtains the susceptibility by integration: χ 0 ( ) = 2 d c 1 ( ) + const. Etc.
The main results in this section are first that the Landau plot is in general not a straight line but can be curved. This curvature is related to a response function χ 0 as given in equation (79). Second, the zero-field phase-shift extracted from a Landau plot not only depends on the Berry phase but also on the orbital magnetic moment (the so-called Wilkinson-Rammal phase) as given in equation (83).

Conclusion
The quantization condition recently proposed by Gao and Niu is a powerful generalization of the Onsager relation. It allows one to recover various known results scattered in the literature but also to establish new results. In this paper, we have explored several consequences of this generalized Onsager condition, which we call Roth-Gao-Niu relation. In particular, we have shown how to analyze magnetic oscillations via Landau plots that are not necessarily linear and how to properly extract geometric quantities such as the Berry phase. Also, on several examples, we have shown how to improve analytically on semiclassical Landau levels. Or in the opposite direction, we have shown how to easily obtain information of magnetic response functions from exactly known Landau levels. We have also illustrated the limitation of the method on several examples in which singular behaviors such as Dirac delta functions or Heaviside step functions can not be captured.
and Niu [58] wrote a quantization condition which includes both the Berry phase and the Wilkinson-Rammal phase. The fifth term was also found by Roth [3], although in a less compact form.
Generally speaking, the semiclassical quantization of matrix-valued classical hamiltonians (for example a Bloch Hamiltonian H(k x , k y )) leads to the appearance of extra phases in the Bohr-Sommerfeld equation. Those phases are the Berry phase (related to the Berry curvature) and the Wilkinson-Rammal phase (related to the orbital magnetic moment). This is nicely discussed in Gat and Avron [51], that refer to the general work of Littlejohn and Flynn [60].
In the semiclassical approximation, the phase of the wavefunction is given by the reduced classical action A cl divided by , where l B = eB is the magnetic length [in this paragraph we restore units of and e]. For a closed cyclotron orbit, A cl = S( )l 2 B , where l B = eB is the magnetic length. Imposing that the wave-function be single-valued along a closed orbit leads to the quantization condition: The first term (nh ∼ A cl ) is classical, i.e. of order 0 . The second (Maslov), third (Berry) and fourth (Wilkinson-Rammal) terms are of order 1 (e.g. Maslov is π ). And the fifth term is of order 2 . The Maslov index term occurs already for a scalar wave-function and is related to turning points (caustics). The other terms appear only for multi-component wave-functions.
C Maxima/minima in magneto-conductance/resistance, density of states and magnetization The longitudinal magneto-conductivity generally has its maxima that coincide with that of the density of states, see e.g. [61,62] and pages 153-155 in Shoenberg's book [7]. Counterintuitively, this also corresponds to maxima in the longitudinal magneto-resistivity. This follows from the fact that the longitudinal resistivity is obtained by inverting the conductivity tensor and that the later has magnetic oscillations both in diagonal (longitudinal) and offdiagonal (Hall) elements [62]. Maxima in the DoS occur when the Fermi energy is exactly at the center of a LL n , i.e. when a LL is half-filled. Therefore = n or in other words when For the magnetization, the minima and maxima are shifted by a quarter period with respect to the density of states. See the discussion in Wright and McKenzie [4].

D Onsager quantization for electrons and for holes
The Onsager quantization condition for electrons reads S( ) = (2π) 2 B(n + 1 2 ) where S( ) is the k-space area of a closed cyclotron orbit of electron type. For holes, it is the k-space area of the hole orbit S h ( ) which is quantized. Therefore it reads: where S tot is the k-space area when the band is full of electrons.