Baptiste Savoie
email: baptiste.savoie@gmail.com

On the atomic orbital magnetism: A rigorous derivation of the Larmor and Van Vleck contributions

Keywords: -2010 number: 75.20.-g, 51.60.+a, 75.20.Ck, 75.30.Cr. MSC-2010 number: 81Q10, 82B10, 82B21, 82D05, 82D20, 82D40 Orbital magnetism, Atomic magnetism, Zero-field susceptibility, Langevin formula, Larmor diamagnetism, Van Vleck paramagnetism, Tight-binding approximation, Geometric perturbation theory, Gauge invariant magnetic perturbation theory

The purpose of this paper is to rigorously investigate the orbital magnetism of core electrons in 3-dimensional crystalline ordered solids and in the zero-temperature regime. To achieve that, we consider a non-interacting Fermi gas subjected to an external periodic potential modeling the crystalline field within the tight-binding approximation (i.e. when the distance between two consecutive ions is large). For a fixed number of particles in the Wigner-Seitz cell and in the zero-temperature limit, we derive an asymptotic expansion for the bulk zero-field orbital susceptibility. We prove that the leading term is the superposition of the Larmor diamagnetic contribution, generated by the quadratic part of the Zeeman Hamiltonian, together with the 'complete' orbital Van Vleck paramagnetic contribution, generated by the linear part of the Zeeman Hamiltonian, and related to field-induced electronic transitions.

The first important contribution to the understanding of diamagnetism of ions (we assimilate an atom to an ion of charge zero) and molecules (including the polar ones) goes back at least to 1905 with the three papers [START_REF] Langevin | Sur la théorie du magnétimse[END_REF][START_REF] Langevin | Sur la théorie du magnétisme[END_REF][START_REF] Langevin | Magnétisme et théorie des électrons[END_REF] of P. Langevin. We mention all the same that, in all likelihood, W. Weber brought between 1852 and 1871 the pioneer ideas through his molecular theory of magnetism, in which he already introduced the idea that the magnetic effects are due to orbiting motion of electric charges around fixed charges of opposite sign, see [START_REF] Assis | Webers's Planetary Model of the Atom[END_REF]. The Langevin's microscopic theory essentially leans on the classical Maxwell equations of electromagnetism. Putting things back into context (the nuclear structure of atoms was experimentally discovered in 1911), Langevin considered that matter is formed by electrons in stable periodic motion, the mutual actions between electrons assuring the mechanical stability. In particular, the molecules contain at least one closed electron orbit with a fixed magnetic moment out of any external field (electrons are assimilated with particulate Ampère's currents), and the different orbits in each molecule have such a moment and such orientations, that their resultant moment may vanish, or not. With these assumptions, Langevin calculated the mean variation of the magnetic moment (orthogonal to the orbit) of electron moving in intramolecular closed orbits under the influence of an external constant magnetic field. This led to the so-called Langevin formula for the diamagnetic susceptibility per unit volume of n electrons, see [38, pp. 89]:

χ dia La = -n e 2 4mc 2 r 2 , (1.1)
where r is the distance from the molecule's centre of mass to the electron (considering its motion in the projected orbit on the plane perpendicular to the magnetic field) and r 2 is the average r 2 over all the molecular electrons. Here and hereafter, e, m, c are universal constants denoting respectively the elementary electric charge, the electron rest mass and the speed of light in vacuum. Note that (1.1) is independent of the temperature (provided that the orbits retain the same size) and independent of whether or not the initial resultant magnetic moment of the molecule is null. For completeness' sake, we mention that Langevin derived also in [START_REF] Langevin | Magnétisme et théorie des électrons[END_REF] the analytic expression of the Curie's empirical law for molecules of paramagnetic substances. The paramagnetic susceptibility per unit volume of N identical molecules with a non-zero resultant magnetic moment with magnitude M reads as, see [38, pp. 119]:

χ para La (β) = N 3 βM 2 , β := 1 k B T , (1.2)
where T is the absolute temperature and k B denotes the Boltzmann constant. Note that, in the derivation of (1.2), Langevin used the Boltzmann factor to determine the spatial distribution of the permanent magnetic moment.

On the 1910's, a series of works came up in response to the Langevin's theory. For a review, see [START_REF] Van Leeuwen | Problèmes de la théorie électronique du magnétisme[END_REF]Sec. X]. Essentially, the purpose consisted in computing, within the framework of classical statistical mechanics, the thermal average magnetic moment of various sorts of molecules taking into account specific kinds of collisions. Among these works, we mention the paper [START_REF] Van Leeuwen | Problèmes de la théorie électronique du magnétisme[END_REF] of J.H. Van Leeuwen which played a crucial role in the understanding of the origin of magnetism. In [START_REF] Van Leeuwen | Problèmes de la théorie électronique du magnétisme[END_REF]Sec. VIII], she investigated the magnetic response of the free electrons in metals within the Drude model, in which the electrons are assimilated to solid spheres and the velocity distribution is given at thermal equilibrium by the Maxwell-Boltzmann distribution. In the presence of a constant magnetic field, and assuming that the electrons collide elastically with the molecules (forming the matter) assumed to be fixed, then she showed that, whenever a state of equilibrium exists, the thermal average magnetic moment of the electron gas in an arbitrary volume element is identically zero. In other words, within the framework of classical statistical mechanics, the free electron gas displays no diamagnetic effect. This result is referred in Physics literature to as the Bohr-Van Leeuwen theorem. This apparent contradiction with the Langevin results (also stated within the classical theory) will be removed by the development of quantum mechanics: the assumptions made by Langevin (stationary of the electron orbits and permanence of the magnetic moment) are actually of a quantum nature.

In 1920, W. Pauli was interested in diamagnetism of monoatomic gases in [START_REF] Pauli | Theoretische Bemerkungen über den Diamagnetismus einatomiger Gase[END_REF]. Still within the framework of the classical theory, his approach slightly differs from the one of Langevin, in the sense that it is based on the Larmor precession theorem. This latter states that, when placed in an external constant magnetic field with intensity B, an atom with any number of electron orbits precesses around an axis through the nucleus and parallel to the magnetic field. It causes the angular velocity of an electron in a periodic orbit to be increased by -Be 2mc without the orbit undergoing any change. This led to the formula for the diamagnetic susceptibility per unit volume of N identical atoms assumed to have random orientation in space, see [42, pp. 203]:

χ dia La * = -N e 2 6mc 2 i r 2 i , (1.3)
where the sum is over all the electrons of a single atom, r 2 i is the square distance from the nucleus to the i-th electron, and r 2 i has to be understood as the time average of r 2 i . Note that the additional 2 3 -factor (see (1.1)) comes from the statistical mean (the orientations of the atoms are random). (1.3) is sometimes referred to as the Langevin formula in the form given by Pauli.

In 1927-1928, J.H. Van Vleck revisited in a series of three papers [START_REF] Vleck | On dielectric constants and magnetic susceptibilities in the new quantum mechanics Part I. A general proof of the Langevin-Debye formula[END_REF][START_REF] Vleck | On dielectric constants and magnetic susceptibilities in the new quantum mechanics. Part II-Application to dielectric constants[END_REF][START_REF] Vleck | On dielectric constants and magnetic susceptibilities in the new quantum mechanics Part III-Application to dia-and paramagnetism[END_REF] the Debye theory on dielectric constant and Langevin theory on magnetic susceptibility of ions/molecules within the framework of quantum mechanics. In particular, he derived a general formula in [51, Eq. ( 13)] for the total magnetic susceptibility per unit volume, at any 'temperature' β > 0 and in the zero-field limit, of N non-interacting randomly oriented identical ions/molecules with a non-zero (time-averaged) resultant magnetic moment. Only the contribution of electrons is considered (the contribution of nuclei is assumed to be negligible), and the interactions between electrons are disregarded. The 'proof' given by Van Vleck requires the assumption that the energy intervals between the various component levels of the low-energy states of the ions/molecules are small compared to β -1 , see also [START_REF] Vleck | On dielectric constants and magnetic susceptibilities in the new quantum mechanics Part I. A general proof of the Langevin-Debye formula[END_REF]Sec. 2] and below pp. 5. This assumption implies that the Van Vleck's results do not cover the zero-temperature regime. The formula [51, Eq. ( 13)] is a generalization of the complete classical Langevin formula (obtained by adding (1.3) and (1.2)) but including quantum effects, and consists of the sum of two contributions:

• A β-dependent contribution (second term in the r.h.s. of [START_REF] Vleck | On dielectric constants and magnetic susceptibilities in the new quantum mechanics Part III-Application to dia-and paramagnetism[END_REF]Eq. (13)]).

It is purely paramagnetic, depends linearly on β and arises from the presence of the non-zero resultant magnetic moment. Van Vleck points out that this term vanishes when the total electronic orbital and spin angular momentum of ions both are null. Moreover, its classical equivalent is (1.2).

• A β-independent contribution (first term in the r.h.s. of [START_REF] Vleck | On dielectric constants and magnetic susceptibilities in the new quantum mechanics Part III-Application to dia-and paramagnetism[END_REF]Eq. (13)]).

It is assumed to represent the diamagnetic effect, since experimentally, diamagnetism in gases is known to be temperature-independent (at constant density). Denoted by N α, it is the sum of two terms, see [START_REF] Vleck | On dielectric constants and magnetic susceptibilities in the new quantum mechanics Part III-Application to dia-and paramagnetism[END_REF]Eq. (15)]: N α = N (X vV + X La ) .

(1.4)

Here, X La is the so-called Larmor contribution which is purely diamagnetic, and N X La reduces to the Langevin formula in (1.3) in the classical limit. As for X vV , it is purely paramagnetic and has no equivalent in the classical theory. It is often referred to as (orbital) Van Vleck paramagnetic susceptibility. Besides, Van Vleck analyzed the origin of each one of these two contributions, see e.g. [START_REF] Vleck | The theory of electric and magnetic susceptibilities[END_REF]Sec. VIII.49]. Restricting to the case of a single ion, he showed that X vV and X La are generated respectively by the the linear part and the quadratic part of the Zeeman Hamiltonian of the electron gas which are defined by:

H Z,lin := µ B (L + g 0 S) • B, H Z,qua := e 2 2mc 2 i B × r i 2 2 , (1.5)
in the case where the magnetic vector potential is chosen to be in the symmetric gauge. Here and hereafter, L and S stand for the total electronic orbital and spin angular momentum respectively, g 0 and µ B for the electronic g-factor and Bohr magneton respectively, B for the external constant magnetic field and r i for the position vector of the i-th electron. From the foregoing, Van Vleck concluded that X La always exists and is in competition with the paramagnetic contribution X vV when the ion has either its total electronic spin angular momentum or orbital angular momentum different from zero in its low-energy states. Moreover, he claimed that X vV does not vanish in great generality in the case of molecules, see [52, Sec. X.69]. Hence N X La is an upper bound limit to the diamagnetism of electrons in all cases (ions/molecules).

Let us outline the 'proof' of [51, Eq. ( 13)], we refer to [51, pp. 594-598] and also [START_REF] Vleck | On dielectric constants and magnetic susceptibilities in the new quantum mechanics Part I. A general proof of the Langevin-Debye formula[END_REF]Sec. 3]. Van Vleck took into account many degrees of freedom of ions/molecules (such that the temperature rotation, the orientations relative to the magnetic field, etc...). We simplify the model and restrict to the case of a single ion. Assume for simplicity that the magnetic field B is constant and parallel to the k-th (k ∈ {1, 2, 3}) direction, i.e. B = Be k . The derivation consists of three steps. Step 1: Expressing the average induced magnetic moment with the Boltzmann distribution. Let E j (B), j ≥ 0 be the energy levels of the ion in the magnetic field. Let M k (j, B) := -dEj dB (B) be the induced magnetic moment per unit volume of the state of energy E j (B). The average value of the induced magnetic moment per unit volume M k (β, B) at 'temperature' β > 0 is defined as the thermal equilibrium average (by the Boltzmann distribution) of the M k (j, B)'s:

M k (β, B) := - 1 j≥0 e -βEj(B) j≥0 dE j (B) dB e -βEj(B) . (1.6)
Here, the energy levels are assumed to be non-degenerate.

Step 2: Expanding the E j (B)'s in powers of B. Provided that B is small enough, the asymptotic perturbation theory allows to compute the changes in the energy levels by treating the magnetic field as a perturbation. Denoting by |j the j-th state, one has, up to the second order correction (in the absence of degeneracy):

E j (B) = E j (0) + j |H Z,lin + H Z,qua | j + l =j | j |H Z,lin + H Z,qua | l | 2 E j (0) -E l (0) + • • • ,
where H Z,lin and H Z,qua are defined in (1.5). This leads to the expansion:

E j (B) = E j (0) + BE (1)
j (0) + B 2 E (2) j (0) + • • • , with : (1.7) E (1) j (0) := µ B j |(L + g 0 S) • e k | j , (1.8)

E

(2)

j (0) := e 2 8mc 2 i j (r i × e k ) 2 j + µ 2 B l =j | j |(L + g 0 S) • e k | l | 2 E j (0) -E l (0) . (1.9)
Step 3: Deriving the zero-field magnetic susceptibility. The rest consists in substituting (1.7) into (1.6), then expanding the Boltzmann factor e -βEj(B) (in the numerator and denominator) in

Taylor series in B. Assuming that the field-dependent corrections to the energy in (1.7) are smaller than β -1 , an expansion in powers of B of M k (β, B) is derived. Within the framework of the linear response theory, the zero-field magnetic susceptibility per unit volume X (β, 0) corresponds to the coefficient of the linear term in B. Since the induced magnetic moment is null in vanishing field, X (β, 0) reduces to:

X (β, 0) = 1 j≥0 e -βEj (0) j≥0 β E (1)
j (0) 2 -2E (2)
j (0) e -βEj(0) .

(1.10)

In the derivation of [51, Eq. ( 13)], the excited states (i.e. j ≥ 1) are assumed to be unoccupied. Its counterpart for the particular case we consider here is:

X (0) (β, 0) = β E (1) 0 (0) 2 -2E (2)
0 (0).

(1.11)

The β-independent part in (1.11) can be identified from (1.9) and reads as:

- e 2 4mc 2 i 0 (r i × e k ) 2 0 + 2µ 2 B l =0 | 0 |(L + g 0 S) • e k | l | 2 E l (0) -E 0 (0) . (1.
12)

The first and second term are the Larmor and Van Vleck contributions respectively. If the ion has its total electronic orbital and spin angular momentum both equal to zero, (1.12) reduces to the Larmor contribution. In accordance with Hund's rules, this is the case of an ion having all its electron shells filled.

In regard to the magnetism of ions/molecules in the zero-temperature regime, another approach is commonly encountered in Physics literature, see e.g. [START_REF] Ashcroft | Solid State Physics[END_REF]Chap. 31]. Let us give the main idea. Consider a system of N non-interacting identical ions subjected to a constant magnetic field at thermal equilibrium. For simplicity, assume that the magnetic field is given by B = Be k , k ∈ {1, 2, 3}. Let E j (B), j ≥ 0 be the energy levels of a single ion in the magnetic field. In the canonical ensemble of statistical mechanics, the zero-field magnetic susceptibility per unit volume is defined by:

X (β, N, B = 0) := - ∂ 2 F ∂B 2 (β, N, B = 0), β > 0, (1.13)
where F (β, N, B) stands for the Helmholtz free energy of the system:

F (β, N, B) := - 1 β ln   j≥0 e -N βEj(B)   . (1.14)
Based on the principle that only the ground-state of the system is occupied in the zero-temperature regime for B small enough, the Helmholtz free energy reduces to the ground-state energy: N E 0 (B).

From (1.13) along with the expansion (1.7) obtained within the asymptotic perturbation theory:

lim β↑∞ X (β, 0) = -2N E (2) 0 (0),
which holds in the absence of degeneracy. If N = 1, this is nothing but the β-independent part of the magnetic susceptibility derived in (1.12).

What are the motivations of this paper?

In the light of the works mentioned in Sec. 1.1, a formula for the zero-field magnetic susceptibility of ions in the low/high-temperature regime can be derived (at least, at the formal level) knowing the energy levels of the ion likely to be occupied, and the associated stationary states.

The Van Vleck's approach and the one in [START_REF] Ashcroft | Solid State Physics[END_REF]Chap. 31] have many similarities. Firstly, an ion is modeled as a (semi)classical object in its own right, characterized by a discrete set of energy levels. In the wake of this model, the Maxwell-Boltzmann distribution is used to derive an expression for the thermal average induced magnetic moment/magnetic susceptibility, see (1.6) and (1.14)- (1.13). The use of the Maxwell-Boltzmann distribution turns out to be essential to write down 'easy-handled' formulae for the quantities of interest. However, such formulae hold only for temperatures large enough or in the semiclassical limit. Quantum mechanics is used only at the second stage in order to compute the changes in the energy levels by treating the magnetic field as a perturbation, see (1.7). Subsequently, when computing the zero-field magnetic susceptibility, only the contribution involving the corrections to the ground-state energy (i.e. which couple the ground-state energy with the excited-states energies) are considered, whereas the contribution involving the corrections to the excited-states energies (i.e. which couple the energy levels between them) are discarded, see the transition from (1.10) to (1.11). Such a simplification is justified by the fact that the excited-states energies of an ion are not occupied in a regime of weak magnetic field or low temperature, and actually dodges the questions of convergence of the series.

Treating an ion as a whole, and retaining only the energy levels which may be occupied dodge the difficulty to compute directly the magnetic susceptibility of the core electrons within the framework of quantum statistical mechanics, whose an expression is more delicate to derive. Indeed, for a wide class of potentials (including Coulomb) used to model the interaction electron-nucleus, the counterpart of the quantity (1.6) for the core electrons is divergent. As a result, some natural question arise: do the results mentioned in Sec. 1.1 still hold true if one computes the magnetic susceptibility of core electrons within the framework of quantum statistical mechanics? Does the susceptibility representing the orbital effects still consist of two contributions? If so, how they differ from the 'usual' Larmor and Van Vleck contributions? This paper takes place in this direction.

The aim of this paper is to rigorously revisit the atomic orbital magnetism (i.e. we focus on the magnetic effects which do not arise from spin effects) in the zero-temperature regime, the most complicated situation. Our approach is substantially different from the ones mentioned in Sec. 1.1. To model the core electrons of an ion, we consider a non-interacting Fermi gas subjected to an external periodic potential (modeling an ideal lattice of fixed nuclei) within the tight-binding approximation. Under this approximation, we suppose that the distance R > 0 between two consecutive ions is sufficiently large so that the Fermi gas 'feels' mainly the potential energy generated by one single nucleus. To investigate the atomic orbital magnetism, our starting-point is the expression derived in [START_REF] Briet | A rigorous approach to the magnetic response in disordered systems[END_REF]Thm. 1.2] for the bulk zero-field orbital susceptibility of the Fermi gas valid for any 'temperature' β > 0 and R > 0. We emphasize that this expression is derived from the usual rules of quantum statistical mechanics. Our main result is Theorem 1.1, and Remark 1.4 makes the connection with the works mentioned in Sec. 1.1.

In the framework of Mathematical-Physics, the rigorous study of orbital magnetism and more generally of diamagnetism, have been the subject of numerous works. Let us list the main ones among them. The first rigorous proof of the Landau susceptibility formula for free electron gases came as late as 1975, due to Angelescu et al. in [START_REF] Angelescu | On the Landau diamagnetism[END_REF]. Then in 1990, Helffer et al. developed for the first time in [START_REF] Helffer | On diamagnetism and the de Haas-van Alphen effect[END_REF] a rigorous theory based on the Peierls substitution and considered the connection with the diamagnetism of Bloch electrons and the de Haas-Van Alphen effect. These and many more results were reviewed in 1991 by Nenciu in [START_REF] Nenciu | Dynamics of band electrons in electric and magnetic fields: Rigorous justification of the effective Hamiltonians[END_REF]. In 2001, Combescure et al. recovered in [START_REF] Combescure | Rigorous semiclassical results for the magnetic response of an electron gas[END_REF] the Landau susceptibility formula in the semiclassical limit. In 2012, Briet et al. gave for the first time in [START_REF] Briet | A rigorous proof of the Landau-Peierls formula and much more[END_REF] a rigorous justification of the Landau-Peierls approximation for the susceptibility of Bloch electron gases. Finally we mention the following papers [START_REF] Fournais | The magnetisation of large atoms in strong magnetic fields[END_REF][START_REF] Fournais | On the semiclassical asymptotics of the current and magnetic moment of a noninteracting electron gas at zero temperature in a strong constant magnetic field[END_REF][START_REF] Fournais | On the total magnetic moment of large atoms in strong magnetic fields[END_REF] in connection with atomic magnetism.

The setting and the main result.

Consider a three-dimensional quantum gas composed of a large number of non-relativistic identical and indistinguishable particles, with charge q = 0 and mass m = 1, obeying the Fermi-Dirac statistics, and subjected to an external constant magnetic field. Furthermore, each particle interacts with an external periodic electric potential modeling the ideal lattice of fixed ions in crystalline ordered solids. The spin of particles is not taken into consideration since we are only interested in orbital effects arising from the 'motion' of particles in the medium (the spin-orbit coupling is disregarded). Moreover, the interactions between particles are neglected (strongly diluted gas assumption) and the gas is at equilibrium with a thermal and particles bath.

Let us make our assumptions more precise. The gas is confined in a cubic box Λ L := (-L 2 , L 2 ) 3 with L > 0, centered at the origin of coordinates. We denote its Lebesgue-measure by |Λ L |. We consider a uniform magnetic field B := Be 3 , e 3 := (0, 0, 1) parallel to the third direction of the canonical basis of R 3 . We choose the symmetric gauge, i.e. the magnetic vector potential is defined by A(x) := 1 2 B × x = Ba(x), a(x) := 1 2 (-x 2 , x 1 , 0) so that B = ∇ × A(x) and ∇ • A(x) = 0. Hereafter, we denote by b := q c B ∈ R the cyclotron frequency. The potential energy modeling the interaction between each particle and the ideal lattice of fixed nuclei is given by:

V R := υ∈Z 3 u(• -Rυ), R > 0, (1.15)
where the single-site potential u satisfies the following assumption:

(A r ) u ∈ C 1 (R 3 ; R) is compactly supported.
We denote by Ω R , R > 0 the Wigner-Seitz cell of the RZ 3 -lattice centered at the origin of coordinates. Its Lebesgue-measure is denoted by |Ω R |. Introduce the one-particle Hamiltonian. On C ∞ 0 (Λ L ), define ∀R > 0:

H R,L (b) := 1 2 (-i∇ -ba) 2 + V R , b ∈ R. (1.16)
By standard arguments, ∀R > 0 and ∀b ∈ R, (1.16) extends to a family of self-adjoint and semibounded operators for any L ∈ (0, ∞), denoted again by H R,L (b), with domain

H 1 0 (Λ L ) ∩ H 2 (Λ L
). This definition corresponds to choose Dirichlet boundary conditions on ∂Λ L . Moreover, since the inclusion

H 1 0 (Λ L ) ֒→ L 2 (Λ L ) is compact, then H R,L (b
) has purely discrete spectrum with an accumulation point at infinity. We denote by {λ When Λ L fills the whole space, on C ∞ 0 (R 3 ) define ∀R > 0:

H R (b) := 1 2 (-i∇ -ba) 2 + V R , b ∈ R. ( 1
N R (E) := lim L↑∞ N R,L (E) |Λ L | = lim L↑∞ Tr L 2 (R 3 ) χ ΛL P (-∞,E) (H R (b))χ ΛL |Λ L | , (1.18)
where χ ΛL is the multiplication operator by the characteristic function of Λ L .

Let us now define some quantities related to the confined Fermi gas introduced above within the framework of quantum statistical mechanics. In the grand-canonical situation, let (β, z, |Λ L |) be the fixed external parameters. Here β := (k B T ) -1 > 0 (k B stands for the Boltzmann constant) is the 'inverse' temperature and z := e βµ > 0 (µ ∈ R stands for the chemical potential) is the fugacity. For any β > 0, z > 0 and b ∈ R, the finite-volume pressure and density are respectively defined ∀R > 0 by, see e.g. [START_REF] Huang | Statistical Mechanics, Second Edition[END_REF][START_REF] Angelescu | On free quantum gases in a homogeneous magnetic field[END_REF][START_REF] Angelescu | On the Landau diamagnetism[END_REF]:

P R,L (β, z, b) := 1 β|Λ L | Tr L 2 (ΛL) ln ½ + ze -βHR,L(b) = 1 β|Λ L | ∞ j=1 ln 1 + ze -βλ (j) R,L (b) , (1.19) ρ R,L (β, z, b) := βz ∂P R,L ∂z (β, z, b) = 1 |Λ L | ∞ j=1 ze -βλ (j) R,L (b) 1 + ze -βλ (j) R,L (b)
.

(

X R,L (β, z, b) := q c 2 ∂ 2 P R,L ∂b 2 (β, z, b), β > 0, z > 0, b ∈ R, R > 0.
When Λ L fills the whole space (i.e. in the limit L ↑ ∞), the thermodynamic limits of the three grand-canonical quantities defined above generically exist, see e.

ρ R (β, z, b) := βz ∂P R ∂z (β, z, b) = lim L↑∞ βz ∂P R,L ∂z (β, z, b), (1.21)
X R (β, z, b) := q c 2 ∂ 2 P R ∂b 2 (β, z, b) = lim L↑∞ q c 2 ∂ 2 P R,L ∂b 2 (β, z, b), (1.22)
and the convergences are compact w.r.t. (β, z, b). The limit commutes with the first derivative (resp. the second derivative) of the pressure w.r.t. the fugacity z (resp. the cyclotron frequency b). Moreover, ∀β > 0 and ∀R > 0, P R (β, • , • ) is jointly smooth in (z, b) ∈ (0, ∞) × R, see e.g. [START_REF] Savoie | [END_REF].

Next, we switch to the canonical conditions and assume that the density of particles ρ 0 > 0 becomes, in addition with the 'inverse' temperature, a fixed external parameter. Seeing the bulk density in (1.21) as a function of the µ-variable, denote by µ

(0) R (β, ρ 0 , b) ∈ R, the unique solution of the equation: ρ R (β, µ, b) = ρ 0 , β > 0, b ∈ R, R > 0.
The inversion of the relation between the bulk density and the chemical potential is ensured by the fact that ∀β > 0, ∀b ∈ R and ∀R > 0, µ → ρ R (β, µ, b) is a strictly increasing function on R, and actually it defines a C ∞ -diffeomorphism of R into (0, ∞), see e.g. [START_REF] Savoie | [END_REF][START_REF] Briet | A rigorous proof of the Landau-Peierls formula and much more[END_REF]. In the following, we consider the situation in which the density of particles is given by:

ρ 0 (R) = n 0 |Ω R | , n 0 ∈ N * , R > 0, (1.23)
where n 0 plays the role of the number of particles in the Wigner-Seitz cell. From (1.22) and under the above conditions, the bulk zero-field orbital susceptibility at fixed β > 0 and density ρ 0 (R), R > 0 is defined by:

X R (β, ρ 0 (R), b = 0) := X R β, µ (0) R (β, ρ 0 (R), b = 0), b = 0 . (1.24)
Before giving our main theorem related to the quantity in (1.24), let us introduce some reference operators which will be involved in its statements.

On C ∞ 0 (R 3 ), define the 'single atom' Schrödinger operator:

H P := 1 2 (-i∇) 2 + u, (1.25)
where u is the function appearing in (1.15) and obeying assumption (A r ). By standard arguments, H P is essentially self-adjoint and its self-adjoint extension, denoted again by H P , is semi-bounded with domain H 2 (R 3 ). Moreover, σ ess (H P ) = [0, ∞) is absolutely continuous, and H P has finitely many eigenvalues in (-∞, 0) if any, see e.g. [START_REF] Reed | Methods of Modern Mathematical Physics, IV : Analysis of Operators[END_REF]Thm. XIII.15]. Throughout, we suppose:

(A m ) H P has at least one eigenvalue in (-∞, 0), together with the non-degeneracy assumption:

(A nd ) All the eigenvalues of H P in (-∞, 0) are non-degenerate.

Hereafter, we denote by {λ l } τ l=1 , τ ∈ N * the set of eigenvalues of H P in (-∞, 0) counting in increasing order, and by {Φ l } τ l=1 the set of corresponding normalized eigenfunctions. Also, we denote by Π l the orthogonal projection onto the eigenvector Φ l , and define Π ⊥ l := ½ -Π l . Since Π l commutes with H P , one has the decomposition in direct sum:

L 2 (R 3 ) = Π l L 2 (R 3 )⊕Π ⊥ l L 2 (R 3
). In the presence of the uniform magnetic field (as in (1.17)), define on C ∞ 0 (R 3 ), the 'single atom' magnetic Schrödinger operator:

H P (b) := 1 2 (-i∇ -ba) 2 + u, b ∈ R. ( 1
H P = -ba • (-i∇) + b 2
2 a 2 provided that it is weak. From the asymptotic perturbation theory in [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF]Sec. VIII] and due to the assumption (A nd ), then there exists b > 0 s.t. ∀|b| ≤ b and for any l ∈ {1, . . . , τ }, H P (b) has exactly one and only one eigenvalue λ l (b) near λ l which reduces to λ l in the limit b → 0. In particular, for any l ∈ {1, . . . , τ }, λ l (• ) can be written in terms of an asymptotic power series in b, see e.g. [START_REF] Briet | Locating the spectrum for magnetic Schrödinger and Dirac operators, Commun[END_REF]Thm. 1.2]. Hereafter we denote by {λ l (b)} τ l=1 , τ ∈ N * and |b| ≤ b, the set of these eigenvalues counting in increasing order.

Getting back to the 'single atom' operator in (1.25), by [47, Thm. B.7.2], for any ξ ∈ C \ [inf σ(H P ), ∞), the resolvent operator (H Pξ) -1 is an integral operator with integral kernel (H Pξ) -1 (• , • ) jointly continuous on R 6 \ D, D := {(x, y) ∈ R 6 : x = y}. Introduce on L 2 (R 3 ) the operators T P,j (ξ), j = 1, 2 generated via their kernel respectively defined on R 6 \ D as:

T P,1 (x, y; ξ) := a(x -y) • (i∇ x )(H P -ξ) -1 (x, y), (1.27)
T P,2 (x, y; ξ) := 1 2 a 2 (xy)(H Pξ) -1 (x, y).

(1.28)

Here is our main result related with the quantity defined in (1.24):

Theorem 1.1 Suppose that the assumptions (A r ), (A m ) and (A nd ) hold.

Assume that the number of particles n 0 ∈ N * in the Wigner-Seitz cell is fixed and satisfies n 0 ≤ τ , while the density is given by (1.23). Then: (i). For any 0 < α < 1, there exists a R-independent constant c > 0 s.t.

X R (ρ 0 (R)) := lim β↑∞ X R (β, ρ 0 (R), b = 0) = 1 |Ω R | X P (n 0 ) + O e -cR α , (1.29)
with:

X P (n 0 ) := - q c 2 i π Tr L 2 (R 3 ) Γn 0 dξ ξ(H P -ξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] , (1.30)
where Γ n0 is any positively oriented simple closed contour enclosing the n 0 smallest eigenvalues of H P , while letting the rest of the spectrum outside.

(ii). The R-independent quantity in (1.30) can be identified with:

X P (n 0 ) = - q c 2 n0 l=1 d 2 λ l db 2 (b = 0). (1.31) (iii).
(1.30) can be rewritten as a sum of two contributions:

1 |Ω R | X P (n 0 ) = 1 |Ω R | X La (n 0 ) + 1 |Ω R | X vV (n 0 ), (1.32)
with:

1 |Ω R | X La (n 0 ) := - q c 2 1 4|Ω R | n0 l=1 Φ l , X 2 1 + X 2 2 Φ l L 2 (R 3 ) , (1.33) 1 |Ω R | X vV (n 0 ) := q c 2 1 2|Ω R | n0 l=1 L 3 Φ l , Π ⊥ l (H P -λ l ) Π ⊥ l -1 L 3 Φ l L 2 (R 3 )
.

(1.34)

Here, X k := X • e k , k ∈ {1, 2, 3} stands for the position operator projected in the k-th direction, and

L k := L • e k the k-th component of the orbital angular momentum operator L := X × (-i∇).
Remark 1.2 (R-behavior in the asymptotic (1. 29)) The leading term decreases as R -3 , the remainder (explicitly identified in the proof ) decreases exponentially in R α , 0 < α < 1.

Remark 1.3 (Rewriting of (1.34) in a sum of two contributions) Denote by Π ac the orthogonal projection onto the absolutely continuous spectrum of the 'single atom' operator H P . Since Π l , l ∈ {1, . . . , τ } and Π ac commute with H P , then one has the decomposition in direct sum:

Π ⊥ l L 2 (R 3 ) = τ m=1 m =l Π m L 2 (R 3 ) ⊕ Π ac L 2 (R 3 ).
Due to this decomposition, the reduced resolvent operator appearing in (1.34) can be rewritten (below, we use the bra-ket notation) as:

Π ⊥ l (H P -λ l ) Π ⊥ l -1 = τ m=1 m =l |Φ m Φ m | λ m -λ l + (Π ac (H P -λ l ) Π ac ) -1 .
In view of (1.34), one has under the same conditions the rewriting:

1 |Ω R | X vV (n 0 ) = 1 |Ω R | X d vV (n 0 ) + 1 |Ω R | X ac vV (n 0 ), (1.35)
with:

X d vV (n 0 ) : = q c 2 1 2 n0 l=1 τ m=1 m =l | Φ m , L 3 Φ l | 2 λ m -λ l =      q c 2 1 2 n0 l=1 τ m=n0+1 | Φ m , L 3 Φ l | 2 λ m -λ l , if n 0 < τ , 0, if n 0 = τ ; (1.36) X ac vV (n 0 ) := q c 2 1 2 n0 l=1 L 3 Φ l , (Π ac (H P -λ l ) Π ac ) -1 L 3 Φ l L 2 (R 3 ) , (1.37)
and the quantity in (1.37) can be expressed in terms of the generalized eigenfunctions of H P by:

(Π ac (H P -λ l ) Π ac ) -1 = [0,∞) dE S 2 dω |Ψ E,ω Ψ E,ω | E -λ l .
Hence, (1.34) is the superposition of two contributions: (1.36) involving the corrections to the isolated eigenvalues which couple the isolated eigenvalues between them, and (1.37) involving the corrections to the isolated eigenvalues coupling with the absolutely continuous part of the spectrum.

Remark 1.4 (Connection with the Van Vleck results

) Under the condition of large separation of ions and in the zero-temperature limit, the leading term of the bulk zero-field orbital susceptibility consists of the superposition of two contributions, see (1.32).

• Contribution in (1.33).

It is purely diamagnetic and generated by the quadratic part of the Zeeman Hamiltonian. Still assuming one single ion in the Wigner-Seitz cell, it reduces in the classical limit to the Langevin formula, see (1.1):

- q c 2 1 4|Ω R | n0 l=1 r 2 l , (1.38)
where r l is the distance from the origin to the l-th particle in the plane orthogonal to e 3 (as a rule, the centre of mass of the ion nucleus is at the origin, and the fixed axis passing through it, is taken parallel to the magnetic field).

• Contribution in (1.34).

In view of the rewriting (1.35) with (1.36)-(1.37), it is purely paramagnetic and generated by the linear part of the Zeeman Hamiltonian. Arising from the quantum description, it has no classical equivalent. It is related to field-induced electronic transitions.

Let us turn to the comparison with the Van Vleck results in Sect. 1.1. The contribution in (1.33) exhibits the same features than X La in (1.4), and its counterpart within the Van Vleck formulation corresponds to the first term in (1.12). Then, it identifies with the diamagnetic Larmor susceptibility. The contribution in (1.34) exhibits the features of X vV in (1.4), but to identify which quantity from the decomposition (1.35) corresponds to the second term in (1.12) within the Van Vleck formulation, we need the following. As pointed out below (1.12), when dealing with ions, the orbital Van Vleck paramagnetic susceptibility vanishes when the ion has all its electron shells filled, in accordance with Hund's rules. In a way, this condition (i.e. the electron shells filled) corresponds to have n 0 = τ within our model, i.e. the number of negative isolated eigenvalues coincides with the number of electrons in the Wigner-Seitz cell. Since (1.36) vanishes when n 0 = τ , then its counterpart within the Van Vleck formulation corresponds to the second term in (1.12). From the foregoing, this means that, within our model, the quantity in (1.37) is an additional contribution to the 'usual' orbital Van Vleck susceptibility. Involving the corrections to the isolated negative eigenvalues coupling with the absolutely continuous part of the spectrum, it has no counterpart within the Van Vleck formulation. Due to this feature, (1.34) represents the 'complete' orbital Van Vleck paramagnetic susceptibility. In that sense, Theorem 1.1 is a generalization of the Van Vleck results.

Outline of the proof of Theorem 1.1.

The starting-point is the expression (3.9) for the bulk zero-field orbital susceptibility (under the grand-canonical conditions) which holds for any β > 0 and R > 0. The derivation of such an expression is the main subject of [START_REF] Briet | A rigorous approach to the magnetic response in disordered systems[END_REF]. Its proof relies on the gauge invariant magnetic perturbation theory applied on the magnetic resolvent operator. It allows to keep a good control on the linear growth of the potential vector, see below for further details.

The difficulties coming up when proving Theorem 1.1 are threefold: [START_REF] Angelescu | On the Landau diamagnetism[END_REF]. Isolating the main R-dependent contribution (which decays polynomially in R) and a remainder term (which decays exponentially in R α ) from (3.9) when switching to the canonical conditions (the density of particles is fixed, and given by (1.23)). [START_REF] Angelescu | On free quantum gases in a homogeneous magnetic field[END_REF]. Performing the zero-temperature limit, while taking into account the location of the Fermi energy defined in (3.1). (3). Identifying the leading term of the asymptotic expansion in the zero-temperature regime with the corresponding relevant physical quantities.

• Outline of the proof of (i).

The proof of the asymptotic expansion in (1.29) requires three steps.

Step 1: Isolating the main R-dependent contribution -Part 1. The formula in (3.9) involves the trace per unit volume:

|Ω R | -1 Tr L 2 (R 3 ) χ ΩR (H R -ξ) -1 [T R,1 (ξ)T R,1 (ξ) -T R,2 (ξ)] χ ΩR , (1.39)
where χ ΩR denotes the operator multiplication by the indicator function of the Wigner-Seitz cell Ω R , and T R,j (ξ), with j = 1, 2 involve the resolvent operator (H Rξ) -1 , see (3.10)- (3.11). In the tight-binding situation, we naturally expect the error made in replacing each resolvent (H Rξ) -1 with (H Pξ) -1 (H P is the 'single atom' operator in (1.25)) inside (1.39) to be 'small' for large values of the R-parameter. Arising from this approximation, the 'corrective term' (3.14) of (3.13) involves the trace per unit volume:

|Ω R | -1 Tr L 2 (R 3 ) χ ΩR (H R -ξ) -1 [T R,1 (ξ)T R,1 (ξ) -T R,2 (ξ)] -(H P -ξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] χ ΩR , (1.40)
where T P,j (ξ), j = 1, 2 involve the resolvent (H Pξ) -1 , see (1.27)-(1.28). In view of the operators entering in the trace, to control the R-behavior of (1.40), it is necessary to control for large values of R the error made in approximating χ ΩR (H Rξ) -1 by χ ΩR (H Pξ) -1 (or with χ ΩR from the right) in the Hilbert-Schmidt sense. To that purpose, we use a geometric perturbation theory in Sec. 2. The key-idea consists in isolating in the Wigner-Seitz cell a region close to the boundary:

B R (κ) := {x ∈ Ω R : dist(x, ∂Ω R ) ≤ κR α }
for some κ > 0 and 0 < α < 1, from the bulk where only the 'single atom' operator H P will act. It comes down to approximate (H Rξ) -1 with:

ĝR (H P -ξ) -1 g R + ĝR (H R -ξ) -1 (1 -g R ),
where g R , ĝR and ĝR are smooth cutoff functions; g R , ĝR are supported in the bulk and satisfy (2.2)-(2.3) and ĝR is supported outside and satisfies (2.4)-(2.5). Starting with this approximation, we arrive at (2.13). From the exponential decay of the resolvent's kernel along with the properties (2.3)-(2.5) on the cutoff functions, then provided that the single-site potential u is compactly supported, one gets that the Hilbert-Schmidt norm of:

χ ΩR (H R -ξ) -1 -(H P -ξ) -1 , (H R -ξ) -1 -(H P -ξ) -1 χ ΩR ,
are exponentially small in R α , 0 < α < 1. We emphasize that the assumption u is compactly supported in (A r ) is crucial to get the exponential decay in R α of the norms. Indeed, (2.13) involves the operator defined in (2.10), and the support of χ ΩR ĝR (or (1g R )χ ΩR ) is disjoint from the one of VR only if u is compactly supported. Based on the foregoing, we prove that the trace in (1.40) is exponentially small in R α , 0 < α < 1, see Lemmas 4.2-4.5.

Step 2: Switching to the canonical conditions and performing the zero-temperature limit. Remind that, under the assumptions (A m )-(A nd ), {λ l } τ l=1 , τ ∈ N * denotes the set of eigenvalues of H P in (-∞, 0) counting in increasing order. Now, we switch to the canonical conditions. We assume that the number of particles n 0 ∈ N * in the Wigner-Seitz cell is fixed and obeys n 0 ≤ τ , while the density is given by (1.23). (The reason for this will become clear below). The starting-point is (3.16) corresponding to the sum of (3.13) with (3.14), but with µ (1.24). With the aim of performing the zero-temperature limit in (3.16), we need beforehand to turn to the location for large values of R of the Fermi energy defined as the limit:

(0) R (β, ρ 0 (R), b = 0) instead of µ, see
E R (ρ 0 (R)) := lim β↑∞ µ (0) R (β, ρ 0 (R), b = 0).
Recall that, since H R commutes with the translations of the RZ 3 -lattice, then σ(H R ) is absolutely continuous and consists of the union of compact intervals E R,l (the Bloch bands): ), or in the interior of a Bloch band (this is the metallic situation). In Proposition 3.1, we state that for large values of R, our assumptions (A m )-(A nd ) along with the fact that ρ 0 (R) obeys (1.23), automatically lead to the isolating situation:

σ(H R ) = ∪ ∞ l=1 E R,l , with min E R,l ≤ min E R,l+1 and max E R,l ≤ max E R,
E R (ρ 0 (R)) = max E R,n0 + min E R,n0+1 2 < 0. ( 1
β -1 f β, µ (0)
R (β, ρ 0 (R), b = 0); ξ , f(β, µ; ξ) := ln(1 + e β(µ-ξ) ),
and

∀µ ∈ R, f(β, µ; • ) is holomorphic on {ζ ∈ C : ℑζ ∈ (-π β , π β )}.
The difficulty encountered is as follows. The contours involved in the integration w.r. (we set λ τ +1 := 0) and s.t. I n0 ∩ σ(H R ) = ∅ and I n0 ∩ σ(H P ) = ∅. In this way, we can decompose the contour (3.15) into three parts: the first one enclosing the first n 0 eigenvalues of H P , the second one enclosing I n0 but no eigenvalue, and the third one enclosing the rest of σ(H P ), see (4.2). Note that, for R large enough, these contours do not intersect σ(H R ). Since the second contour lies in the holomorphic domain of the integrand, the integral is null. Since for a given µ 0 ∈ R, f(β, µ 0 ; • ) is holomorphic on {ζ ∈ C : ℜζ = µ 0 }, then the first and third contours can be deformed in some β-independent contours. It remains to use the Lebesgue's dominated theorem knowing (4.5). From (3.16), Propositions 3.6-3.7 together lead to (3.18). The behavior for large R of the remainder in (3.18) follows from the one of (1.40) discussed in Step 1. We mention that, if (A nd ) is disregarded, then the insulating situation can still occur if one sets some restrictions on the n 0 , see Sec. 1.5 and Sec. 1.6.

Step 3: Isolating the main R-dependent contribution -Part 2. We start by removing the Rdependence arising from the indicator functions χ ΩR inside the trace of (3.19). Getting the trace out the integration w.r.t. the ξ-variable, it remains to control the R-behavior of the error made in extending the trace to the whole space. In Proposition 3.9, we state that it is exponentially small in R. To do so, the crucial ingredient is the exponential localization of the eigenfunctions associated with the eigenvalues of H P in (-∞, 0). Next, the Fermi energy can be removed from (3.22) without changing the value of the trace, see Proposition 3.10. Thus, the main R-dependent term behaves like R -3 .

• Outline of the proof of (ii) and (iii).

The proofs heavily rely on the stability of the eigenvalues {λ

l } τ l=1 , τ ∈ N * of H P in (-∞, 0) under the perturbation W (b) := H P (b) -H P .
Outline of the proof of (ii). For b small enough, the key-idea consists in expressing the sum over the n 0 eigenvalues λ l (b) by means of the Riesz integral formula for orthogonal projections, see (3.23). In the rest of the proof, we show that the second derivative w.r.t. b at b = 0 of the r.h.s. of (3.23) is nothing but (1.30), see Proposition 3.11. From the rewriting of the trace in (3.23) as the integral of the diagonal kernel (3.34), it comes down to prove that the integral kernel of (H P (b)-ξ) -1 , far away the diagonal, is twice differentiable in a neighborhood of b = 0, and to write down a formula for its second derivative, see Lemma 3.13. To that purpose, the crucial ingredient is the so-called gauge invariant magnetic perturbation theory, see e.g. [START_REF] Nenciu | On asymptotic perturbation theory for quantum mechanics: Almost invariant subspaces and gauge invariant magnetic perturbation theory[END_REF][START_REF] Briet | Locating the spectrum for magnetic Schrödinger and Dirac operators, Commun[END_REF][START_REF] Cornean | The Faraday effect revisited: Thermodynamic limit[END_REF][START_REF] Briet | A rigorous proof of the Landau-Peierls formula and much more[END_REF][START_REF] Briet | A rigorous approach to the magnetic response in disordered systems[END_REF][START_REF] Savoie | On the zero-field orbital magnetic susceptibility of Bloch electrons in graphene-like solids: Some rigorous results[END_REF]. The idea behind is to isolate the singularity of the magnetic perturbation (arising from the linear growth of the vector potential) via an exponential factor involving the magnetic phase in (3.26). Then (H P (b)ξ) -1 can be approximated by RP (b, 0, ξ) generated by the kernel in (4.21), while the 'corrective term' behaves like O(|b|), see (4.27). Iterating twice (4.27) in the kernels sense, then expanding the exponential factor in Taylor series, one obtains the beginning of an expansion in powers of b for the kernel of (H P (b)ξ) -1 . We mention that the exponential localization of the eigenfunctions associated with the λ j (b)'s for b small enough (see Lemma 3.14) plays a crucial role to control the involved quantities, see Corollary 3.15.

Outline of the proof of (iii). From (1.31), the aim is to derive an expression for the second derivative w.r.t. b at b = 0 of the perturbed eigenvalue λ l (b), l ∈ {1, . . . , n 0 }. To do that, we use the Feshbach formula to write down 'the corrections' to the unperturbed eigenvalue λ l under the perturbation

W (b) = -ba • (-i∇) + b 2
2 a 2 for small values of b, see (3.40). Here, the exponential localization of the eigenfunction associated with λ l allows to control the linear growth of a. Iterating (3.40), we obtain the beginning of an expansion in powers of b for λ l (• ), see (3.41). To conclude, it remains to use that λ l (• ) has an asymptotic expansion series, see e.g. [START_REF] Briet | Locating the spectrum for magnetic Schrödinger and Dirac operators, Commun[END_REF]Thm. 1.2]. Note that the asymptotic perturbation theory and the gauge-invariant magnetic perturbation theory both are used to control the b-behavior of the remainder in (3.41).

Discussion on the assumptions.

Let us first discuss the assumption (A r ) on the single-site potential u. The physical modeling requires that the 'single atom' operator H P possesses finitely many eigenvalues (at least one) below the essential spectrum. Then u has to be chosen accordingly. As emphasized in Sec. 1.4, choosing u compactly supported plays a crucial role in our analysis since it gives rise to the exponentially decreasing behavior in R α , 0 < α < 1 of the remainder in (1.29), see Proposition 3.7. However, we believe that the optimal α is α = 1, i.e. the remainder should behave like O(e -cR ), c > 0. We think that it could be obtained from our analysis by using a more refined geometric perturbation theory to approximate the resolvent (H Rξ) -1 . Furthermore, we point out that the leading term in the asymptotic expansion (1.29) is unchanged if the single-site potential is not compactly supported. Consider the assumption:

(A r * ) u ∈ C 1 (R 3 ; R) with u = O(|x| -(3+ǫ) ) for |x| sufficiently large.
From [44, Thm. XIII.6], σ ess (H P ) = [0, ∞) and H P has a finite number of bound states in (-∞, 0). Replacing (A r ) with (A r * ) in Theorem 1.1, then under the same conditions, one may expect the behavior of the remainder to be only polynomially decreasing with R owing to the 'tail' of the potential. Finally, let us mention that u is chosen continuously differentiable to ensure some regularities for the eigenvectors of H P , see (4.17) and Lemma 3.14.

Let us discuss the non-degeneracy assumption (A nd ). Our analysis is based on the insulating situation which occurs when the Fermi energy lies in the middle of a spectral gap of H R , see [START_REF] Briet | A rigorous proof of the Landau-Peierls formula and much more[END_REF]Thm 1.1]. When the number of particles n 0 in the Wigner-Seitz cell is any integer lesser than the number of negative eigenvalues of H P , while the density is given by (1.23), this together with (A nd ) automatically lead to the insulating situation for R sufficiently large, see Proposition 3.1. Nonetheless, we stress the point that, when disregarding the assumption (A nd ), then the insulating condition can still occur for R sufficiently large provided that one sets some restrictions on the n 0 's in (1.23). It has to obey (see proof of Proposition 3.1 and Remark 3.5):

n 0 ≤ τ and ∃κ ∈ {1, . . . , ν} s.t. n 0 = dim E 1 + • • • + dim E κ . (1.42)
Here, τ is the number of eigenvalues of H P in (-∞, 0) counting multiplicities, ν ≤ τ is the number of distinct eigenvalues and E l , l ∈ {1, . . . , ν} stands for the eigenspace associated with the (possibly degenerate) eigenvalue λ l of H P , λ l = {λ (m) l } dimE l m=1 . Supposing that the assumptions (A r )-(A m ) hold, and assuming that the number of particles n 0 ∈ N * in the Wigner-Seitz cell is fixed and obeys (1.42), while the density of particles is given by (1.23), then the formulae in Theorem 1.1 (ii)-(iii) have to be modified accordingly (the statement in (i) is unchanged). Indeed, (1.31) becomes:

X P (n 0 ) = - q c 2 κ l=1 dimE l m=1 d 2 λ (m) l db 2 (b = 0);
and (1.33)-(1.34) respectively become (with intuitive notations):

1 |Ω R | X La (n 0 ) := - q c 2 1 4|Ω R | κ l=1 dimE l m=1 Φ (m) l , X 2 1 + X 2 2 Φ (m) l L 2 (R 3 ) , 1 |Ω R | X vV (n 0 ) := q c 2 1 2|Ω R | κ l=1 dimE l m=1 L 3 Φ (m) l , Π (m),⊥ l H P -λ (m) l Π (m),⊥ l -1 L 3 Φ (m) l L 2 (R 3 )
.

Note that the contribution X vV (n 0 ) can be decompose into two contributions as in Remark 1.3.

An open problem.

From the foregoing, a natural question arises: does the insulating situation occur when n 0 is any integer less than the number of negative eigenvalues of H P counting multiplicities if some of the eigenvalues are degenerate, while the density is given by (1.23)? Such a problem comes up when the single-site potential is chosen spherically symmetric. To tackle it, we need to know precisely the behavior of the negative spectral bands of H R near the degenerate eigenvalues of H P in (-∞, 0) for large values of the R-parameter. For instance, suppose that one of the negative eigenvalue of H P , say λ c , is two-fold degenerate. For R sufficiently large, it is well-known that the spectrum of H R in a neighborhood of λ c consists of the union of two Bloch bands, see Lemma 3.4 and also [START_REF] Daumer | Equation de Schrödinger avec champ électrique périodique et champ magnétique constant dans l'approximation du tight-binding[END_REF]Thm. 2.1]. But we need to know much more; in particular we need to control how the Bloch bands behave the one relative to the other. Especially, do they always overlap for R sufficiently large or do they overlap only in the limit R ↑ ∞? How fast each of the Bloch bands reduces to λ c ? This remains a challenging spectral problem.

The content of the paper.

Our current paper is organized as follows. In Sec. 2, we use a geometric perturbation theory to approximate the resolvent operator (H Rξ) -1 in the tight-binding situation. Sec. 3 is devoted to the proof of Theorem 1.1. In Sec. 3.1.1, we show that under our assumptions, only the insulating situation can occur in the tight-binding situation. Subsequently, in Sec. 3.1.2-3.1.3 we prove the asymptotic expansion in (1.29). In Sec. 3.2, we prove the identity (1.31). In Sec. 3.3 we prove the formula (1.32). In Sec. 4, we have collected all the proofs of the technical intermediary results needed in Sec. 3.

2 An approximation of the resolvent via a geometric perturbation theory.

The method we use below is borrowed from [START_REF] Cornean | The Faraday effect revisited: Thermodynamic limit[END_REF][START_REF] Cornean | Sharp trace asymptotics for a class of 2D-magnetic operators[END_REF]. For any 0 < α < 1, 0 < κ ≤ 10 and R > 0 define:

B R (κ) := x ∈ Ω R : dist (x, ∂Ω R ) ≤ κR α .
Hence, for R sufficiently large, B R (κ) models a 'thin' compact subset of Ω R near the boundary with Lebesgue-measure |B R (κ)| of order O(R 2+α ). Let 0 < α < 1 be fixed. Below, by R sufficiently large, we mean:

R ≥ R 0 with R 0 = R 0 (α) ≥ 1 s.t. B R0 (8) Ω R0 . (2.1)
Let us introduce some well-chosen families of smooth cutoff functions. Let {g R } and {ĝ R } satisfying for R ≥ R 0 :

Supp(g R ) ⊂ (Ω R \ B R (3)) , g R = 1 if x ∈ (Ω R \ B R (4)) , 0 ≤ g R ≤ 1; Supp(ĝ R ) ⊂ (Ω R \ B R (1)) , ĝR = 1 if x ∈ (Ω R \ B R (2)) , 0 ≤ ĝR ≤ 1.
Moreover, there exists a constant C > 0 s.t.

∀R ≥ R 0 , max { D s g R ∞ , D s ĝR ∞ } ≤ CR -|s|α , ∀|s| ≤ 2, s ∈ N 3 . Also, let { ĝR } satisfying for R ≥ R 0 : Supp ĝR ⊂ Ω R \ B R (6) ∁ , ĝR = 1 if x ∈ Ω R \ B R (5) ∁ , 0 ≤ ĝR ≤ 1.
Moreover, there exists a constant C > 0 s.t.

∀R ≥ R 0 , D s ĝR ∞ ≤ CR -|s|α , ∀|s| ≤ 2, s ∈ N 3 .
By gathering all these properties together, one straightforwardly gets:

ĝR g R = g R , (2.2)
dist (Supp(D s ĝR ), Supp(g R ))

≥ CR α , ∀1 ≤ |s| ≤ 2, (2.3) ĝR (1 -g R ) = (1 -g R ), (2.4) dist Supp D s ĝR , Supp (1 -g R ) ≥ CR α , ∀1 ≤ |s| ≤ 2, (2.5)
for some R-independent constant C > 0.

Let us now define a series of operators. At first, introduce ∀R ≥ R 0 and ∀ξ ∈ ̺(H R ) ∩ ̺(H P ) (here ̺(• ) denotes the resolvent set) on L 2 (R 3 ):

R R (ξ) := ĝR (H P -ξ) -1 g R + ĝR (H R -ξ) -1 (1 -g R ).
(2.6)

Due to the support of the cutoff functions, one has:

(H R -ξ)ĝ R = (H P -ξ)ĝ R .
Using that Ran(R R (ξ)) ⊂ Dom(H R ), then from (2.2) along with (2.4):

(H R -ξ)R R (ξ) = ½ + W R (ξ),
where, ∀R ≥ R 0 and ∀ξ ∈ ̺(H R ) ∩ ̺(H P ):

W R (ξ) := - 1 2 (∆ĝ R ) -(∇ĝ R ) • ∇ (H P -ξ) -1 g R + - 1 2 ∆ ĝR -∇ ĝR • ∇ (H R -ξ) -1 (1 -g R ) .
Since W R (ξ) is bounded on L 2 (R 3 ), see e.g. [START_REF] Briet | A rigorous approach to the magnetic response in disordered systems[END_REF]Lem. 5.1], this means that:

(H R -ξ) -1 = R R (ξ) -(H R -ξ) -1 W R (ξ). (2.7)
Next, R R (ξ) in (2.6) can be rewritten by the second resolvent equation as:

R R (ξ) = R R (ξ) -W R (ξ), (2.8)
where, ∀R ≥ R 0 and ∀ξ ∈ ̺(H R ) ∩ ̺(H P ):

R R (ξ) := ĝR (H P -ξ) -1 g R + ĝR (H P -ξ) -1 (1 -g R ), (2.9)
W R (ξ) := ĝR (H R -ξ) -1 VR (H P -ξ) -1 (1 -g R ), (2.10)
with: VR :=

υ∈Z 3 \{0} u(• -Rυ). (2.11)
Finally, R R (ξ) in (2.9) can be rewritten as:

R R (ξ) = (H P -ξ) -1 + W R (ξ), (2.12)
where, ∀R ≥ R 0 and ∀ξ ∈ ̺(H R ) ∩ ̺(H P ):

W R (ξ) :=(H P -ξ) -1 - 1 2 (∆ĝ R ) -(∇ĝ R ) • ∇ (H P -ξ) -1 g R + (H P -ξ) -1 - 1 2 ∆ ĝR -∇ ĝR • ∇ (H P -ξ) -1 (1 -g R ) .
Gathering (2.7), (2.8) and (2.12) together, ∀R ≥ R 0 and ∀ξ ∈ ̺(H R ) ∩ ̺(H P ):

(H R -ξ) -1 = (H P -ξ) -1 + W R (ξ) -W R (ξ) -(H R -ξ) -1 W R (ξ), (2.13)
and (2.13) holds in the bounded operators sense on L 2 (R 3 ). From (2.13) and taking into account the assumption (A r ), let us now prove that for R ≥ R 0 sufficiently large, the error made in approximating χ ΩR (H Rξ) -1 with χ ΩR (H Pξ) -1 (or with the indicator function χ ΩR from the right) is exponentially small in R α , 0 < α < 1 in the Hilbert-Schmidt norm sense. To do that, define R 1 ≥ 1 so that:

Supp(u) ⊂ (Ω R1 \ B R1 (7)).
(2.14)

Such R 1 exists since the support of u is compact, see assumption (A r ). Now, look at the r.h.s. of (2.13). For any R ≥ R 0 , (H Rξ) -1 W R (ξ) and W R (ξ) have their operator norm exponentially small in R α , see (2.25)-(2.24) below. However, this is not the case for the operator norm of W R (ξ) even for R ≥ R 0 and large enough, see (2.24). This comes from the fact that the support of ĝR (or (1g R )) and the one of VR in (2.11) are not disjoint. But for any R ≥ max{R 0 , R 1 }, the Hilbert-Schmidt (H-S) norm of W R (ξ) when multiplied by the indicator function χ ΩR from the left (or from the right) is exponentially small in R α , see (2.27). The same holds true for the H-S norms of (H Rξ) -1 W R (ξ) and W R (ξ) when multiplied by the indicator function χ ΩR , see (2.26). This last feature, which results from the fact that u is compactly supported, will turn out to be decisive to get the exponential decay of the remainder term in the asymptotic expansion (1.29).

We end this section by giving a series of estimates we will use throughout. For any ξ ∈ C and real number ℓ > 0, we use the shorthand notation:

ℓ ξ := ℓ(1 + |ξ|) -1 .
(2.15)

Lemma 2.1 Let Ξ = P or R. For every η > 0 (and ∀R > 0 if Ξ = R), there exists a constant ϑ = ϑ(η) > 0 and a polynomial p(• ) s.t. ∀ξ ∈ C satisfying dist(ξ, σ(H Ξ )) ≥ η and ∀(x, y) ∈ R 6 \D:

(H Ξ -ξ) -1 (x, y) ≤ p(|ξ|) e -ϑ ξ |x-y| |x -y| , (2.16
) Lemma 2.2 Let 0 < α < 1 and R 0 = R 0 (α) ≥ 1 as in (2.1). For every η > 0, there exists a constant ϑ = ϑ(η) > 0 and a polynomial p(• ) s.t.

∇ x (H Ξ -ξ) -1 (x, y) ≤ p(|ξ|) e -ϑ ξ |x-y| |x -y| 2 . ( 2
(i). ∀R ≥ R 0 , ∀ξ ∈ C obeying dist(ξ, σ(H R ) ∩ σ(H P )) ≥ η, ∀(x, y) ∈ R 6 \ D: max {|(R R (ξ))(x, y)| , |(R R (ξ))(x, y)|} ≤ p(|ξ|) e -ϑ ξ |x-y| |x -y| , (2.18) max {|∇ x (R R (ξ))(x, y)| , |∇ x (R R (ξ))(x, y)|} ≤ p(|ξ|) e -ϑ ξ |x-y| |x -y| 2 . (2.19) (ii). ∀R ≥ R 0 , ∀ξ ∈ C obeying dist(ξ, σ(H R ) ∩ σ(H P )) ≥ η and ∀(x, y) ∈ R 6 : |(W R (ξ))(x, y)| ≤ p(|ξ|)e -ϑ ξ |x-y| , (2.20
)

|∇ x (W R (ξ))(x, y)| ≤ p(|ξ|) e -ϑ ξ |x-y| |x -y| , x = y, (2.21) max {|(W R (ξ))(x, y)| , |(W R (ξ))(x, y)| , |∇ x (W R (ξ))(x, y)|} ≤ p(|ξ|)e -ϑ ξ R α e -ϑ ξ |x-y| . (2.22) (iii). ∀R ≥ max{R 0 , R 1 } (R 1 is defined through (2.14)), ∀ξ ∈ C obeying dist(ξ, σ(H R ) ∩ σ(H P )) ≥ η and ∀(x, y) ∈ R 6 : max χ ΩR ĝR (x)(H R -ξ) -1 (x, y) VR (y) , VR (x)(H P -ξ) -1 (x, y) (χ ΩR (1 -g R )) (y) ≤ p(|ξ|)e -ϑ ξ R α e -ϑ ξ |x-y| .
Proof. We start by (i). In view of (2.6) and (2.9), (2.18) directly follows from (2.16) knowing that 0 ≤ g R , ĝR , ĝR ≤ 1. Next, one has ∀(x, y) ∈ R 6 \ D:

∇ x (R R (ξ))(x, y) = (∇ĝ R ) (x)(H P -ξ) -1 (x, y) + ĝR (x)∇ x (H P -ξ) -1 (x, y) g R (y) + ∇ ĝR (x)(H R -ξ) -1 (x, y) + ĝR (x)∇ x (H R -ξ) -1 (x, y) (1 -g R (y)),
and by replacing (H Rξ) -1 with (H Pξ) -1 above, we get

∇ x (R R (ξ))(• , • ).
∇ x (W R (ξ))(x, y) = ∇ ĝR (x) R 3 dz (H R -ξ) -1 (x, z) VR (z)(H P -ξ) -1 (z, y)(1 -g R )(y) + ĝR (x) R 3 dz ∇ x (H R -ξ) -1 (x, z) VR (z)(H P -ξ) -1 (z, y)(1 -g R )(y).
Then

max ∇ ĝR (x)∇ x (H Ξ -ξ) -1 (x, y)(1 -g R )(y) , (∇ĝ R ) (x)∇ x (H Ξ -ξ) -1 (x, y)g R (y) , (∆ĝ R ) (x)(H Ξ -ξ) -1 (x, y)g R (y) , ∆ ĝR (x)(H Ξ -ξ) -1 (x, y)(1 -g R )(y) ≤ p(|ξ|)e -ϑ ξ R α e -ϑ ξ |x-y| , (2.23)
for another R-independent ϑ > 0 and polynomial p(• ). This leads to (2.22) for the kernel of W R (ξ) and W R (ξ). As for

∇ x (W R (ξ))(• , • ),
max (H R -ξ) -1 , R R (ξ) , R R (ξ) , W R (ξ) , ∇R R (ξ) , ∇R R (ξ) , ∇W R (ξ) ≤ p(|ξ|), (2.24) max { W R (ξ) , W R (ξ) ∇W R (ξ) } ≤ p(|ξ|)e -ϑ ξ R α , (2.25)
for another R-independent constant ϑ > 0 and polynomial p(• ).

Remark 2.4 Let (I 2 (L 2 (R 3 )), • I2 ) be the Banach space of Hilbert-Schmidt operators. From Lemma 2.2 (ii) and the * -ideal property of I 2 (L 2 (R 3 )):

∀R ≥ R 0 , max χ ΩR W R (ξ) I2 , χ ΩR (H R -ξ) -1 W R (ξ) I2 , W R (ξ)χ ΩR I2 , (H R -ξ) -1 W R (ξ)χ ΩR I2 ≤ p(|ξ|)e -ϑ ξ R α , (2.26)
and from Lemma 2.2 (iii) one has ∀R ≥ max{R 0 , R 1 }:

max { χ ΩR W R (ξ) I2 , W R (ξ)χ ΩR I2 } ≤ p(|ξ|)e -ϑ ξ R α , (2.27)
for another R-independent constant ϑ > 0 and polynomial p(• ).

3 Proof of Theorem 1.1.

This section is organized as follows. The first part is devoted to the proof of the asymptotic expansion (1.29) in the tight-binding situation and in the zero-temperature limit. The second and third part are respectively concerned with the proof of (1.31) and (1.32). For reader's convenience, the proofs of technical intermediary results are collected in Appendix, see Sec. 4.

Proof of (i).

The location of the Fermi energy.

Here, we are interested in the location of the Fermi energy in the tight-binding situation when the number of particles n 0 ∈ N * in the Wigner-Seitz cell is fixed, while the density is given by (1.23). We recall that under our conditions, the Fermi energy defined as:

E R (ρ 0 (R)) := lim β↑∞ µ (0) R (β, ρ 0 (R), b = 0), R > 0, (3.1)
always exists, see [START_REF] Briet | A rigorous proof of the Landau-Peierls formula and much more[END_REF]Thm 1.1].

Before stating the main result, let us introduce some basic elements of Bloch-Floquet theory. We refer to [6, Sec. 3.5] and also [START_REF] Wilcox | Theory of Bloch waves[END_REF]. We mention that the results we give below hold true for any R > 0. Denote by Ω * R the unit cell of the dual lattice (2π/R)Z 3 (the so-called first Brillouin zone) of the Bravais lattice RZ 3 . Denoting by S (R 3 ) the Schwartz space of rapidly decreasing functions on R 3 , the Bloch-Floquet(-Zak) transformation is defined by:

U : S (R 3 ) → L 2 Ω * R , L 2 (Ω R ) ∼ = ⊕ Ω * R dk L 2 (Ω R ) (Uφ)(x; k) = 1 |Ω * R | υ∈RZ 3 e -ik•(x+υ) φ(x + υ), k ∈ Ω * R , x ∈ Ω R ,
which can be continued in a unitary operator on L 2 (R 3 ). The unitary transformation of H R is decomposable into a direct integral:

UH R U * = ⊕ Ω * R dk h R (k), h R (k) := 1 2 (-i∇ + k) 2 + V R ,
and the Bloch Hamiltonian h R lives in L 2 (R 3 /RZ 3 ). By standard arguments, h R is essentially self-adjoint on C ∞ (R 3 /RZ 3 ) and the domain of its closure is H 2 (R 3 /RZ 3 ). For each k ∈ Ω * R , h R (k) has purely discrete spectrum with an accumulation point at infinity. We denote by {E R,l (k)} l≥1 the set of eigenvalues counting multiplicities and in increasing order. Due to this choice of labeling, the E R,l 's are periodic and Lipschitz continuous on Ω * R . Indeed, when crossing-points occur, the E R,l 's are not differentiable on a zero Lebesgue-measure subset of Ω * R corresponding to crossing-points. The spectrum of H R is absolutely continuous and given (as a set of points) by σ(

H R ) = ∞ l=1 E R,l where ∀l ≥ 1, E R,l := [min k∈Ω * R E R,l (k), max k∈Ω * R E R,l (k)
] stands for the l-th Bloch band function. Note that the sets E R,l can overlap each other in many ways, and some of them can coincide. The energy bands are disjoint unions of E R,l 's. If max E R,l < min E R,l+1 for some l ≥ 1, then a spectral gap occurs. Since the Bethe-Sommerfeld conjecture holds true under our conditions, see e.g. [28, Coro. 2.3], then the number of spectral gaps is finite, if not zero.

Our main result below states that, under our conditions and in the tight-binding situation, the Fermi energy always lies in the middle of a spectral gap of H R (in other words, only the insulating situation can occur, see [11, Thm. 1.1]), and moreover, it provides an asymptotic expansion: Proposition 3.1 Suppose that assumptions (A r ), (A m ) and (A nd ) hold, and assume that the number of particles n 0 ∈ N * in the Wigner-Seitz cell is fixed and satisfies n 0 ≤ τ , while the density is given by (1.23). (i). For any β > 0, let µ

(0) R (β, ρ 0 (R), b = 0) ∈ R be the unique solution of the equation ρ R (β, µ, b = 0) = ρ 0 (R). Then for R sufficiently large, E R (ρ 0 (R)) := lim β↑∞ µ (0) R (β, ρ 0 (R), b = 0) = max E R,n0 + min E R,n0+1 2 < 0. (3.2)
(ii). Define:

E P (n 0 ) := 1 2 × (λ n0 + λ n0+1 ) if n 0 < τ , λ τ if n 0 = τ . (3.3)
Then, under the additional assumption that n 0 < τ , one has:

E R (ρ 0 (R)) = E P (n 0 ) + O e - √ |λn 0 +1|R .

Remark 3.2 (i).

We stress the point that the non-degeneracy assumption (A nd ) along with the fact that ρ 0 (R) is given as in (1.23), together imply the insulating situation for R sufficiently large.

(ii). We will see in Proposition 3.10 that the Fermi energy plays in fact no role in the statements of Theorem 1.1 since it can be removed from the quantities of interest without changing their values.

The rest of this paragraph is devoted to the proof of Proposition 3.1. Let us start by writing down an expression for the bulk density of particles. Under the grandcanonical conditions, let β > 0 and µ ∈ R. ∀R > 0, let C (R) β be the positively oriented simple contour around [inf σ(H R ), ∞) defined as:

C (R) β := ℜξ ∈ [δ R , ∞), ℑξ = ± π 2β ∪ ℜξ = δ R , ℑξ ∈ - π 2β , π 2β , (3.4)
δ R := inf σ(H R ) -1.
Let us mention that, for any R > 0, the closed subset surrounding by

C (R) β
is a strict subset of the holomorphic domain D := {ζ ∈ C : ℑζ ∈ (-π β , π β )} of the Fermi-Dirac distribution function f F D (β, µ; ξ) := e β(µ-ξ) (1 + e β(µ-ξ) ) -1 . From (1.21), and seen as a function of the µ-variable, the bulk zero-field density of particles reads ∀β > 0, ∀µ ∈ R, ∀R > 0 as, see e.g. [12, Eq. (6.3)]:

ρ R (β, µ, b = 0) := 1 |Ω R | i 2π Tr L 2 (R 3 ) χ ΩR C (R) β dξ f F D (β, µ; ξ)(H R -ξ) -1 χ ΩR . (3.5)
Another way to express the bulk zero-field density consists in bringing into play the integrated density of states of H R . Under the conditions of (3.5),

ρ R (β, µ, b = 0) = - ∞ -∞ dt ∂f F D ∂t (β, µ; t)N R (t),
where N R is the integrated density of states of the operator H R = H R (b = 0) defined in (1.18). We recall that, when the magnetic field vanishes, N R is a positive, continuous and non-decreasing function, and N R is piecewise constant when the energy parameter belongs to a spectral gap. Next, in order to write down a convenient expression for the bulk zero-field density of particles in the zero-temperature limit, we need to rewrite (3.5) by the use of the Bloch-Floquet decomposition, see [START_REF] Briet | A rigorous proof of the Landau-Peierls formula and much more[END_REF]Sec. 2]. Using the notation introduced previously, all the needed results are collected in:

Lemma 3.3 (i). Let β > 0 and µ ∈ R. Then for any R > 0: ρ R (β, µ, b = 0) = 1 |Ω R ||Ω * R | ∞ j=1 Ω * R dk f F D (β, µ; E R,j (k)).
(ii). For any R > 0, let µ ≥ inf σ(H R ) be fixed. One has the identity:

lim β↑∞ ρ R (β, µ, b = 0) = 1 |Ω R ||Ω * R | ∞ j=1 Ω * R dk χ [inf σ(HR),µ] (E R,j (k)) = N R (µ), (3.6)
where χ [inf σ(HR),µ] denotes the indicator function of the compact interval [inf σ(H R ), µ], and N R the integrated density of states of the operator H R .

Now, let us get back to the location of the Fermi energy in the tight-binding situation when the density is given by (1.23). We need first to know how the negative spectral bands of the operator H R are localized at negative eigenvalues of the operator H P for large values of the R-parameter. For completeness' sake, in the below lemma, we disregard assumption (A nd ) and we allow the negative eigenvalues of H P to be (possibly) finitely degenerate: Lemma 3.4 Let τ ∈ N * be the number of the eigenvalues of H P in (-∞, 0) counting multiplicities. Denote by {λ l } ν l=1 , ν ≤ τ the set of distinct eigenvalues counting in increasing order. Then, for R sufficiently large, there exist 2ν + 1 real numbers: c R,l , d R,l < 0 with l = 1, . . . , ν and C R,ν+1 satisfying:

-λ 1 - 1 2 < c R,1 < d R,1 < • • • < c R,ν < d R,ν < C R,ν+1 , s.t. (i) σ(H R ) ↾(-∞,0) ⊂ ν l=1 [c R,l , d R,l ], (ii) [c R,l , d R,l ] ∩ [c R,m , d R,m ] = ∅ if l = m, (iii) λ l ∈ (c R,l , d R,l ), (iv) d R,ν + C R,ν+1 < 0;
together with the properties that c R,l , d R,l → λ l and C R,ν+1 → 0 when R ↑ ∞. Furthermore, one has the following relations with the Bloch bands of H R :

[c R,l , d R,l ] ∩ σ(H R ) = dimE l m=1 E R,m , l = 1, . . . , ν, C R,ν+1 = min E R,τ +1 .
Here, E l denotes the eigenspace associated with the degenerate eigenvalue λ l .

Proof. The lemma follows from [21, Thm. 2.1] taking into account the labeling in increasing order for the E R,l 's. Note that [21, Thm. 2.1] is stated under the assumptions that V R is smooth and sufficiently fast decaying at infinity. But the statements still hold true under our conditions on the potential V R , see [START_REF] Galloonov | Estimations for negative spectral bands of threedimensional periodical Schrödinger operator[END_REF]Thm. 2]. Now we are ready for: Proof of Proposition 3.1. Let us first prove (i). Consider the equation:

1 |Ω * R | l≥1 Ω * R dk χ [inf σ(HR),E] (E R,l (k)) = n 0 .
Due to the assumption (A nd ), Lemma 3.4 ensures that the Bloch bands E R,l , l = 1, . . . , τ are simple, isolated from each other and from the rest of the spectrum for R sufficiently large. Hence,

if n 0 ∈ N * satisfies n 0 ≤ τ , then E must belong to [max k∈Ω * R E R,n0 (k), min k∈Ω * R E R,n0+1 (k)
]. This comes from the fact that the Lebesgue-measure of the set {k ∈ Ω

* R : E R,l (k) ≤ E} equals |Ω * R | if and only if E ≥ max k∈Ω * R E R,l (
k). Getting back to (3.6), this means that for R sufficiently large: Let us turn to (ii). Assume that n 0 < τ . It is enough to use that:

∀E ∈ [max E R,n0 , min E R,n0+1 ], N R (E) = ρ 0 (R). ( 3
E R (ρ 0 (R)) -E P (n 0 ) = max E R,n0 -λ n0 2 + min E R,n0+1 -λ n0+1 2 ,
along with the estimate which holds uniformly in k, see e. We start by writing down an expression for the bulk zero-field orbital susceptibility. In the grand-canonical situation, let β > 0 and µ ∈ R. For any R > 0, let C (R) β be the positively oriented simple contour around [inf σ(H R ), ∞) defined in (3.4). The closed subset surrounding by

C (R) β is a strict subset of the holomorphic domain D := {ζ ∈ C : ℑζ ∈ (-π β , π β )} of ξ → f(β, µ; ξ) := ln(1 + e β(µ-ξ) ) satisfying -1 β (∂ ξ f)(β, µ; ξ) = f F D (β, µ; ξ). Note that f(β, µ; • ) admits an exponential decay on C (R)
β , i.e. there exists a constant c > 0 s.t.

∀β > 0, ∀ξ ∈ C (R) β , |f(β, µ; ξ)| ≤ ce βµ e -βℜξ .
(3.8)

From (1.22), and seen as a function of µ, the bulk zero-field orbital susceptibility reads ∀β > 0, ∀µ ∈ R and ∀R > 0 as, see e.g. [START_REF] Briet | A rigorous approach to the magnetic response in disordered systems[END_REF]Eq. (1.21)]:

X R (β, µ, b = 0) := q c 2 2 β|Ω R | i 2π C (R) β dξ f(β, µ; ξ) × Tr L 2 (R 3 ) χ ΩR (H R -ξ) -1 [T R,1 (ξ)T R,1 (ξ) -T R,2 (ξ)] χ ΩR , (3.9)
where T R,j (ξ), j = 1, 2 are bounded operators on L 2 (R 3 ) generated via their kernel respectively defined on R 6 \ D (D := {(x, y) ∈ R 6 : x = y}) as:

T R,1 (x, y; ξ) := a(x -y) • (i∇ x )(H R -ξ) -1 (x, y), (3.10) T R,2 (x, y; ξ) := 1 2 a 2 (x -y)(H R -ξ) -1 (x, y). (3.11)
Remind that we have introduced the operators T P,j (ξ), j = 1, 2 via their kernel defined similarly to (3.10)-(3.11) but with (H P -ξ) -1 instead of (H R -ξ) -1 , see (1.27)-(1.28). Since |a(x-y)| ≤ |x-y|, then under the conditions of Lemma 2.1 (below Ξ := R or P ), one has on R 6 \ D:

|T Ξ,j (x, y; ξ)| ≤ p(|ξ|) e -ϑ ξ |x-y| |x -y| , ϑ ξ := ϑ 1 + |ξ| , j = 1, 2, (3.12)
for another constant ϑ > 0 and polynomial p(• ). Due to (3.12) and (2.16), the operators

(H Ξ -ξ) -1 T Ξ,1 (ξ)T Ξ,1 (ξ) and (H Ξ -ξ) -1 T Ξ,2 (
ξ), Ξ = R or P are locally trace class on L 2 (R 3 ). Furthermore, both are integral operators with a jointly continuous integral kernel on R 6 , whose diagonal part is bounded above by some polynomial in |ξ| uniformly in the spatial variable, see e.g. [START_REF] Briet | A rigorous approach to the magnetic response in disordered systems[END_REF]Lem. A.1]. This along with (3.8) ensure that (3.9) is well-defined.

Now, let us turn to the proof of (1.29). We point out that the main difficulty consists in isolating the main R-dependent contribution from (3.9) in the tight-binding situation, while keeping a good control on the behavior of the 'remainder' term, even in the zero-temperature limit.

The starting point is the approximation of the resolvent (H Rξ) -1 derived in Sec. 2. By replacing in (3.9) each resolvent (H Rξ) -1 (including the ones appearing in (3.10)-(3.11)) with the r.h.s. of (2.13), and taking into account the features of the three last operators in the r.h.s. of (2.13) we discussed in Sec. 2, we naturally expect the main R-dependent contribution from (3.9) to be obtained by replacing each (H Rξ) -1 with (H Pξ) -1 . In this way, define ∀β > 0, ∀µ ∈ R and ∀R > 0, the following quantities:

X R (β, µ, b = 0) := q c 2 1 β 1 |Ω R | i π C (P ) β dξ f(β, µ; ξ) × Tr L 2 (R 3 ) χ ΩR (H P -ξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] χ ΩR , (3.13) ∆ R (β, µ) := X R (β, µ, b = 0) -X R (β, µ, b = 0), (3.14)
where

C (P ) β
appearing in (3.13) denotes the counter-clockwise oriented simple contour around the interval [inf σ(H P ), ∞) defined by:

C (P ) β := ℜξ ∈ [δ P , ∞), ℑξ = ± π 2β ∪ ℜξ = δ P , ℑξ ∈ - π 2β , π 2β
(3.15)

δ P := inf σ(H P ) -1 = λ 1 -1.
Now, we switch to the canonical conditions. Assume that the number of particles n 0 ∈ N * in the Wigner-Seitz cell is fixed and obeys n 0 ≤ τ , while the density is given by (1.23). Let µ

X R (β, ρ 0 (R), b = 0) = X R β, µ (0) R (β, ρ 0 (R), b = 0), b = 0 + ∆ R β, µ (0) R (β, ρ 0 (R), b = 0) .
(3.16) The next step of the proof consists in performing the zero-temperature limit in (3.16) in the tightbinding situation. Here, the crucial point is the insulating situation: for R sufficiently large, the Fermi energy lies outside the spectrum of H R and is located in a neighborhood of the middle of the interval (λ n0 , λ n0+1 ) if n 0 < τ , (λ τ , 0) otherwise; see Proposition 3.1 along with Lemma 3.4. Remind that the insulating situation results from the non-degeneracy assumption (A nd ) together with the condition (1.23). To perform the zero-temperature limit in (3.16), we need the two following propositions. Their proof are placed in Appendix, see Sec. 4.1. Recall that {λ l } τ l=1 , τ ∈ N * denotes the set of eigenvalues of H P in (-∞, 0) counting in increasing order.

lim β↑∞ X R (β, µ, b = 0) = 1 |Ω R | X R (µ, b = 0),
uniformly in µ ∈ K, with:

X R (µ, b = 0) := q c 2 i π Γς dξ (µ -ξ) × Tr L 2 (R 3 ) χ ΩR (H P -ξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] χ ΩR .
Γ ς is any positively oriented simple closed contour surrounding the ς smallest eigenvalues of H P in (-∞, 0) while letting outside the rest of the spectrum. Let us emphasize that the exponentially decaying estimate appearing in (3.17) arises from the fact that u is compactly supported, see assumption (A r ). The proof of Proposition 3.7 is essentially based on the features of the three last operators in the r.h.s. of (2.13) we mentioned in Sec. 2.

Subsequently to Propositions 3.6 and 3.7, we are in a position to isolate a first main Rdependent contribution from (3.16) in the tight-binding situation and in the zero-temperature regime. Under the conditions of (3.16), let us show that ∀0 < α < 1, there exists a R-independent

constant c > 0 s.t. lim β↑∞ X R (β, ρ 0 (R), b = 0) = 1 |Ω R | X R (E R (ρ 0 (R)), b = 0) + O e -cR α , (3.18)
with:

X R (E R (ρ 0 (R)), b = 0) := q c 2 i π Γn 0 dξ (E R (ρ 0 (R)) -ξ) ×Tr L 2 (R 3 ) χ ΩR (H P -ξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] χ ΩR . (3.19)
To achieve that, let I n0 be an open interval s.t. I n0 (λ n0 , λ n0+1 ) and Remark 3.8 In the case of n 0 < τ , due to the asymptotic expansion in Lemma 3.1 (ii) along with (3.17), then one obtains instead of (3.19):

λn 0 +λn 0 +1 2 ∈ I n0 if n 0 < τ ; otherwise I n0 (λ τ ,
lim β↑∞ X R (β, ρ 0 (R), b = 0) = 1 |Ω R | X R (E P (n 0 ), b = 0) + O e -cR α , with X R (E P (n 0 ), b = 0) as in (3.19) but with E P (n 0 ) instead of E R (ρ 0 (R)).
3.1.3 Isolating the main R-dependent contribution at zero-temperature -Continuation and end.

The continuation of the proof of (1.29) consists in removing the R-dependence arising from the presence of the indicator functions χ ΩR inside the trace in (3.19).

The key-ingredient is Proposition 3.9 below whose proof lies in Appendix, see Sec. 4.2. The exponential localization of the eigenfunctions associated with the eigenvalues of H P in (-∞, 0) plays a crucial role in the proof.

Introduce the family of functions g θ,w : C → C, θ ∈ C, w = 0, 1, defined as:

g θ,1 (ξ) := θ -ξ, g θ,0 (ξ) := θ, (ξ, θ) ∈ C 2 .
(3.20) Proposition 3.9 ∀θ ∈ C, ∀w ∈ {0, 1} and ∀ς ∈ {1, . . . , τ }, there exist two constants C = C(θ, ς) > 0 and c > 0 s.t. ∀x ∈ R 3 and for any j ∈ {1, 2}:

max i 2π Γς dξ g θ,w (ξ) (H P -ξ) -1 T P,1 (ξ)T P,1 (ξ) (x, x) , i 2π Γς dξ g θ,w (ξ) (H P -ξ) -1 T P,j (ξ) (x, x) ≤ Ce -c|x| . (3.21)
Γ ς is any positively oriented simple closed contour surrounding the ς smallest eigenvalues of H P in (-∞, 0), while letting outside the rest of the spectrum.

As a result of Proposition 3.9, the error made in getting the trace out the integration w.r.t. ξ in (3.19) and extending the trace to the whole space behaves like O(e -cR ) for some R-independent constant c > 0. Hence, under the conditions of (3.19), there exists a R-independent c > 0 s.t.:

X R (E R (ρ 0 (R)), b = 0) = q c 2 i π Tr L 2 (R 3 ) Γn 0 dξ (E R (ρ 0 (R)) -ξ)
× (H Pξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] + O e -cR .

(3.22)

Gathering (3.18), (3.19) and (3.22) together, then to complete the proof of (1.29)-(1.30), it remains to show that the quantity containing the Fermi energy inside the trace of the leading term in the r.h.s. of (3.22) plays no role (i.e. we can get rid of it without changing the value of the trace). This is contained in the below result, whose proof lies in Appendix, see Sec. 4.2: Proposition 3.10 With the notation of Proposition 3.9, one has:

i 2π Tr L 2 (R 3 ) Γς dξ (H P -ξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] = 0.

Proof of (ii).

Before beginning, let us recall and introduce some notation. Under the assumptions (A r )-(A m )-(A nd ), let {λ l } τ l=1 , τ ∈ N * be the set of eigenvalues of H P in (-∞, 0) counting in increasing order. Let γ l , l = 1, . . . , τ be positively oriented simple closed contours assumed to be two by two disjoint, chosen in such a way that ∀l ∈ {1, . . . , τ }, γ l surrounds the eigenvalue λ l , while letting outside the rest of the spectrum of H P . Let H P (b), b ∈ R be the magnetic 'single atom' operator in (1.26). From [7, Thm. 1.1] (see also [START_REF] Nenciu | Stability of energy gaps under variations of the magnetic field[END_REF][START_REF] Helffer | Equation de Schrödinger avec champ magnétique et équation de Harper[END_REF]), there exists b 0 > 0 s. Let us turn to the proof of (1.31). Under the above conditions, the Riesz integral formula allows to write for any l ∈ {1, . . . , τ }:

∀|b| ≤ b, λ l (b)Π l (b) = H P (b)Π l (b) = i 2π γ l dξ ξ(H P (b) -ξ) -1 .
Since dim Ran(Π l (b)) = 1 by stability of the λ l 's, then for any n 0 ∈ {1, . . . , τ }:

∀|b| ≤ b, n0 l=1 λ l (b) = i 2π Tr L 2 (R 3 ) ∪ n 0 l=1 γ l dξ ξ(H P (b) -ξ) -1 . (3.23)
Next, the proof of (1.31) follows from the following result which is concerned with the quantity in the r.h.s. of (3.23) (seen as a function of b): Proposition 3.11 Let {g θ,w , θ ∈ C}, w = 0, 1 be the families of functions defined in (3.20). Then, there exists a neighborhood I of b = 0 s.t. ∀θ ∈ C and ∀w ∈ {0, 1}, the following map:

b → F θ,w (b) := i 2π Tr L 2 (R 3 ) ∪ n 0 l=1 γ l dξ g θ,w (ξ)(H P (b) -ξ) -1 , (3.24)
is twice differentiable on I. Moreover, its second derivative at b = 0 reads as: The rest of this section is devoted to the proof of Proposition 3.11. It is essentially based on the so-called gauge invariant magnetic perturbation theory, see [START_REF] Cornean | On the essential spectrum of two-dimensional periodic magnetic Schrödinger operators[END_REF][START_REF] Cornean | Two dimensional magnetic Schrödinger operators: width of minibands in the tight-binding approximation[END_REF][START_REF] Cornean | On the magnetization of a charged Bose gas in the canonical ensemble[END_REF][START_REF] Nenciu | On asymptotic perturbation theory for quantum mechanics: Almost invariant subspaces and gauge invariant magnetic perturbation theory[END_REF][START_REF] Briet | Locating the spectrum for magnetic Schrödinger and Dirac operators, Commun[END_REF][START_REF] Briet | Diamagnetic expansions for perfect quantum gases[END_REF][START_REF] Cornean | The Faraday effect revisited: General theory[END_REF][START_REF] Briet | Diamagnetic expansions for perfect quantum gases II: Uniform bounds[END_REF][START_REF] Cornean | The Faraday effect revisited: Thermodynamic limit[END_REF][START_REF] Cornean | On the Lipschitz continuity of spectral bands of Harper-like and magnetic Schrödinger operators[END_REF][START_REF] Briet | A rigorous proof of the Landau-Peierls formula and much more[END_REF][START_REF] Briet | A rigorous approach to the magnetic response in disordered systems[END_REF][START_REF] Savoie | On the zero-field orbital magnetic susceptibility of Bloch electrons in graphene-like solids: Some rigorous results[END_REF] and references therein.

d 2 F θ,w db 2 (b = 0) := i π Tr L 2 (R 3 ) ∪ n 0 l=1 γ l dξ g θ,w (ξ)(H P -ξ) -1 [T P,1 (ξ)
Before starting, let us introduce some notation. Define ∀(x, y) ∈ R 6 the magnetic phase φ as:

φ(x, y) := 1 2 e 3 • (y × x) = -φ(y, x), e 3 := (0, 0 , 1). (3.26)
By [47, Thm. B.7.2], ∀b ∈ R and ∀ξ ∈ ̺(H P (b)), the resolvent (H P (b)-ξ) -1 is an integral operator with integral kernel (H

P (b) -ξ) -1 (• , • ) jointly continuous on R 6 \ D. Introduce on L 2 (R 3
) the operators T P,j (b, ξ), j = 1, 2 via their kernel respectively defined on R 6 \ D by:

T P,1 (x, y; b, ξ) := a(x -y) • (i∇ x + ba(x))(H P (b) -ξ) -1 (x, y), T P,2 (x, y; b, ξ) := 1 2 a 2 (x -y)(H P (b) -ξ) -1 (x, y).
From [12, Eq. (2.9)] together with [12, Lem. 2.4], ∀η > 0 there exists a constant ϑ = ϑ(η) > 0 and a polynomial p(• ) s.t. ∀b ∈ R, ∀ξ ∈ C satisfying dist(ξ, σ(H P (b))) ≥ η and ∀(x, y) ∈ R 6 \ D:

(H P (b) -ξ) -1 (x, y) ≤ p(|ξ|) e -ϑ ξ |x-y| |x -y| , ϑ ξ := ϑ(1 + |ξ|) -1 , (3.27)
|T P,j (x, y; b, ξ)| ≤ p(|ξ|)(1 + |b|) 3 e -ϑ ξ |x-y| |x -y| , j = 1, 2. (3.28)
By the Shur-Holmgren criterion, the T P,j (ξ)'s are bounded on L 2 (R 3 ):

(H P (b) -ξ) -1 ≤ p(|ξ|), T P,j (b, ξ) ≤ p(|ξ|)(1 + |b|) 3 , (3.29)
for another polynomial p(• ). Let j ∈ N * and i := (i 1 , . . . , i j ) ∈ {1, 2} j . For any k ∈ N * with k ≥ j ≥ 1, let χ k j (i) be the characteristic function given by:

χ k j (i) := 1 if i 1 + • • • + i j = k, 0 otherwise , 1 ≤ j ≤ k.
Restricting to k ∈ {1, 2}, define ∀m ∈ {0, 1} and ∀b ∈ R the function on R 6 :

T m P,k (x, y; b, ξ) := k j=1 (-1) j i∈{1,2} j χ k j (i) R 3 dz 1 • • • R 3 dz j (iφ(z j , y) -iφ(z j , x)) m × (H P (b) -ξ) -1 (x, z 1 )T P,i1 (z 1 , z 2 ; b, ξ) • • • T P,ij (z j , y; b, ξ), (3.30)
where by convention, we set 0 0 = 1. Note that, due to the antisymmetry of φ, the terms in the r.h.s. of (3.30) containing the magnetic phases vanish when x = y. Moreover, ∀η > 0, ∀b ∈ R and ∀ξ ∈ C satisfying dist(ξ, σ(H

P (b))) ≥ η, T m P,k (• , • ; b, ξ) is jointly continuous on R
|T m P,k (x, y; b, ξ)| ≤ p(|ξ|)(1 + |b|) 6 × (|x| m + |y| m ) if x = y, 1 if x = y, (3.31)
where, in the case of x = y, we used the rough estimate:

∀(x, y) ∈ R 6 , |φ(x, y)| ≤ |y||x -y|.
Remark 3.12 In view of (3.30), one has on R 3 :

T 0 P,1 (x, x; b, ξ) = - R 3
dz (H P (b)ξ) -1 (x, z)T P,1 (z, x; b, ξ), The first one deals with the regularity of the integral kernel of (H P (b)-ξ) -1 , seen as a function of the b-variable, in a neighborhood of b = 0:

T 1 P,1 (x, x; b, ξ) = 0, T 0 P,2 (x, x; b, ξ) = - R 3 dz (H P (b) -ξ) -1 (x, z)T P,2 (z, x; b, ξ) + R 3 dz 1 R 3 dz 2 (H P (b) -ξ) -1 (x, z 1 )T P,1 (z 1 , z 2 ; b, ξ)T P,1 (z 2 , x; b, ξ).
Lemma 3.13 Let K ⊂ (̺(H P ) ∩ {ζ ∈ C : ℜζ < 0}) be a compact subset. Let b K > 0 s.t. ∀|b| ≤ b K , K ⊂ (̺(H P (b)) ∩ {ζ ∈ C : ℜζ < 0}). Then, ∀ξ ∈ K and ∀(x, y) ∈ R 6 \ D, the map b → (H P (b) -ξ) -1 (x, y) is twice differentiable on the interval (-b K , b K ). Furthermore, its first two derivatives at b 0 ∈ (-b K , b K ) read as: ∂ s ∂b s (H P (b) -ξ) -1 (x, y) b=b0 = (iφ(x, y)) s (H P (b 0 )-ξ) -1 (x, y)+s s k=1 T s-k P,k (x, y; b 0 , ξ), s = 1, 2,
where the functions

T m P,k (• , • ; b 0 , ξ), k = 1, 2, m = 0, 1 are defined in (3.30).
The second one deals with the 'uniform' exponential localization of the eigenfunctions associated with the λ l (b)'s, l = 1, . . . , τ for b sufficiently small: Lemma 3.14 Let λ l , l ∈ {1, . . . , τ } be a simple eigenvalue of H P in (-∞, 0). Let b sufficiently small s.

∀|b| ≤ b, ∀x ∈ R 3 , max {|Φ l (x; b)| , |∂ x k Φ l (x; b)|} ≤ Ce -c|x| . (3.32)
The exponential decay for the Φ l (• ; b)'s is a well-known result, see e.g. [START_REF] Iftimie | Décroissance exponentielle des fonctions propres des opérateurs de Schrödinger avec des champs magnétiques[END_REF]Thm. 4.4] and also [START_REF] Iftimie | Eigenfunctions decay for magnetic pseudodifferential operators[END_REF]Thm. 1.10], [START_REF] Fournais | Spectral Methods in Surface Superconductivity[END_REF]Sec. 7.2]. For b sufficiently small, we stress the point that all the constants can be chosen b-independent. As a result of Lemma 3.14, one directly gets as a corollary of Proposition 3.9: Then ∀θ ∈ C, ∀w ∈ {0, 1} and ∀ς ∈ {1, . . . , τ }, there exist two constants C = C(θ, ς) > 0 and c > 0 s.t. ∀|b| ≤ b, ∀x ∈ R 3 and for any j ∈ {1, 2}:

max i 2π Γς dξ g θ,w (ξ) (H P (b) -ξ) -1 T P,j (b, ξ) (x, x) , i 2π Γς dξ g θ,w (ξ) (H P (b) -ξ) -1 T P,1 (b, ξ)T P,1 (b, ξ) (x, x) ≤ C(1 + |b|) 6 e -c|x| , (3.33)
where g θ,w are the maps in (3.20) and Γ ς a contour as in Proposition 3.9.

Note that the presence of the factor (1 + |b|) 6 in the upper bound in (3.33) comes from the estimate (3.28) (the estimate in (3.27) is b-independent).

We are now ready for: Proof of Proposition 3.11. Let θ ∈ C and w ∈ {0, 1}. Let b > 0 s.t. ∀|b| ≤ b, (3.23) holds. Let us first prove that:

∀|b| ≤ b, F θ,w (b) := i 2π ∪ n 0 l=1 γ l dξ g θ,w (ξ)(H P (b) -ξ) -1 ,
has an integral kernel jointly continuous on R 6 . By the Dunford-Pettis theorem, F θ,w (b) is an integral operator since (H

P (b) -ξ) -1 is bounded from L 2 (R 3 ) to L ∞ (R 3
) by some polynomial in |ξ|, see (3.27). Let ξ 0 < inf σ(H P ) and large enough s.t. ξ 0 < min 1≤l≤τ {γ l ∩ R}. By the first resolvent equation:

F θ,w (b) = i 2π ∪ n 0 l=1 γ l dξ g θ,w (ξ) (H P (b) -ξ 0 ) -1 + i 2π ∪ n 0 l=1 γ l dξ g θ,w (ξ)(ξ -ξ 0 )(H P (b) -ξ) -1 (H P (b) -ξ 0 ) -1 .
By the Cauchy-Goursat theorem, the first term in the above r.h.s. is identically zero. The second term has a jointly continuous integral kernel on R 6 due to [12, Lem. A.1] along with (3.27).

Denoting by F θ,w (• , • ; b) the integral kernel of F θ,w (b), its diagonal part reads ∀b ∈ (-b, b) as:

∀x ∈ R 3 , F θ,w (x, x; b) := i 2π ∪ n 0 l=1 γ l dξ g θ,w (ξ)(H P (b) -ξ) -1 (x, y) y=x . (3.34)
Next, from Lemma 3.13 together with (3.31), one can prove that for any x ∈ R 3 , the map b → F θ,w (x, x; b) is twice differentiable on (-b, b). Moreover, its first two derivatives at b 0 ∈ (-b, b) satisfy (below s = 1, 2):

∀x ∈ R 3 , ∂ s F θ,w ∂b s (x, x; b 0 ) := s i 2π ∪ n 0 l=1 γ l dξ g θ,w (ξ) s k=1
T s-k P,k (x, x; b 0 , ξ).

Here, we used that φ(x, x) = 0. Finally, let 0 < b ≤ b s.t. ∀|b| ≤ b and for any l ∈ {1, . . . , n 0 }, the (normalized) eigenfunction associated with λ l (b) obeys an estimate of type (3.32). From the explicit expressions in Remark 3.12 together with the estimate (3.33), then for any compact subset K ⊂ (-b, b), there exist two constants c > 0 and

C = C(n 0 , θ, K) > 0 s.t. ∀x ∈ R 3 , sup b∈K ∂ s F θ,w ∂b s (x, x; b) ≤ Ce -c|x| , s = 1, 2.
Since the upper bound belongs to L 1 (R 3 ), the proposition follows.

Proof of (iii).

Let us recall some notation. Under the assumptions (A r )-(A m )-(A nd ), let {λ l } τ l=1 be the set of eigenvalues of H P in (-∞, 0) counting in increasing order. For any l ∈ {1, . . . , τ }, denote by Φ l the normalized eigenfunction associated with λ l , and by Π l = |Φ l Φ l | the orthogonal projection onto Φ l . We define Π ⊥ l := ½ -Π l . From [5, Thm. Let us turn to the proof of (iii). Let l ∈ {1, . . . , τ } and K be a compact neighborhood of λ l . Let 0 < b ≤ b 1 s.t. ∀|b| ≤ b, K ∩ σ(H P (b)) \ λ l (b) = ∅ and λ l (b) ∈ K. From the Feshbach formula in [START_REF] Feshbach | Unified theory of nuclear reactions[END_REF] and under our conditions, ∀|b| ≤ b and small enough, λ l (b) is the unique number ζ near λ l for which:

(λ l -ζ)Π l + Π l W (b)Π l -Π l W (b) Π ⊥ l (H P + W (b) -ζ)Π ⊥ l -1 W (b)Π l , (3.35)
is not invertible. Note that W (b)Π l is bounded on L 2 (R 3 ), and its operator norm behaves like O(|b|). This follows from the exponential localization of Φ l and ∇Φ l in (3.32). We point out that W (b)Π l remains bounded even with an exponential weight. More precisely, there exists ε 0 > 0 and a constant C > 0 s.t.

∀0 < ε ≤ ε 0 , W (b)Π l e ε • ≤ C|b|, • := 1 + | • | 2 . (3.36)
Now, let us justify that the operator Π ⊥ l (H P (b)ξ)Π ⊥ l is invertible ∀|b| ≤ b and small enough, and ∀ξ ∈ K. To do so, introduce the Sz-Nagy transformation in [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF]Sec. I.4.6] corresponding to the pair of projections Π l (b), Π l :

U (b) = ½ -(Π l (b) -Π l ) 2 -1 2 Π l (b)Π l + Π ⊥ l (b)Π ⊥ l . ( 3
U (b)Π ⊥ l (H P (b) -ξ)Π ⊥ l U * (b) = Π ⊥ l (b)U (b)(H P (b) -ξ)U * (b)Π ⊥ l (b), (3.39)
i.e. Π ⊥ l (H P (b)ξ)Π ⊥ l is unitarily equivalent with the operator in the r.h.s. of (3.39). Next, from the asymptotic perturbation theory in [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF]Sec. VIII.2.4]:

U (b) = ½ + Π l W (b)(H P -ξ) -1 (½ -Π l ) -(½ -Π l )(H P -ξ) -1 W (b)Π l + u 1 (b), U * (b) = ½ + (½ -Π l )(H P -ξ) -1 W (b)Π l -Π l W (b)(H P -ξ) -1 (½ -Π l ) + u 2 (b),
where the u j (b)'s, j = 1, 2 are operators satisfying b -1 u j (b) → 0 in the strong sense. Putting these asymptotic expansions into (3.39), then we obtain:

U (b)Π ⊥ l (H P (b) -ξ)Π ⊥ l U * (b) = Π ⊥ l (b)(H P (b) -ξ)Π ⊥ l (b) [½ + Υ(ξ, b)] ,
where Let us get back to the quantity in (3.35). By using the scalar product, ∀l ∈ {1, . . . , τ } and ∀|b| ≤ b and small enough, λ l (b) has to obey the equation:

λ l (b) = λ l + Φ l , W (b)Φ l -W (b)Φ l , Π ⊥ l (H P (b) -λ l (b)) Π ⊥ l -1 W (b)Φ l . (3.40)
By iterating the identity in (3.40), one obtains under the same conditions:

λ l (b) = λ l + Φ l , W (b)Φ l -W (b)Φ l , Π ⊥ l (H P -λ l ) Π ⊥ l -1 W (b)Φ l + O |b| 3 . (3.41)
To control the behavior in b of the remainder term, we used that:

e -ε • Π ⊥ l (H P (b) -λ l (b)) Π ⊥ l -1 -Π ⊥ l (H P -λ l ) Π ⊥ l -1 e -ε • = O(|b|),
with 0 < ε ≤ ε 0 as in (3.36). It can be proved by using an asymptotic expansion for the projection Π l (b) together with the application of the gauge invariant magnetic perturbation theory for the kernel of the unperturbed projector as we did for the kernel of the resolvent in Sec. 3.2. Note that the fourth term involved in the expansion (3.41) can be identified with:

-Φ l , W (b)Φ l W (b)Φ l , Π ⊥ l (H P -λ l ) Π ⊥ l -2 W (b)Φ l + W (b)Φ l , Π ⊥ l (H P -λ l ) Π ⊥ l -1 W (b) Π ⊥ l (H P -λ l ) Π ⊥ l -1 W (b)Φ l .
Since λ l (• ) has an asymptotic series expansion for b small enough, then one obtains from (3.41):

d 2 λ l db 2 (b = 0) = Φ l , a 2 Φ l -2 a • (-i∇)Φ l , Π ⊥ l (H P -λ l ) Π ⊥ l -1 a • (-i∇)Φ l . (3.42)
From (3.42), the proof of (1.32) directly follows from (1.31).

4 Appendix.

C (P ) β = Γς,β ∪ γ (2) ς,β ∪ γ (1) ς,β , ς ∈ {1, . . . , τ }, (4.2) γ (1)
ς,β := ℜξ ∈ [ω + ς , ∞), ℑξ = ± π 2β ∪ ℜξ = ω + ς , ℑξ ∈ - π 2β , π 2β , γ (2)
ς,β := ℜξ ∈ [ω - ς , ω + ς ], ℑξ = ± π 2β ∪ ℜξ = ω ± ς , ℑξ ∈ - π 2β , π 2β , Γς,β := ℜξ ∈ [δ P , ω - ς ], ℑξ = ± π 2β ∪ ℜξ = δ P , ℑξ ∈ - π 2β , π 2β
∪ ℜξ = ω - ς , ℑξ ∈ - π 2β , π 2β
,
where ω - ς := 19λς +5λς+1 24 and ω + ς := 5λς +19λς+1

24

, with the convention λ τ +1 = 0.

Let us turn to the proof of Proposition 3.6: Proof of Proposition 3.6. Let ς ∈ {1, . . . , τ }. From (3.13) and (4.2), ∀µ ∈ I ς : ς,β is a strict subset of the holomorphic domain of the integrand in (4.3). Therefore, the Cauchy-Goursat theorem yields:

X R (β, µ, b = 0) = q c 2 1 β 1 |Ω R | i π Γς,β ∪γ (2) ς,β ∪γ (1) ς,β dξ f(β, µ; ξ)Θ R (ξ), (4.3)
γ (2) ς,β dξ f(β, µ; ξ)Θ R (ξ) = 0.
Now, the contours Γς,β and γ (1) ς,β can be deformed respectively in Γς,1 and γ (1) ς,1 (set β = 1 in their definition) due to the location of the interval I ς . In the wake of the deformation of the contours, then under the conditions of Lemma 2.1, from (2.16) and (3.12) with Ξ = P and by setting η = 1, there exists a polynomial p(• ) independent of β (and R) s.t.

|Θ R (ξ)| ≤ p(|ξ|)|Ω R |. (4.4)
It remains to use the Lebesgue's dominated convergence theorem which yields:

lim β↑∞ 1 β Γς,1 dξ f(β, µ; ξ)Θ R (ξ) = Γς,1 dξ (µ -ξ)Θ R (ξ), lim β↑∞ 1 β γ (1) ς,1 dξ f(β, µ; ξ)Θ R (ξ) = 0,
where we used the pointwise convergence:

lim β↑∞ 1 β f(β, µ; ξ) = (µ -ξ)χ (-∞,µ) (ℜξ). ( 4.5)
The uniform convergence on compact subsets K ⊂ I ς is straightforward. ς,β in some β-independent contours is crucial to make the estimate in (4.4) β-independent.

Next, let us turn to the proof of Proposition 3.7. To do that, we need beforehand to introduce some operators. Recall that R R (ξ) and R R (ξ), R ≥ R 0 are respectively defined in (2.6) and (2.12) (R 0 is defined through (2.1)). Introduce ∀R ≥ R 0 , the bounded operators T R,j (ξ) and T R,j (ξ), j = 1, 2 on L 2 (R 3 ) generated via their kernel respectively defined on R 6 \ D by: T R,1 (x, y; ξ) := a(xy)

• (i∇ x )(R R (ξ))(x, y), T R,2 (x, y; ξ) := 1 2 a 2 (x -y)(R R (ξ))(x, y), T R,1 (x, y; ξ) := a(x -y) • (i∇ x )(R R (ξ))(x, y), T R,2 (x, y; ξ) := 1 2 a 2 (x -y)(R R (ξ))(x, y).
The proof relies on the following four lemmas whose proofs can be found in Appendix, Sec. 4.4. Recall that R 0 , R 1 ≥ 1 are respectively defined through (2.1)- (2.14). With the shorthand notation introduced in (2.15): Lemma 4.2 Let 0 < α < 1. For every η > 0, there exists a constant ϑ = ϑ(η) > 0 and a polynomial p(• ) s.t. ∀R ≥ R 0 and ∀ξ ∈ C satisfying dist(ξ, σ(H R ) ∩ σ(H P )) ≥ η:

1 |Ω R | Tr L 2 (R 3 ) χ ΩR (H R -ξ) -1 [T R,1 (ξ)T R,1 (ξ) -T R,2 (ξ)] χ ΩR -Tr L 2 (R 3 ) {χ ΩR R R (ξ) [T R,1 (ξ)T R,1 (ξ) -T R,2 (ξ)] χ ΩR } ≤ p(|ξ|)e -ϑ ξ R α . (4.6)
Lemma 4.3 Let 0 < α < 1. For every η > 0, there exists a constant ϑ = ϑ(η) > 0 and a polynomial p(• ) s.t. ∀R ≥ max{R 0 , R 1 } and ∀ξ ∈ C satisfying dist(ξ, σ(H R ) ∩ σ(H P )) ≥ η:

1 |Ω R | Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,2 (ξ)χ ΩR } -Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,2 (ξ)χ ΩR } ≤ p(|ξ|)e -ϑ ξ R α .
(4.7)

Lemma 4.4 Let 0 < α < 1. For every η > 0, there exists a constant ϑ = ϑ(η) > 0 and a polynomial p(• ) s.t. ∀R ≥ max{R 0 , R 1 } and ∀ξ ∈ C satisfying dist(ξ, σ(H R ) ∩ σ(H P )) ≥ η:

1 |Ω R | Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,1 (ξ)T R,1 (ξ)χ ΩR } -Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,1 (ξ)T R,1 (ξ)χ ΩR } ≤ p(|ξ|)e -ϑ ξ R α . (4.8)
Lemma 4.5 Let 0 < α < 1. For every η > 0, there exists a constant ϑ = ϑ(η) > 0 and a polynomial p(• ) s.t. ∀R ≥ R 0 and ∀ξ ∈ C satisfying dist(ξ, σ(H R ) ∩ σ(H P )) ≥ η:

1 |Ω R | Tr L 2 (R 3 ) {χ ΩR R R (ξ) [T R,1 (ξ)T R,1 (ξ) -T R,2 (ξ)] χ ΩR } -Tr L 2 (R 3 ) χ ΩR (H P -ξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] χ ΩR ≤ p(|ξ|)e -ϑ ξ R α .
(4.9)

We are now ready for: Proof of Proposition 3.7. In view of (3.9) and (3.13), then ∀β > 0 and ∀µ ∈ R the quantity ∆ R (β, µ) in (3.14) can be rewritten for R sufficiently large as:

∆ R (β, µ) = q c 2 1 β 1 |Ω R | i π C (P ) β dξ f(β, µ; ξ)K R (ξ), (4.10) K R (ξ) := Tr L 2 (R 3 ) χ ΩR (H R -ξ) -1 [T R,1 (ξ)T R,1 (ξ) -T R,2 (ξ)] -(H P -ξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] χ ΩR . (4.11)
By R sufficiently large, we mean R ≥ R 2 ≥ 1 chosen in such a way that ∀R ≥ R 2 , inf σ(H R ) ≥ λ 1 -1 2 . Such R 2 exists since inf σ(H R ) has to coincide with inf σ(H P ) when R ↑ ∞, see Lemma 3.4. Next, let us estimate the quantity in (4.11). In view of the definition of C (P ) β in (3.15), set η := min{ 1 2 , π 2β } > 0. Pick 0 < α < 1 and set R 3 := max{R 0 , R 1 , R 2 } ≥ 1. From Lemmas 4.2-4.5, there exists a constant ϑ = ϑ(η) > 0 and a polynomial p(• ) s.t. ∀R ≥ R 3 :

∀ξ ∈ C (P ) β , |K R (ξ)| ≤ p(|ξ|)|Ω R |e -ϑ ξ R α , ϑ ξ := ϑ(1 + |ξ|) -1 .
(4.12)

We point out that ϑ and p(• ) in (4.12) are β-dependent (at least for large β).

Let ς ∈ {1, . . . , τ }. From now on, we limit the domain of µ to the interval I ς . Without loss of generality, we give the explicit form (4.1) to I ς . Justified by Lemma 3.4, we assume that R ≥ R 3 and large enough so that [ 5λς +λς+1 6

, λς +5λς+1 6 ] ∩ σ(H R ) = ∅ (we set λ τ +1 := 0). Let ω ± ς be the reals defined below (4.2). Since

I ς ⊂ [ω - ς , ω + ς ] ⊂ [ 5λς +λς+1 6 , λς +5λς+1 6
], then ∀R ≥ R 3 and large enough,

I ς ∩ σ(H R ) = ∅ and Γς,β , γ (1)
ς,β ∩ σ(H R ) = ∅. Here Γς,β and γ

(1) ς,β are the contours coming from the decomposition of C (P ) β in (4.2). From the foregoing, it remains to mimic the proof of Proposition 3.6 to perform the zero-temperature limit. One arrives, for any compact subset K ⊂ I ς , at:

∆ R (µ) := lim β↑∞ ∆ R (β, µ) = q c 2 1 |Ω R | i π Γς,1 dξ (µ -ξ)K R (ξ),
uniformly in µ ∈ K. We emphasize that, as in the proof of Proposition 3.6, the deformation of Γς,β and γ

ς,β in some β-independent contours makes the estimate in (4.12) β-independent on these contours, see above the definition of the η. Finally, it remains to use (4.6)-(4.9) along with (4.12) which lead to the existence of two constants c, C > 0 s.t. ∀R ≥ R 3 and large enough:

∀µ ∈ I ς , |∆ R (µ)| ≤ C(1 + |µ|)e -cR α .
4.2 Proof of Propositions 3.9 and 3.10.

Proof of Proposition 3.9. Let {λ l } τ l=1 , τ ∈ N * be the set of eigenvalues of H P in (-∞, 0) counting in increasing order. Let ς ∈ N * with 1 ≤ ς ≤ τ . For any l ∈ {1, . . . , ς}, denote by Π l the orthogonal projection onto the eigenvector corresponding to the eigenvalue λ l . Recall that Π l is given by a Riesz integral:

Π l = i 2π γ l dξ (H P -ξ) -1 ,
where γ l is any positively oriented closed simple contour surrounding the eigenvalue λ l , while letting outside the rest of the spectrum of H P . Denote by Π ℵ := ς l=1 Π l . Introduce the following decomposition of the resolvent:

∀ξ ∈ ̺(H P ), (H P -ξ) -1 = (H P -ξ) -1 ℵ + (H P -ξ) -1 ⊥ , (4.13)
(H P -ξ) -1 ℵ := Π ℵ (H P -ξ) -1 = ς l=1 1 λ l -ξ Π l , (4.14)
(H P -ξ) -1 ⊥ := (½ -Π ℵ )(H P -ξ) -1 = - i 2π ∪ ς l=1 γ l dζ ζ -ξ (H P -ζ) -1 . (4.15)
Due to (4.14), we have in the kernels sense on R 6 :

(H P -ξ) -1 ℵ (x, y) = ς l=1 1 λ l -ξ Φ l (x)Φ l (y), (4.16)
where Φ l denotes the normalized eigenfunction associated with the eigenvalue λ l . Note that under our conditions, the Φ l 's with = 1, . . . , ς satisfy the following, see e.g. [START_REF] Daumer | Equation de Schrödinger avec champ électrique périodique et champ magnétique constant dans l'approximation du tight-binding[END_REF]Eq. (5.8)]. There exists a constant C > 0 s.t.

∀x ∈ R 3 , max {|Φ l (x)|, |∂ x k Φ l (x)|} ≤ Ce - √ |λ l ||x| , k ∈ {1, 2, 3}. (4.17)
Moreover, denote by T ℘ P,j (ξ), j = 1, 2 and ℘ = ℵ, ⊥, the operators on L 2 (R 3 ) defined via their kernel as in (1.27)-(1.28) but with (H Pξ) -1 ℘ instead of (H Pξ) -1 . Due to (4.15), the kernel of T ⊥ P,j (ξ), j = 1, 2 still obey (3.12). Let θ ∈ C, w ∈ {0, 1} and j ∈ {1, 2}. Due to the decomposition in (4.13), introduce for any (℘ 1 , ℘ 2 ) ∈ {ℵ, ⊥} 2 the bounded operator on L 2 (R 3 ):

M (j),℘1,℘2 θ,w := i 2π Γς dξ g θ,w (ξ)(H P -ξ) -1 ℘1 T ℘2 P,j (ξ). (4.18)
It is clearly an integral operator with a jointly continuous integral kernel on R 6 . We denote it by M (j),℘1,℘2 θ,w (• , • ). The first part of the proof consists in showing the existence of two constants c > 0 and C = C(θ, ς) > 0 s.t.

∀x ∈ R 3 , (℘1,℘2)∈{ℵ,⊥} 2 M (j),℘1,℘2 θ,w (x, x) ≤ Ce -c|x| .
At first, due to the location of Γ ς , M (j),⊥,⊥ θ,w (x, x) = 0 on R 3 by the Cauchy-Goursat theorem. Secondly, from (4.16) followed by the residue theorem:

∀x ∈ R 3 , M (j),ℵ,⊥ θ,w (x, x) = ς l=1 g θ,w (λ l )Φ l (x) R 3 dz Φ l (z)T ⊥ P,j (z, x; λ l ).
Due to (4.17) and (3.12), use next that there exists a constant C > 0 s.t.

∀1 ≤ l ≤ ς, sup

x∈R 3 R 3 dz Φ l (z)T ⊥ P,j (z, x; λ l ) ≤ C.
Then, it follows that |M

∀ν ≥ 0, ∀µ > 0, t ν e -µt ≤ Ce -µ 2 t , t ≥ 0, (4.19)
for another C = C(ν, µ) > 0. The last term we have to treat reads as:

M (1),ℵ,ℵ θ,1 (x, x) := ς l1,l2=1 Φ l1 (x) R 3 dz Φ l1 (z)a(z -x) • (i∇ z Φ l2 )(z)Φ l2 (x), M (2),ℵ,ℵ θ,1 (x, x) := ς l1,l2=1 Φ l1 (x) R 3 dz Φ l1 (z) 1 2 a 2 (z -x)Φ l2 (z)Φ l2 (x),
if w = 1 and j = 1 or j = 2 respectively; and if w = 0 regardless of j = 1, 2:

M (j),ℵ,ℵ θ,0 (x, x) = 0;
where we used the following identity provided by the residue theorem:

ς l1,l2=1 Γς dξ g θ,w (ξ) (λ l1 -ξ)(λ l2 -ξ) = -2iπ ς l1,l2=1 if w = 1, 0 if w = 0 .
In the case of w = 1, it is enough to use (4.17), the rough estimate |a(x-y)| ≤ (|x|+|y|) combined with (4.19). This leads to |M (j),ℵ,ℵ θ,1 (x, x)| ≤ Ce -c|x| ∀x ∈ R 3 and j = 1, 2, for another constant C = C(θ, ς) > 0 and c > 0.

Afterwards, we turn to the second part of the proof. Under the conditions of (4.18), introduce ∀(℘ 1 , ℘ 2 , ℘ 3 ) ∈ {ℵ, ⊥} 3 the operator on L 2 (R 3 ):

N ℘1,℘2,℘3 θ,w := i 2π Γς dξ g θ,w (ξ)(H P -ξ) -1 ℘1 T ℘2 P,1 (ξ)T ℘3 P,1 (ξ).
It is an integral operator with a jointly continuous integral kernel on R 6 . We denote it by N ℘1,℘2,℘3 θ,w (• , • ). Let us show that there exist c, C > 0 s.t.

∀x ∈ R 3 , (℘1,℘2,℘3)∈{ℵ,⊥} 3 N ℘1,℘2,℘3 θ,w (x, x) ≤ Ce -c|x| .
Firstly, N ⊥,⊥,⊥ θ,w (x, x) = 0 on R 3 by the Cauchy-Goursat theorem. Also, by a direct calculation, N ℵ,ℵ,ℵ θ,w (x, x) = 0 since the residue theorem leads to:

ς l1,l2,l3=1 Γς dξ g θ,w (ξ) (λ l1 -ξ)(λ l2 -ξ)(λ l3 -ξ) = 0.
It remains to treat six terms. Let us treat the trickiest one (we make it clear hereafter) which is N ⊥,ℵ,⊥ θ,w (x, x). From the residue theorem:

N ⊥,ℵ,⊥ θ,w (x, x) = ς l=1 g θ,w (λ l ) R 3 dz 1 R 3 dz 2 (H P -λ l ) -1 ⊥ (x, z 1 ) × a(z 1 -z 2 ) • (i∇ z1 Φ l )(z 1 )Φ l (z 2 )T ⊥ P,1 (z 2 , x; λ l ). Define for any 1 ≤ l ≤ ς, k ∈ {1, 2} and ∀(x, z 1 , z 2 ) ∈ R 9 with x = z 1 = z 2 : J (k) l (x, z 1 , z 2 ) := (H P -λ l ) -1 ⊥ (x, z 1 )(i∇ z1 Φ l )(z 1 ) |z k | Φ l (z 2 )T ⊥ P,1 (z 2 , x; λ l ) .
From (2.16), (3.12) and (4.17), there exist two constants C, c > 0 s.t.

J

(1)

l (x, z 1 , z 2 ) ≤ C e -c|x-z1| |x -z 1 | |z 1 |e - √ |λ l ||z1| e - √ |λ l ||z2| e -c 2 |z2-x| e -c 2 |z2-x| |z 2 -x| .
Using that e -min{ c 2 ,

√ |λ l |}(|z2|+|z2-x|) ≤ e -min{ c 2 , √ |λς |}|x| ∀x ∈ R 3 , along
∀x ∈ R 3 , R 3 dz 1 R 3 dz 2 J (1) l (x, z 1 , z 2 ) ≤ Ce -min{ c 2 , √ |λς |}|x| . (4.20)
Clearly, the upper bound in (4.20) still holds true when replacing J

(1) l with J

l . Hence, one arrives at |N ⊥,ℵ,⊥ θ,w (x, x)| ≤ Ce -c|x| ∀x ∈ R 3 , for another constant C = C(θ, ς) > 0 and c > 0. We do not treat the other terms which are simpler. Note that they have all the form N ℵ,℘2,℘3 θ,w

(x, x) with (℘ 2 , ℘ 3 ) ∈ {ℵ, ⊥} 2 and ℘ 2 , ℘ 3 = ℵ, or N ℘1,℘2,ℵ θ,w (x, x) with (℘ 1 , ℘ 2 ) ∈ {ℵ, ⊥} 2 and ℘ 1 , ℘ 2 = ℵ.
Therefore, they can always be written on R 3 , via the residue theorem, as ς l=1 Φ l (x)A l (x) with sup x∈R 3 |A l (x)| ≤ cste. Thus, the announced exponential decay in |x| results only from (4.17).

Proof of Proposition 3.10. The proof we give requires the notation introduced in Sec. 3.2. We do not recall them, and we refer to the beginning of Sec. 3.2. Let ς ∈ {1, . . . , τ }. From the Riesz integral formula (for orthogonal projections) together with the definition (3.24), one has:

∀|b| ≤ b, ς l=1 Tr L 2 (R 3 ) {Π l (b)} = i 2π Tr L 2 (R 3 ) ∪ ς l=1 γ l dξ (H P (b) -ξ) -1 =: F 1,0 (b).
Here, the b and the γ l 's are the same as the ones appearing in (3.23). Since the map b → F 1,0 (b) is twice differentiable in a neighborhood of b = 0 (see Proposition 3.11), then by using the expression for the second derivative at b = 0 given in (3.25) (with θ = 1 and w = 0), one gets the identity:

d 2 db 2 ς l=1 Tr L 2 (R 3 ) {Π l (b)} b=0 = i π Tr L 2 (R 3 ) ∪ ς l=1 γ l dξ (H P -ξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] .
But the quantity in the above l.h.s. is zero. Indeed, by stability of the eigenvalues of H P in (-∞, 0), one has dim RanΠ l (b) = 1 and thus the sum is a b-independent quantity (equals to ς). (-1) j i∈{1,2} j χ k j (i)

R 3 dz 1 • • • R 3
dz j e iδb(φ(x,z1)+•••+φ(zj ,y))

× (H P (b 0 )ξ) -1 (x, z 1 )T P,i1 (z 1 , z 2 ; b 0 , ξ) • • • T P,ij (z j , y; b 0 , ξ), where the function T k-m P,m (• , • ; b 0 , ξ), 2 ≥ k ≥ m ≥ 1 is defined in (3.30), and the last term stands for the remainder term. We mention that we have used the explicit expressions in Remark 3.12 to rewrite the second term in the r.h.s. of (4.31) coming from (4.30). Note also that, by construction, T This can be extended to (-b K , b K ). The lemma follows by induction.

4.4 Proof of Lemmas 4.2-4.5.

Throughout this section, (I 2 (L 2 (R 3 )), • I2 ) and (I 1 (L 2 (R 3 )), • I1 ) denote the Banach space of Hilbert-Schmidt (H-S) and trace class operators on L 2 (R 3 ) respectively.

Proof of Lemma 4.2. Let us denote:

Y R,1 (ξ) := |Ω R | -1 Tr L 2 (R 3 ) χ ΩR (H R -ξ) -1 T R,1 (ξ)T R,1 (ξ)χ ΩR , Y R,2 (ξ) := |Ω R | -1 Tr L 2 (R 3 ) χ ΩR (H R -ξ) -1 T R,2 (ξ)χ ΩR .
By replacing (H Rξ) -1 with the r.h.s. of (2.7) in Y R,j (ξ), j = 1, 2 then:

Y R,1 (ξ) = |Ω R | -1 Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,1 (ξ)T R,1 (ξ)χ ΩR } + Q R,1 (ξ), Y R,2 (ξ) = |Ω R | -1 Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,2 (ξ)χ ΩR } + Q R,2 (ξ),
where Q R,1 (ξ) and Q R,2 (ξ) consist of seven and three terms respectively. Let 0 < α < 1 and η > 0 be fixed. Let us show that, under the conditions of Lemma 4.2, the quantities |Q R,j (ξ)|, j = 1, 2 obey an estimate of type (4.6). To do that, take one generical term from Q R,1 (ξ) and one from Q R,2 (ξ):

q R,1 (ξ) := 1 |Ω R | ΩR dx R 3 dz 1 R 3 dz 2 (R R (ξ))(x, z 1 ) × a(z 1 -z 2 ) • ∇ z1 (R R (ξ))(z 1 , z 2 )a(z 2 -x) • ∇ z2 (H R -ξ) -1 W R (ξ) (z 2 , x), q R,2 (ξ) := - 1 |Ω R | ΩR dx R 3 dz (R R (ξ))(x, z) 1 2 a 2 (z -x) (H R -ξ) -1 W R (ξ) (z, x).
We need the following estimates. From (2.16)-(2.17 By replacing R R (ξ) with the r.h.s. of (2.8) in Y R,2 (ξ), then we have:

Y R,2 (ξ) = |Ω R | -1 Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,2 (ξ)χ ΩR } + Q R,2 (ξ),
where Q R,2 (ξ) consists of three terms. Let 0 < α < 1 and η > 0 be fixed. Let us show that, under the conditions of Lemma 4.3, the quantity |Q R,2 (ξ)| obeys an estimate of type (4.7). To do so, take a generical term from Q R,2 (ξ):

q R,2 (ξ) := - 1 |Ω R | ΩR dx R 3 dz (R R (ξ))(x, z) 1 2 a 2 (z -x)(W R (ξ))(z, x).
Introduce on L 2 (R 3 ) the operators Z R (ξ) and T

Ξ (ξ), with Ξ := R or P and R ≥ R 0 , generated via their kernel respectively defined by: ∀(x, y) ∈ R 6 , Z R (x, y; ξ) := 1 2 a 2 (xy)(W R (ξ))(x, y),

∀(x, y) ∈ R 6 \ D, T (1)
Ξ (x, y; ξ) := a(xy)(H Ξξ) -1 (x, y).

Ξ (ξ) are bounded on L 2 (R 3 ). We stress the point that Z R (ξ) can be rewritten as: Z R (ξ) = ĝR T R,2 (ξ) VR (H Pξ) -1 + (H Rξ) -1 VR T P,2 (ξ) + T

R (ξ) VR T

(1)

P (ξ) (1 -g R ),
where we used that a 2 (xy) = {a(xz) + a(zy)} 2 for any x, y, z ∈ R 3 . Now, it comes down to prove that there exists a ϑ > 0 and a polynomial p(• ) s.t. ∀R ≥ max{R 0 , R 1 } and ∀ξ ∈ C obeying dist(ξ, σ(H R ) ∩ σ(H P )) ≥ η:

|Ω R | -1 χ ΩR R R (ξ)Z R (ξ)χ ΩR I1 ≤ p(|ξ|)e -ϑ ξ R α . (4.35)
To do so, it is enough to use these estimates which hold ∀R ≥ max{R 0 , R 1 }:

max χ ΩR ĝR (H Ξ -ξ) -1 I2 , χ ΩR ĝR (H Ξ -ξ) -1 I2 ≤ p(|ξ|) |Ω R |, (4.36)
max VR (H Ξξ) -1 (1g R )χ ΩR I2 , VR T where Q R,1 (ξ) consists of seven terms. Let 0 < α < 1 and η > 0 be fixed. Let us show that, under the conditions of Lemma 4.4, the quantity |Q R,1 (ξ)| obeys an estimate of type (4.8). To do so, take a generical term from Q R,1 (ξ):

q R,1 (ξ) := 1 |Ω R | ΩR dx R 3 dz 1 R 3 dz 2 (R R (ξ))(x, z 1 ) × a(z 1 -z 2 ) • ∇ z1 (R R (ξ))(z 1 , z 2 )a(z 2 -x) • ∇ z2 (W R (ξ))(z 2 , x).
Let us note that from (2.9) and (2.10):

∇R R (ξ) = (∇ĝ R ) (H Pξ) -1 + ĝR ∇(H Pξ) -1 g R + ∇ ĝR (H Pξ) -1 + ĝR ∇(H Pξ) -1 (1g R ),

∇W R (ξ) = ∇ ĝR (H Rξ) -1 VR (H Pξ) -1 + ĝR ∇(H Rξ) -1 VR (H Pξ) -1 (1g R ).

By using the definitions (4.34), the quantity q R,1 (ξ) can be rewritten as:

q R,1 (ξ) = q

(1) R,1 (ξ) + q

(2) R,1 (ξ), where:

q (1) R,1 (ξ) := 1 |Ω R | ΩR dx R 3 dz 1 R 3 dz 2 (R R (ξ))(x, z 1 ) (∇ĝ R )(z 1 )T (1)
P (z 1 , z 2 ; ξ)g R (z 2 ) + ∇ ĝR (z 1 )T

(1)

P (z 1 , z 2 ; ξ)(1 -g R )(z 2 ) a(z 2 -x) • ∇ z2 (W R (ξ))(z 2 , x), q (2)
R,1 (ξ) := -

i |Ω R |
Tr L 2 (R 3 ) χ ΩR R R (ξ) ĝR T P,1 (ξ)g R + ĝR T P,1 (ξ)(1g R )

× ∇ ĝR T

(1) R (ξ) VR (H Pξ) -1 + (H Rξ) -1 VR T

(1)

P (ξ)
+ ĝR -iT R,1 (ξ) VR (H Pξ) -1 + ∇(H Rξ) -1 VR T

(1)

P (ξ) (1 -g R )χ ΩR .
Let us first estimate |q

Y R,1 (ξ) := |Ω R | -1 Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,1 (ξ)T R,1 (ξ)χ ΩR } , Y R,2 (ξ) := |Ω R | -1 Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,2 (ξ)χ ΩR } .
By replacing R R (ξ) with the r.h.s. of (2.12) in Y R,j (ξ), j = 1, 2 then: Y R,1 (ξ) = |Ω R | -1 Tr L 2 (R 3 ) χ ΩR (H Pξ) -1 T P,1 (ξ)T P,1 (ξ)χ ΩR + Q R,1 (ξ), Y R,2 (ξ) = |Ω R | -1 Tr L 2 (R 3 ) χ ΩR (H Pξ) -1 T P,2 (ξ)χ ΩR + Q R,2 (ξ), where Q R,1 (ξ) and Q R,2 (ξ) consist of seven and three terms respectively. Let 0 < α < 1 and η > 0 be fixed. Let us show that, under the conditions of Lemma 4.5, the quantities |Q R,j (ξ)|, j = 1, 2 obey an estimate of type (4.9). To do so, take one generical term from Q R,1 (ξ) and one from Q R,2 (ξ):

q R,2 (ξ) := 1 |Ω R | ΩR dx R 3 dz (H P -ξ) -1 (x, z) 1 2 a 2 (z -x)(W R (ξ))(z, x), q R,1 (ξ) := - 1 |Ω R | ΩR dx R 3 dz 1 R 3 dz 2 (H P -ξ) -1 (x, z 1 )a(z 1 -z 2 ) × ∇ z1 (H P -ξ) -1 (z 1 , z 2 )a(z 2 -x) • ∇ z2 (W R (ξ))(z 2 , x).
From (2.16)-(2.17), (2.22) combined with [12, Lem. A.2 (ii)], then one obtains ∀R ≥ R 0 , |q R,j (ξ)| ≤ p(|ξ|)e -ϑ ξ R α j = 1, 2, for another constant ϑ > 0 and polynomial p(• ) both R-independent. The other terms coming from Q R,j (ξ), j = 1, 2 can be treated by using similar arguments.

5 Acknowledgments.

The author was partially supported by the Lundbeck Foundation, and the European Research Council under the European Community's Seventh Framework Program (FP7/2007-2013)/ERC grant agreement 202859. B.S. warmly thanks Horia D. Cornean and Søren Fournais for many fruitful and stimulating discussions. As well, the author thanks the referee for a number of relevant comments and helpful suggestions that helped improve the present manuscript.

1

  Introduction & the main result. 1.1 An historical review. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 What are the motivations of this paper? . . . . . . . . . . . . . . . . . . . . . . . . 1.3 The setting and the main result. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Outline of the proof of Theorem 1.1. . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Discussion on the assumptions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 An open problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 The content of the paper. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 An approximation of the resolvent via a geometric perturbation theory. 3 Proof of Theorem 1.1. 3.1 Proof of (i). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 The location of the Fermi energy. . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Isolating the main R-dependent contribution at zero-temperature. . . . . . 3.1.3 Isolating the main R-dependent contribution at zero-temperature -Continuation and end. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Proof of (ii). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Proof of (iii). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Introduction & the main result.1.1 An historical review.

  (b)} j≥1 the set of eigenvalues of H R,L (b) counting multiplicities and in increasing order. As well, ∀R > 0 and ∀b ∈ R, denote by N R,L (E), E ∈ R the number of eigenvalues (counting multiplicities) of H R,L (b) in the interval (-∞, E).

  g. [12, Thms. 1.1 & 1.2] and [11, Sec. 3.1]. Denoting ∀β > 0, ∀z > 0, ∀b ∈ R and ∀R > 0 the bulk pressure by P R (β, z, b) := lim L↑∞ P R,L (β, z, b), then under our conditions, one has the following pointwise convergences:

  t. ξ in (3.13)-(3.14) depend on β, see (3.15)-(3.4). When β ↑ ∞, they reduce to a half-line on the real axis. Here the isolating situation in (1.41) plays a crucial role allowing to dodge this difficulty. In fact, our assumptions have been chosen to force this situation. Indeed, for R sufficiently large, E R (ρ 0 (R)) lies in the interior of an open set I n0 containing λn 0 +λn 0 +1 2

. 7 )

 7 In view of (3.7), it remains to use[START_REF] Briet | A rigorous proof of the Landau-Peierls formula and much more[END_REF] Thm. 1.1] what leads to (3.2).

Remark 3 . 5 dim E l . 3 . 1 . 2

 35312 g. [27, Thm. 2]: |E R,l (k)| -|λ l | = O e - √ |λ l |R R , l = 1, . . . , τ. From the foregoing, one can infer a sufficient condition which leads to the insulating situation when allowing the eigenvalues of H P in (-∞, 0) to be finitely degenerate. From Lemma 3.4, then for R sufficiently large, (3.7) holds ∀E belonging to a spectral gap of H R provided that: ∃κ ∈ {1, . . . , ν} s.t. n 0 = κ l=1 Isolating the main R-dependent contribution at zero-temperature.

R

  (β, ρ 0 (R), b = 0) ∈ R be the unique solution of the equation ρ R (β, e βµ , b = 0) = ρ 0 (R). Then, from (1.24) together with (3.13)-(3.14), one has ∀β > 0 and ∀R > 0:

Proposition 3 . 6

 36 Let I ς , ς ∈ {1, . . . , τ } be an open interval s.t. I ς (λ ς , λ ς+1 ) and λς +λς+1 2 ∈ I ς if ς < τ ; otherwise I τ (λ τ , 0) and λτ 2 ∈ I τ . Then, for any R > 0 and for any compact subset K ⊂ I ς :

Proposition 3 . 7

 37 Let I ς , ς ∈ {1, . . . , τ } be an open interval as in Proposition 3.6. Then, for R sufficiently large, lim β↑∞ ∆ R (β, µ) exists uniformly on compact subsets K ⊂ I ς . Denote it by ∆ R (µ). Furthermore, ∀0 < α < 1 there exist two constants c, C > 0 s.t. ∀µ ∈ I ς and for R sufficiently large: |∆ R (µ)| ≤ C(1 + |µ|)e -cR α .(3.17)

0) and λτ 2 ∈

 2 I n0 . From Lemma 3.4 and Proposition 3.1, for R sufficiently large, I n0 ∩ σ(H R ) = ∅ and the Fermi energy E R (ρ 0 (R)) := lim β↑∞ µ (0) R (β, ρ 0 (R), b = 0) ∈ I n0 respectively. Then, (3.18) follows from Propositions 3.6 and 3.7 together.

  t. ∀|b| ≤ b 0 , ∪ τ l=1 γ l ∈ ̺(H P (b)) (the resolvent set). Moreover, by [5, Thm. 6.1], the eigenvalues λ l , l = 1, . . . , τ are stable under the perturbation H P (b) -H P for b sufficiently small. Due to the assumption (A nd ), this means that there exists b 1 > 0 s.t. ∀|b| ≤ b 1 , H P (b) has exactly one and only one eigenvalue λ l (b) nearby λ l , l = 1, . . . , τ . Furthermore, by the asymptotic perturbation theory, ∀l ∈ {1, . . . , τ } λ l (b) can be written in the first order in b as λ l (b) = λ l + be l + o(b), see e.g. [35, Sec. VIII, Thm. 2.6]. Actually, the λ l (b)'s can be written in terms of an asymptotic power series in b, see e.g. [7, Thm. 1.2]. From the foregoing, then there exists 0 < b ≤ min{b 0 , b 1 } s.t. ∀|b| ≤ b and ∀l ∈ {1, . . . , τ }, λ l (b) lies inside the closed contour γ l surrounding λ l . Under these conditions, we denote hereafter by Π l (b) the orthogonal projection onto the eigenvector corresponding to the eigenvalue λ l (b).

  T P,1 (ξ) -T P,2 (ξ)] . (3.25) Since the λ l (• )'s are twice differentiable w.r.t. b in a neighborhood of b = 0, then (1.31) follows directly from (3.23) along with Proposition 3.11.

  6 . This follows by applying j-times [12, Lem. A.1] together with (3.27)-(3.28). Furthermore, from (3.26), (3.27)-(3.28) along with [12, Lem. A.2 (ii)], there exists a b-independent polynomial p(• ) s.t. for any k ∈ {1, 2}, m ∈ {0, 1} and ∀(x, y) ∈ R 6 :

Now, let us turn to the proof of Proposition 3 . 11 .

 311 It requires the two following technical results whose proofs lie in Appendix, see Sec. 4.3.

  t. H P (b) has exactly one and only one eigenvalue λ l (b) near λ l . Denote by Φ l (• ; b) the associated (normalized) eigenfunction. Then, there exists b > 0 and two constants c, C > 0 s.t. for any k ∈ {1, 2, 3}:

Corollary 3 . 15

 315 Let {λ l } τ l=1 be the set of eigenvalues of H P in (-∞, 0). Let b > 0 s.t. ∀|b| ≤ b and ∀l ∈ {1, . . . , τ }: (i). H P (b) has exactly one and only one eigenvalue λ l (b) located nearby λ l . (ii). The (normalized) eigenfunction associated with λ l (b) obeys estimate (3.32).

6 . 1 ]

 61 , there exists b 1 > 0 s.t. ∀|b| ≤ b 1 and ∀l ∈ {1, . . . , τ }, λ l is stable under the perturbation W (b) := H P (b) -H P = -ba • (-i∇) + 1 2 b 2 a 2 . This means that ∀|b| ≤ b 1 , H P (b) has exactly one and only one eigenvalue λ l (b) nearby λ l which reduces to λ l in the limit b → 0. For such b's and any l ∈ {1, . . . , τ }, denote by Φ l (b) the normalized eigenfunction associated with λ l (b), and by Π l (b) = |Φ l (b) Φ l (b)| the orthogonal projection onto the eigenvector Φ l (b). We define Π ⊥ l (b) := ½ -Π l (b).

. 37 )

 37 Since Π l (b) converges to Π l in norm as b → 0, see e.g.[START_REF] Kato | Perturbation Theory for Linear Operators[END_REF] Sec. VIII.2], the square-root in (3.37) is well-defined by a binomial series. U (b) is a unitary operator, and it intertwines both projections:Π l (b) = U (b)Π l U * (b). (3.38) From (3.38), one gets Π ⊥ l (b) = U (b)Π ⊥ l U * (b) what leads to U * (b)Π ⊥ l (b) = Π ⊥ l U * (b).As a result, ∀|b| ≤ b and small enough, and for any ξ ∈ K:

  Υ(ξ, b) is an operator s.t. Υ(ξ, b) = O(|b|) uniformly in ξ. It follows that Π ⊥ l (H P (b)-ξ)Π ⊥ l is invertible ∀|b| ≤ b and small enough, and ∀ξ ∈ K.

4. 1 7 . 2 ∈ 3 , λ ς + 2λ ς+1 3 ,

 17233 Proof of Propositions 3.6 and 3.We start by introducing some notation. Recall that under the assumptions (A r )-(A m )-(A nd ), {λ l } τ l=1 with τ ∈ N * denotes the set of eigenvalues of H P in (-∞, 0) counting in increasing order. For simplicity's sake, we set λ τ +1 := 0. Let I ς , ς ∈ {1, . . . , τ } be an open interval s.t. I ς (λ ς , λ ς+1 ) and λς +λς+1 I ς . Without loss of generality, we give an explicit form to I ς : I ς := 2λ ς + λ ς+1 with the convention λ τ +1 := 0. (4.1) Introduce the following decomposition of the contour C (P ) β defined in (3.15):

Θ

  R (ξ) := Tr L 2 (R 3 ) χ ΩR (H Pξ) -1 [T P,1 (ξ)T P,1 (ξ) -T P,2 (ξ)] χ ΩR , and (4.3) holds ∀R > 0. Since [ω - ς , ω + ς ] ∩ σ(H P ) = ∅, then ∀µ ∈ I ς the closed subset surrounding by γ (2)

Remark 4 . 1

 41 We emphasize that the deformation of Γς,β and γ[START_REF] Angelescu | On the Landau diamagnetism[END_REF]

  x)| ≤ Ce -c|x| ∀x ∈ R 3 , for another c > 0 and C = C(θ, ς) > 0. A similar estimate still holds true for |M (j),⊥,ℵ θ,w (x, x)|. This comes from the previous arguments, the upper bound (4.17) and the rough estimate |a(xy)| ≤ (|x| + |y|), combined with:

4. 3 2 - 3 + 2 k=1( 2 k=1(

 32322 Proof of Lemma 3.13.Introduce ∀(b, b 0 ) ∈ R 2 and ∀ξ ∈ ̺(H P (b 0 )), the operators RP (b, b 0 , ξ) and TP,j (b, b 0 , ξ), j = 1, 2 on L 2 (R 3 ) via their kernel respectively defined on R 6 \ D by: RP (x, y; b, b 0 , ξ) := e iδbφ(x,y) (H P (b 0 )ξ) -1 (x, y), δb := bb 0 , (4.21) TP,j (x, y; b, b 0 , ξ) := e iδbφ(x,y) T P,j (x, y; b 0 , ξ),(4.22) where φ is the magnetic phase defined in(3.26). Also, set:TP (b, b 0 , ξ) := δb TP,1 (b, b 0 , ξ) + (δb) 2 TP,2 (b, b 0 , ξ). (4.23)Except for a gauge phase factor, the kernel of RP (b, b 0 , ξ) and TP,j (b, b 0 , ξ), j = 1, 2 is the same as the one of (H P (b 0 )ξ) -1 and T P,j (b 0 , ξ) respectively. Therefore, ∀η > 0 and ∀ξ ∈ C satisfying dist(ξ, σ(HP (b 0 ))) ≥ η, then ∀b ∈ R RP (b, b 0 , ξ)and TP,j (b, b 0 , ξ) are bounded with operator norm obeying (3.29) (with b 0 instead of b). Under the same conditions, introduce on L 2 (R 3 ): T(1) P (b, b 0 , ξ) := -RP (b, b 0 , ξ) TP,1 (b, b 0 , ξ), (4.24) T(2) P (b, b 0 , ξ) := RP (b, b 0 , ξ) TP,1 (b, b 0 , ξ) TP,2 (b, b 0 , ξ) , (4.25) as well as, with the additional condition ξ ∈ ̺(H P (b)) ∩ ̺(H P (b 0 )): T(3) P (b, b 0 , ξ) := -(H P (b)ξ) -1 TP (b, b 0 , ξ) RP (b, b 0 , ξ) 2 r=1 TP,ir (b, b 0 , ξ). (4.26) Now, we turn to the actual proof. Let K ⊂ (̺(H P ) ∩ {ζ ∈ C : ℜζ < 0}) be a compact subset. From [7, Thm. 1.1], then there exists b K > 0 s.t. ∀|b| ≤ b K , K ⊂ (̺(H P (b)) ∩ {ζ ∈ C : ℜζ < 0}). From now on, let b 0 ∈ (-b K , b K ) be fixed. The starting point of the gauge invariant magnetic perturbation theory is the following identity which holds ∀ξ ∈ K in the bounded operators sense on L 2 (R 3 ), see [18, Proof of Prop. 3.2] and also [12, Lem. 3.2]: ∀|b| ≤ b K , (H P (b)ξ) -1 = RP (b, b 0 , ξ) -(H P (b)ξ) -1 TP (b, b 0 , ξ). (4.27) It results from (4.27) that, for b sufficiently close to b 0 , the operator norm of (H P (b)ξ) -1 -RP (b, b 0 , ξ) behaves like O(|δb|). This comes from (4.23), the definitions (4.21)-(4.22) and the estimates (3.27)-(3.28) which yield: ∀|b| ≤ b K , max sup ξ∈K RP (b, b 0 , ξ) , sup ξ∈K TP,j (b, b 0 , ξ) ≤ C, j = 1, 2, for some constant C = C(|b 0 |, K) > 0. Next, by iterating twice (4.27) and in view of (4.24)-(4.26), one arrives ∀|b| ≤ b K and ∀ξ ∈ K on L 2 (R 3 ) at: (H P (b)ξ) -1 = RP (b, b 0 , ξ) + δb) k T(k) P (b, b 0 , ξ) + T(3) P (b, b 0 , ξ). (4.28) In terms of corresponding integral kernels, one has on R 6 \ D: (H P (b)ξ) -1 (x, y) = RP (x, y; b, b 0 , ξ) + δb) k T(k) P (x, y; b, b 0 , ξ) + T(3) P (x, y; b, b 0 , ξ), (4.29) where, for any k ∈ {1, 2}, ∀(x, y) ∈ R 6 and ∀|b| ≤ b K : T(k) P (x, y; b, b 0 , ξ) := k j=1

2 k=1( 2 k=0( 2 k=1(

 222 P (• , • ; b, b 0 , ξ) in (4.29) stands for the kernel of T(3) P (b, b 0 , ξ). Let us note that, in view of (4.26) along with (4.21)-(4.22) and the estimates (3.27)-(3.28), the kernel T(3) P (• , • ; b, b 0 , ξ) behaves like O(|δb| 3 ) uniformly in ξ ∈ K. Afterwards, we remove the b-dependence in the first two terms of the r.h.s. of (4.29). To achieve that, we expand in Taylor power series the exponential phase factor appearing in (4.21) and (4.30) up to the second order in δb. In this way, one obtains ∀|b| ≤ b K and ∀(x, y) ∈ R 6 \ D: RP (x, y; b, b 0 , ξ) + δb) k T(k) P (x, y; b, b 0 , ξ) = δb) k (iφ(x, y)) k k! (H P (b 0 )ξ) -1 (x, y) + (x, y; b 0 , ξ) + T (4) P (x, y; b, b 0 , ξ), (4.31)

( 4 )

 4 P (• , • ; b, b 0 , ξ) satisfies the property that its first two derivatives at b 0 are identically zero. From (4.31), it holds onR 6 \ D for b ∈ [-b K , b K ] sufficiently close to b 0 : (H P (b)ξ) -1 (x, y) -(H P (b 0 )ξ) -1 (x, y) = δb iφ(x, y)(H P (b 0 )ξ) -1 (x, y) -(H P (b 0 )ξ) -1 T P,1 (b 0 , ξ) (x, y) + o (δb) . It follows that the map b → (H P (b)ξ) -1 (x, y) is differentiableat b 0 , and: ∂ ∂b (H P (b)ξ) -1 (x, y) b=b0 := iφ(x, y)(H P (b 0 )ξ) -1 (x, y) -(H P (b 0 )ξ) -1 T P,1 (b 0 , ξ) (x, y).

3 dz

 3 ),(2.22) combined with[START_REF] Briet | A rigorous approach to the magnetic response in disordered systems[END_REF] Lem. A.2], then under the conditions of Lemma 2.2 (i), one has ∀R ≥ R 0 :R 3 dz (H Rξ) -1 (x, z)(W R (ξ))(z, y) ≤ p(|ξ|)e -ϑ ξ R α e -ϑ ξ |x-y| ,(4.32)R ∇ x (H Rξ) -1 (x, z)(W R (ξ))(z, y) ≤ p(|ξ|)e -ϑ ξ R α e -ϑ ξ |x-y| |x -y| , (4.33) for another ϑ > 0 and polynomial p(• ) both R-independent. From (2.18)-(2.19), (4.32)-(4.33) combined with [12, Lem. A.2 (ii)], |q R,j (ξ)| ≤ p(|ξ|)e -ϑ ξ R α for another ϑ > 0 and polynomial p(• ) both R-independent. The other terms coming from Q R,j (ξ), j = 1, 2 can be treated by using similar arguments. Proof of Lemma 4.3. Let us denote: Y R,2 (ξ) := |Ω R | -1 Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,2 (ξ)χ ΩR } .

(4. 34 )

 34 Due to the estimates (2.20) and (2.16), the operators Z R (ξ) and T

Ξ

  (ξ)(1g R )χ ΩR I2 , VR T Ξ,2 (ξ)(1g R )χ ΩR I2 ≤ p(|ξ|) |Ω R |e -ϑ ξ R α , (4.37)with Ξ = R or P , and for another R-independent polynomial p(• ). Here, we used the definitions (4.34), (3.11)-(1.28) with (2.16) and Lemma 2.2 (iii). The two other terms of Q R,2 (ξ) can be treated by using similar arguments.Proof of Lemma 4.4. Let us denote:Y R,1 (ξ) := |Ω R | -1 Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,1 (ξ)T R,1 (ξ)χ ΩR } .By replacing R R (ξ) with the r.h.s. of (2.8) in Y R,1 (ξ), then we have:Y R,1 (ξ) = |Ω R | -1 Tr L 2 (R 3 ) {χ ΩR R R (ξ)T R,1 (ξ)T R,1 (ξ)χ ΩR } + Q R,1 (ξ),

R, 1

 1 (ξ)|. By mimicking the proof of (2.23), then under the same conditions, there exists another ϑ > 0 and polynomial p(• ) s.t. ∀R ≥ R 0 : max (∇ĝ R ) (x)T (1) P (x, y; ξ)g R (y) , ∇ ĝR (x)T (1) P (x, y; ξ)(1g R )(y) ≤ p(|ξ|)e -ϑ ξ R α e -ϑ ξ |x-y| .From (2.18), (2.21) and the above estimate, one concludes by[12, Lem. A.2] that |q (1) R,1 (ξ)| ≤ p(|ξ|)e -ϑ ξ R α ∀R ≥ R 0 ,for another constant ϑ > 0 and polynomial p(• ) both R-independent. Next, let us estimate |q

( 2 )

 2 R,1 (ξ)|. By mimicking the arguments leading to (4.35), the trace norm of the operator inside the trace in the expression of q

( 2 )

 2 R,1 (ξ) is bounded above ∀R ≥ max{R 0 , R 1 } by polynomial × |Ω R |e -ϑ ξ R α dueto the H-S norms in (4.36)-(4.37) and the operator norms in (2.24)-(2.25). Then ∀R ≥ max{R 0 , R 1 }, |q

( 2 )

 2 R,1 (ξ)| ≤ p(|ξ|)e -ϑ ξ R α for another constant ϑ > 0 and polynomial p(• ) both R-independent. The other terms coming from Q R,1 (ξ) can be treated by similar arguments.Proof of Lemma 4.5. Let us denote:

  R,l 's within the Bloch-Floquet theory). If max E R,l < min E R,l+1 for some l ≥ 1 then a spectral gap occur. Under the assumption (A r ), we know that the Fermi energy always exists, see[START_REF] Briet | A rigorous proof of the Landau-Peierls formula and much more[END_REF] Thm. 1.1]. Moreover, only two situations can occur: either the Fermi energy lies in the middle of a spectral gap (this corresponds to the semiconducting/isolating situation if max E R,s < min E R,s+1 , the semimetal situation if max E R,s = min E R,s+1

l+1 (we refer to Sec. 3.1.1 for the definition of the E

  .[START_REF] Nenciu | On asymptotic perturbation theory for quantum mechanics: Almost invariant subspaces and gauge invariant magnetic perturbation theory[END_REF] The proof leans on two ingredients. The first one concerns the location of the negative Bloch bands of H R at negative eigenvalues of H P for large R, see Lemma 3.4. Since the λ l 's (l = 1, . . . , τ ) are simple, then for R sufficiently large, the Bloch bands E R,l (l = 1, . . . , τ ) are simple, isolated from each other and from the rest of the spectrum. Moreover, ∀l ∈ {1, . . . , τ } E R,l is localized in a neighborhood of λ

l , and reduces to λ l when R ↑ ∞. The second one is a relation between the I.D.S. of H R and the zero-temperature limit of the grand-canonical density expressed in terms of the eigenvalues of the Bloch Hamiltonian, see Lemma 3.3. Then, it remains to use the criterion in

[START_REF] Briet | A rigorous proof of the Landau-Peierls formula and much more[END_REF] Thm 1.1]

. Let us get back to the zero-temperature limit. Performing it needs special attention. In (3.13) and (3.14), the dependence in β is contained in

  .17) Proof. See [47, Thm. B.7.2] and [12, Lem. 2.4] respectively.

  , (2.21) follows from (2.16)-(2.17) combined with [12, Lem. A.2]. Next, we prove (2.22). The properties (2.3)-(2.5) are the key-ingredient to make appear the exponential decay in R α . Under the conditions of the lemma:

  Remark 2.3 From Lemma 2.1 and Lemma 2.2 (i)-(ii), together with the Shur-Holmgren criterion, one has ∀R ≥ R 0 :

it is enough to use (2.23), (2.17) along with [12, Eq. (A.12)]. Finally (iii) follows from (2.16) with (2.14) assuring that dist(Supp(u), Supp(χ ΩR ĝR )) ≥ R α .

  with (4.19) leading to the estimate |z 1 |e - √ |λς ||z1| ≤ |λ ς | -1 2 ∀z 1 ∈ R 3 , then by [12, Lem. A.2 (ii)], there exist another constant C > 0 s.t. for any 1 ≤ l ≤ ς: