Philippe Briet
email: briet@univ-tln.fr

Horia D Cornean
email: cornean@math.aau.dk

Baptiste Savoie
email: bap-tiste.savoie@gmail.com

W Kohn

A rigorous proof of the Landau-Peierls formula and much more

We present a rigorous mathematical treatment of the zero-field orbital magnetic susceptibility of a non-interacting Bloch electron gas, at fixed temperature and density, for both metals and semiconductors/insulators. In particular, we obtain the Landau-Peierls formula in the low temperature and density limit as conjectured by T.

Introduction and the main results

Understanding the zero-field magnetic susceptibility of a Bloch electron gas is one of the oldest problems in quantum statistical mechanics.

The story began in 1930 with a paper by L. Landau [START_REF] Landau | Diamagnetismus der Metalle[END_REF], in which he computed the diamagnetic susceptibility of a free degenerate gas. (Note that the rigorous proof of Landau's formula for free electrons was given by Angelescu et al. [START_REF] Angelescu | On the Landau Diamagnetism[END_REF] and came as late as 1975.) For Bloch electrons (which are subjected to a periodic background electric potential), the problem is much harder and -to our best knowledge-it has not been solved yet in its full generality.

The first important contribution to the periodic problem came in 1933, when R. Peierls [START_REF] Peierls | Zur Theorie des Diamagnetismus von Leitungselektronen[END_REF] introduced his celebrated Peierls substitution and constructed an effective band Hamiltonian which permitted to reduce the problem to free electrons. Needless to say that working with only one energy band instead of the full magnetic Schrödinger operator is an important simplification. Under the tight-binding approximation he claimed that the dominant contribution to the zero-field orbital susceptibility of a Bloch electron gas in metals (at zero temperature) is purely diamagnetic and is given by the so-called Landau-Peierls formula which consists of replacing in the Landau formula the mass by the effective mass of the electron. He showed as well the existence of another contribution which has no simple interpretation and whose magnitude and sign are uncertain.

In 1953, E.N. Adams [START_REF] Adams | Magnetic susceptibility of a diamagnetic electron gas -The role of small effective electron mass[END_REF] claimed that the Landau-Peierls susceptibility is not always the dominant contribution to the zero-field orbital susceptibility. By considering the case of 'simple metals' (for which the tight-binding approximation is not appropriate), he showed that there exists others contributions (certain have even positive sign!) coming from the bands not containing the Fermi energy. Besides in special cases these contributions are of the same order of magnitude than the Landau-Peierls formula. However any general formulation of these contributions were stated.

In 1957, T. Kjeldaas and W. Kohn [START_REF] Kjeldaas | Theory of the Diamagnetism of Bloch Electrons[END_REF] were probably the first ones who suggested that for 'simple metals' the Landau-Peierls approximation is only valid in the limit of weak density of electrons, moreover, the Landau-Peierls formula (see below (1.17) and (1.18)) has to be corrected with some higher order terms in the particle density, and these terms must come from the bands not containing the Fermi energy.

These three papers generated a lot of activity, where the goal was to write down an exact expression for the zero-field magnetic susceptibility of a Bloch electron gas in metals at zero temperature. In what follows we comment on some of the most important works.

The first attempt to address the full quantum mechanical problem -even though the carriers were boltzons and not fermions-was made by J.E. Hebborn and E.H. Sondheimer [START_REF] Hebborn | Diamagnetism of Conduction Electrons in Metals[END_REF][START_REF] Hebborn | The Diamagnetism of Conduction Electrons in Metals[END_REF] in 1960. Unlike the previous authors, they developed a magnetic perturbation theory for the trace per unit volume defining the pressure. The biggest problem of their formalism is that they assumed that all Bloch energy bands are not overlapping (this is generically false; for a proof of the Bethe-Sommerfeld conjecture in dimension 3 see e.g. [START_REF] Helffer | Asymptotics of the density of states for the Schrödinger operator with periodic potential[END_REF]), and that the Bloch basis is smooth in the quasi-momentum variables. This assumption can fail at the points where the energy bands cross each other. Not to mention that no convergence issues were addressed in any way.

In 1962, L.M. Roth [START_REF] Roth | Theory of Bloch Electrons in a Magnetic Field[END_REF] developed a sort of magnetic pseudodifferential calculus starting from the ideas of Peierls, Kjeldaas and Kohn. She used this formalism in order to compute local traces and magnetic expansions. Similar results are obtain by E.I. Blunt [START_REF] Blount | Bloch Electrons in a Magnetic Field[END_REF]. Their formal computations can most probably be made rigorous in the case of simple bands.

In 1964, Hebborn et al. [START_REF] Hebborn | The Orbital Diamagnetic Susceptibility of Bloch Electrons[END_REF] simplified the formalism developed in [START_REF] Hebborn | The Diamagnetism of Conduction Electrons in Metals[END_REF] and gave for the first time a formula for the zero-field susceptibility of a boltzon gas. Even though the proofs lack any formal rigor, we believe that their derivation could be made rigorous for systems where the Bloch bands do not overlap. But this is generically not the case.

The same year, G.H. Wannier and U.N. Upadhyaya [START_REF] Wannier | Zero-Field Susceptibility of Bloch Electrons[END_REF] go back to the method advocated by Peierls, and replace the true magnetic Schrödinger operator with a (possibly infinite) number of bands modified with the Peierls phase factor. They claim that their result is equivalent with that one of Hebborn and Sondheimer [START_REF] Hebborn | The Diamagnetism of Conduction Electrons in Metals[END_REF], but no details are given. Anyhow, the result uses in an essential way the non-overlapping of Bloch bands. At the same time, L. Glasser [START_REF] Glasser | Magnetic properties of nearly free electrons, nonoscillatory magnetic susceptibility[END_REF] gave an expression of the bulk zero-field susceptibility in terms of effective mass by the usual nearly free electron approximation.

In 1969, P.K. Misra and L.M. Roth [START_REF] Misra | Theory of Diamagnetic Susceptibility of Metals[END_REF] combined the method of [START_REF] Roth | Theory of Bloch Electrons in a Magnetic Field[END_REF] with the ideas of Wannier in order to include the core electrons in the computation.

In 1972, P.K. Misra and L. Kleinman [START_REF] Misra | Theory of the Diamagnetic Susceptibility of Bloch Electrons[END_REF] had the very nice idea of using sum-rules in order to replace derivatives with respect to the quasi-momentum variables, with matrix elements of the "true" momentum operator. They manage in this way to rewrite the formulas previously derived by Misra and Roth (which only made sense for non-overlapping bands) in a form which might also hold for overlapping bands.

As we have already mentioned, the first serious mathematical approach on the zero-field susceptibility appeared as late as 1975, due to Angelescu et al. [START_REF] Angelescu | On the Landau Diamagnetism[END_REF]. Then in 1990, B. Helffer and J. Sjöstrand [START_REF] Helffer | On diamagnetism and the de Haas-Van Alphen effect[END_REF] developed for the first time a rigorous theory based on the Peierls substitution and considered the connection with the de Haas-Van Alphen effect. These and many more results were reviewed by G. Nenciu in 1991 [START_REF] Nenciu | Dynamics of band electrons in electric and magnetic fields: rigorous justification of the effective Hamiltonians[END_REF]. A related problem in which the electron gas is confined by a trapping potential was considered by M. Combescure and D. Robert in 2001 [START_REF] Combescure | Rigorous semiclassical results for the magnetic response of an electron gas[END_REF]. They obtained the Landau formula in the limit → 0.

Finally we mention that the magnetic response can be described using the semiclassical theory of the orbital magnetism and the Berry-phase formula, see [START_REF] Resta | Electric polarization and orbital magnetisation: the modern theories[END_REF] for further details. The link between this approach and our work has yet to be clarified.

Our current paper is based on what we call magnetic perturbation theory, as developed by the authors and their collaborators in a series of papers starting with 2000 (see [START_REF] Cornean | On eigenfunction decay for two dimensional magnetic Schrödinger operators[END_REF][START_REF] Cornean | On the magnetization of a charged Bose gas in the canonical ensemble[END_REF][START_REF] Cornean | The Faraday effect revisited: general theory[END_REF][START_REF] Cornean | The Faraday effect revisited: Thermodynamic limit[END_REF][START_REF] Cornean | Faraday effect revisited: sum rules and convergence issues[END_REF][START_REF] Briet | Locating the spectrum for magnetic Schrödinger and Dirac operators[END_REF][START_REF] Briet | Diamagnetic expansions for perfect quantum gases[END_REF][START_REF] Briet | Generalized susceptibilities for a perfect quantum gas. Markov Process[END_REF][START_REF] Briet | Diamagnetic expansions for perfect quantum gases II: uniform bounds[END_REF][START_REF] Briet | Diamagnetism of quantum gases with singular potentials[END_REF] and references therein). The results we obtain in Theorem 1.2 give a complete answer to the problem of zero-field susceptibility. Let us now discuss the setting and properly formulate the mathematical problem.

The setting

Consider a confined quantum gas of charged particles obeying the Fermi-Dirac statistics. The spin is not considered since we are only interested in orbital magnetism. Assume that the gas is subjected to a constant magnetic field and an external periodic electric potential. The interactions between particles are neglected and the gas is at thermal equilibrium.

The gas is trapped in a large cubic box, which is given by Λ

L = -L 2 , L 2
3 , L ≥ 1.
Let us introduce our one-body Hamiltonian. We consider a uniform magnetic field B = (0, 0, B) with B ≥ 0, parallel to the third direction of the canonical basis of R 3 . Let a(x) be the symmetric (transverse) gauge a(x) := 1 2 (-x 2 , x 1 , 0) which generates the magnetic field (0, 0, 1). We consider that the background electric potential V is smooth, i.e. V ∈ C ∞ (R 3 ) is a realvalued function and periodic with respect to a (Bravais) lattice Υ with unit cell Ω. Without loss of generality, we assume that Υ is the cubic lattice Z 3 , thus Ω is the unit cube centered at the origin of coordinates.

When the box is finite i.e. 1 ≤ L < ∞, the dynamics of each particle is determined by a Hamiltonian defined in L 2 (Λ L ) with Dirichlet boundary conditions on ∂Λ L :

H L (ω) = 1 2 -i∇ x -ωa(x) 2 + V L (x) (1.1)
where V L stands for the restriction of V to the box Λ L . Here ω := e c B ∈ R denotes the cyclotron frequency. The operator

H L (ω) is self-adjoint on the domain D(H L (ω) = H 1 0 (Λ L ) ∩ H 2 (Λ L ).
It is well known (see [START_REF] Reed | Methods of Modern Mathematical Physics, II : Fourier Analysis and Self-Adjointness[END_REF]) that H L (ω) is bounded from below and has compact resolvent. This implies that its spectrum is purely discrete with an accumulation point at infinity. We denote the set of eigenvalues (counting multiplicities and in increasing order) by {e j (ω)} j≥1 .

When L = ∞ we denote by H ∞ (ω) the unique self-adjoint extension of the operator

1 2 -i∇ x -ωa(x) 2 + V (x) (1.2)
initially defined on C ∞ 0 (R 3 ). Then H ∞ (ω) is bounded from below and only has essential spectrum (see e.g. [START_REF] Briet | Do Bosons Condense in a Homogeneous Magnetic Field ?[END_REF]). Now let us define some quantum statistical quantities related to the quantum gas introduced above. For the moment we use the grand canonical formalism. The finite volume pressure and density of our quantum gas at inverse temperature β := (k B T ) -1 > 0 (k B stands for the Boltzmann constant), at fugacity z := e βµ > 0 (µ ∈ R stands for the chemical potential) and at cyclotron frequency ω ∈ R are given by (see e.g. [START_REF] Huang | Statistical Mechanics[END_REF]):

P L (β, z, ω) := 1 β|Λ L | Tr L 2 (ΛL) ln 1 + ze -βHL(ω) = 1 β|Λ L | ∞ j=1 ln 1 + ze -βej(ω) (1.3) ρ L (β, z, ω) := βz ∂P L ∂z (β, z, ω) = 1 |Λ L | ∞ j=1 ze -βej(ω) 1 + ze -βej(ω) .
(1.4)

As the semi-group e -βHL(ω) is trace class, the series in (1.3) and (1.4) are absolutely convergent.

Since the function R ∋ ω → P L (β, z, ω) is smooth (see [START_REF] Briet | Diamagnetism of quantum gases with singular potentials[END_REF]), we can define the finite volume orbital susceptibility as the second derivative of the pressure with respect to the intensity B of the magnetic field at B = 0 (see e.g. [START_REF] Angelescu | On the Landau Diamagnetism[END_REF]):

X GC L (β, z) := e c 2 ∂ 2 P L ∂ω 2 (β, z, 0). (1.5)
When Λ L fills the whole space, we proved in [START_REF] Savoie | [END_REF] that the thermodynamic limits of the three grand canonical quantities defined above exist. By denoting P ∞ (β, z, ω) := lim L→∞ P L (β, z, ω), we proved moreover the following pointwise convergence:

ρ ∞ (β, z, ω) := βz ∂P ∞ ∂z (β, z, ω) = lim L→∞ βz ∂P L ∂z (β, z, ω) (1.6) X GC ∞ (β, z) := e c 2 ∂ 2 P ∞ ∂ω 2 (β, z, 0) = lim L→∞ e c 2 ∂ 2 P L ∂ω 2 (β, z, 0) (1.7)
and the limit commutes with the first derivative (resp. the second derivative) of the grand canonical pressure with respect to the fugacity z (resp. to the external magnetic field B). Now assume that our fixed external parameter is the density of particles ρ 0 > 0. We prefer to see ρ ∞ as a function of the chemical potential µ instead of the fugacity z; the density is a strictly increasing function with respect to both µ and z. Denote by µ ∞ (β, ρ 0 ) ∈ R the unique solution of the equation:

ρ 0 = ρ ∞ β, e βµ∞(β,ρ0) , 0 . (1.8)
The bulk orbital susceptibility at β > 0 and fixed density ρ 0 > 0 defined from (1.7) is defined as: ,ρ0) .

X (β, ρ 0 ) := X GC ∞ β, e βµ∞(β
(1.9)

In fact one can also show that

X (β, ρ 0 ) = -e c 2 ∂ 2 f∞ ∂ω 2 (β, ρ 0 , 0) where f ∞ (β, ρ 0 , ω
) is the thermodynamic limit of the reduced free energy defined as the Legendre transform of the thermodynamic limit of the pressure (see e.g. [START_REF] Ruelle | Statistical Mechanics -Rigorous Results[END_REF]). Note that for a perfect quantum gas and in the limit of low temperatures, (1.9) leads to the so-called Landau diamagnetic susceptibility, see e.g. [START_REF] Angelescu | On the Landau Diamagnetism[END_REF].

In order to formulate our main result, we need to introduce some more notation. In the case in which ω = 0, the Floquet theory for periodic operators (see e.g. [START_REF] Berezin | The Schrödinger Equation[END_REF], [START_REF] Kuchment | Floquet Theory for Partial Differential Equations[END_REF] and section 3) allows one to use the band structure of the spectrum of H ∞ (0). Denote by Ω * = 2πΩ the Brillouin zone of the dual lattice Υ * ≡ 2πZ 3 .

If j ≥ 1, the jth Bloch band function is defined by E j := [min k∈Ω * E j (k), max k∈Ω * E j (k)] where {E j (k)} j≥1 is the set of eigenvalues (counting multiplicities and in increasing order) of the fiber Hamiltonian h(k) := 1 2 (-i∇ + k) 2 + V living in L 2 (T 3 ) with T 3 := R 3 /Z 3 the 3-dimensional torus. With this definition, the Bloch energies E j (•) are continuous on the whole of Ω * , but they are differentiable only outside a zero Lebesgue measure subset of Ω * corresponding to cross-points.

In the following we make the assumption that the E j 's are simple eigenvalues for k in a subset of Ω * with full measure. Note that this assumption is not essential for our approach but it simplifies the presentation, see Remark 3 below the Theorem 1.2.

The spectrum of H ∞ (0) is absolutely continuous and given (as a set of points) by σ(H ∞ (0)) = ∞ j=1 E j . Note that the sets E j can overlap each other in many ways, and some of them can even coincide even though they are images of increasingly ordered functions. The energy bands are disjoint unions of E j 's. Moreover, if max E j < min E j+1 for some j ≥ 1 then we have a spectral gap. Since the Bethe-Sommerfeld conjecture holds true under our conditions [START_REF] Helffer | Asymptotics of the density of states for the Schrödinger operator with periodic potential[END_REF], the number of spectral gaps is finite, if not zero.

It remains to introduce the integrated density of states of the operator H ∞ (0). Recall its definition. For any E ∈ R, let N L (E) be the number of eigenvalues of H L (0) not greater than E. The integrated density of states of H ∞ (0) is defined by the limit (see [START_REF] Doi | The Uniqueness of the Integrated Density of States for the Schrödinger Operators with Magnetic Fields[END_REF]):

n ∞ (E) := lim L→∞ N L (E) |Λ L | = lim L→∞ Tr χ (-∞,E] (H L (0)) |Λ L | (1.10)
and n ∞ (•) is a positive continuous and non-decreasing function (see e.g. [START_REF] Berezin | The Schrödinger Equation[END_REF]). In this case one can express n ∞ (E) with the help of the Bloch energies in the following way:

n ∞ (E) = 1 (2π) 3 j≥1 Ω * χ [E0,E] (E j (k)) dk (1.11)
where

χ [E0,E] (•) is the characteristic function of the interval [E 0 , E].
Thus n ∞ is clearly continuous in E due to the continuity of the Bloch bands. Moreover, this function is piecewise constant when E belongs to a spectral gap.

The statements of our main results

The first theorem is not directly related to the magnetic problem, and it deals with the rigorous definition of the Fermi energy for Bloch electrons. Even though these results are part of the 'physics folklore', we have not found a serious mathematical treatment in the literature.

Theorem 1.1. Let ρ 0 > 0 be fixed. If µ ∞ (β, ρ 0 ) is the unique real solution of the equation ρ ∞ β, e βµ , 0 = ρ 0 (see (1.8)), then the limit:

E F (ρ 0 ) := lim β→∞ µ ∞ (β, ρ 0 ) (1.12)
exists and defines an increasing function of ρ 0 called the Fermi energy. There can only occur two cases: SC (semiconductor/insulator/semimetal): Suppose that there exists some

N ∈ N * such that ρ 0 = n ∞ (E) for all E ∈ [max E N , min E N +1 ]. Then: E F (ρ 0 ) = max E N + min E N +1 2 .
(1.13) M (metal): Suppose that there exists a unique solution E M of the equation n ∞ (E M ) = ρ 0 which belongs to (min E N , max E N ) for some (possibly not unique) N . Then :

E F (ρ 0 ) = E M . ( 1
:= inf σ(H ∞ (0)).
(i). Assume that the Fermi energy is in the middle of a non-trivial gap (see (1.13)). Then there exist 2N functions c j (• ), d j (• ), with 1 ≤ j ≤ N , defined on Ω * outside a set of Lebesgue measure zero, such that the integrand in (1.15) can be extended by continuity to the whole of Ω * and:

X SC (ρ 0 ) := lim β→∞ X (β, ρ 0 ) = e c 2 1 2 1 (2π) 3 Ω * dk N j=1 c j (k) + E j (k) -E F (ρ 0 ) d j (k) . (1.15)
(ii). Suppose that there exists a unique N ≥ 1 such that E F (ρ 0 ) ∈ (min E N , max E N ). Assume that the Fermi surface S F := {k ∈ Ω * : E N (k) = E F (ρ 0 )} is smooth and non-degenerate. Then there exist 2N + 1 functions F N (•), c j (•), d j (•) with 1 ≤ j ≤ N , defined on Ω * outside a set of Lebesgue measure zero, in such a way that they are all continuous on S F while the second integrand in (1.16) can be extended by continuity to the whole of Ω * :

X M (ρ 0 ) := lim β→∞ X (β, ρ 0 ) = - e c 2 1
12 1 (2π) 3 (1.16)

SF dσ(k) ∇E N (k) ∂ 2 E N (k) ∂k 2 1 ∂ 2 E N (k) ∂k 2 2 - ∂ 2 E N (k) ∂k 1 ∂k 2 2 -3F N (k) -6 Ω * dk N j=1 χ [E0,EF (ρ0)] E j (k) c j (k) + E j (k) -E F (ρ 0 ) χ [E0,EF (ρ0)] E j (k) d j (k) .
Here χ [E0,EF (ρ0)] (• ) denotes the characteristic function of the interval E 0 ≤ t ≤ E F (ρ 0 ). (iii). Let k F := (6π 2 ρ 0 ) 1 3 be the Fermi wave vector. Then in the limit of small densities, (1.16) gives the Landau-Peierls formula:

X M (ρ 0 ) = - e 2 24π 2 c 2 (m * 1 m * 2 m * 3 ) 1 3 m * 1 m * 2 k F + o(k F ); (1.
X M (ρ 0 ) ∼ - e 2 24π 2 m * c 2 k F when k F → 0. (1.18)
Note that our expression is twice smaller than the one in [START_REF] Peierls | Zur Theorie des Diamagnetismus von Leitungselektronen[END_REF] since we do not take into account the degeneracy related to the spin of the Bloch electrons. Remark 5 . The assumption V ∈ C ∞ (T 3 ) can be relaxed to V ∈ C r (T 3 ) with r ≥ 23. The smoothness of V plays an important role in the absolute convergence of the series defining X (β, ρ 0 ) in Theorem 3.1, before the zero-temperature limit; see [START_REF] Cornean | Faraday effect revisited: sum rules and convergence issues[END_REF] for a detailed discussion on sum rules and local traces for periodic operators. Remark 6 . The role of magnetic perturbation theory (see Section 3) is crucial when one wants to write down a formula for X (β, ρ 0 ) which contains no derivatives with respect to the quasimomentum k. Remember that the Bloch energies ordered in increasing order and their corresponding eigenfunctions are not necessarily differentiable at crossing points. Remark 7 . We do not treat the semi-metal case, in which the Fermi energy equals E F (ρ 0 ) = max E N = min E N +1 for some N ≥ 1 (see (1.13)). This remains as a challenging open problem.

The content of the paper

Let us briefly discuss the content of the rest of this paper:

• In Section 2 we thoroughly analyze the behavior of the chemical potential µ ∞ when the temperature goes to zero defining the Fermi energy. These results are important for our main theorem.

• In Section 3 we give the most important technical result. Applying the magnetic perturbation theory we arrive at a general formula for X (β, ρ 0 ) which contains no derivatives with respect to k. The strategy is somehow similar to the one used in [START_REF] Cornean | The Faraday effect revisited: general theory[END_REF] for the Faraday effect.

• In Section 4 we perform the zero temperature limit and separately analyze the situations in which the Fermi energy is either in an open spectral gap or inside the spectrum. It contains the proofs of Theorem 1.2 (i) and (ii).

• In Section 5 we obtain the Landau-Peierls formula by taking the low density limit. It contains the proof of Theorem 1.2 (iii).

The Fermi energy

This section, which can be read independently of the rest of the paper, is only concerned with the location of the Fermi energy when the intensity of the magnetic field is zero (i.e. ω = 0). In particular, we prove Theorem 1.1. Although we assumed in the introduction that V ∈ C ∞ (T 3 ), all results of this section can be extended to the case V ∈ L ∞ (T 3 ).

Some preparatory results

Let ξ → f(β, µ; ξ) := ln 1 + e β(µ-ξ) be a holomorphic function on the domain {ξ ∈ C : ℑξ ∈ (-π/β, π/β)}. Let Γ the positively oriented simple contour included in the above domain defined by: Γ

:= ℜξ ∈ [δ, ∞), ℑξ = ± π 2β ∪ ℜξ = δ, ℑξ ∈ - π 2β , π 2β , (2.1)
where δ is any real number smaller than

E 0 := inf σ(H ∞ (0)) ≤ inf σ(H ∞ (ω)).
In the following we use δ := E 0 -1.

The thermodynamic limit of the grand-canonical density at β > 0, µ ∈ R and ω ≥ 0 is given by (see e.g. [START_REF] Briet | Do Bosons Condense in a Homogeneous Magnetic Field ?[END_REF]):

ρ ∞ (β, e βµ , ω) = 1 |Ω| i 2π Tr L 2 (R 3 ) χ Ω Γ dξ f F D (β, µ; ξ)(H ∞ (ω) -ξ) -1 (2.2)
where

f F D (β, µ; ξ) := -β -1 ∂ ξ f(β, µ; ξ) = (e β(ξ-µ) + 1) -1
is the Fermi-Dirac distribution function and χ Ω denotes the characteristic function of Ω. We prove in [START_REF] Savoie | [END_REF] (even for singular potentials) that ρ ∞ (β, •, ω) can be analytically extended to the domain C \ (-∞, -e βE0(ω) ]. Now assume that the intensity of the magnetic field is zero (ω = 0). The following proposition (stated without proof since the result is well known), allows us to rewrite (2.2) only using the Bloch energy functions k → E j (k) of H ∞ (0): Proposition 2.1. Let β > 0 and µ ∈ R. Denote by Ω * the first Brillouin zone of the dual lattice 2πZ 3 . Then:

ρ ∞ (β, e βµ , 0) = 1 (2π) 3 ∞ j=1 Ω * dk f F D (β, µ; E j (k)).
(2.3)

Note that another useful way to express the grand-canonical density at zero magnetic field consists in bringing into play the integrated density of states (IDS) of the operator H ∞ (0) (see (1.10) for its definition):

ρ ∞ (β, e βµ , 0) = - ∞ -∞ dλ ∂f F D ∂λ (β, µ; λ)n ∞ (λ). (2.4)
When the density of particles ρ 0 > 0 becomes the fixed parameter, the relation between the fugacity and density can be inverted. This is possible since for all β > 0, the map ρ ∞ (β, •, 0) is strictly increasing on (0, ∞) and defines a C ∞ -diffeomorphism of this interval onto itself. Then there exists an unique z ∞ (β, ρ 0 ) ∈ (0, ∞) and therefore an unique µ ∞ (β, ρ 0 ) ∈ R satisfying: ρ ∞ (β, e βµ∞(β,ρ0) , 0) = ρ 0 .

(2.5)

We now are interested in the zero temperature limit. The following proposition (again stated without proof) is a well known, straightforward consequence of the continuity of n ∞ (•): Proposition 2.2. Let µ ≥ E 0 := inf σ(H ∞ (0)) be fixed. We have the identity:

lim β→∞ ρ ∞ (β, e βµ , 0) = 1 (2π) 3 ∞ j=1 Ω * dk χ [E0,µ] (E j (k)) = n ∞ (µ), (2.6)
where

χ [E0,µ] (•) denotes the characteristic function of the interval [E 0 , µ].
We end this paragraph with another preparatory result concerning the behavior of n ∞ near the edges of a spectral gap. This result is contained in the following lemma:

Lemma 2.3. Let ρ 0 > 0 be fixed. Assume that there exists N ≥ 1 such that n ∞ (E) = ρ 0 for all E satisfying max E N ≤ E ≤ min E N +1 . We set a N := max E N and b N := min E N +1 . Assume that the gap is open, i.e. a N < b N . Then for δ > 0 sufficiently small, there exists a constant C = C δ > 0 such that: n ∞ (a N ) -n ∞ (λ) ≥ C(a N -λ) 3 whenever λ ∈ [a N -δ, a N ] (2.7) and n ∞ (λ) -n ∞ (b N ) ≥ C(λ -b N ) 3 whenever λ ∈ [b N , b N + δ]. (2.8)
Proof . We only prove (2.7), since the proof of the other inequality (2.8) is similar. Since a N = max k∈Ω * E N (k), the maximum is attained in a (possibly not unique) point k 0 , i.e. a N = E N (k 0 ). This means that a N is a discrete eigenvalue of finite multiplicity 1 ≤ M ≤ N of the fiber operator

h(k 0 ) = 1 2 (-i∇ + k 0 ) 2 + V .
In particular, a N is isolated from the rest of the spectrum since we assumed that a N < b N ≤ E N +1 (k 0 ). Now when k slightly varies around k 0 , the eigenvalue a N will split into at most M different eigenvalues, the largest of which being E N (k). Thus from the second equality in (2.6) we obtain:

n ∞ (a N ) -n ∞ (λ) ≥ 1 (2π) 3 Vol{k ∈ Ω * : λ ≤ E N (k) ≤ a N }.
We now choose δ small enough such that

σ(h(k 0 )) ∩ [a N -δ, a N + δ] = {E N (k 0 )}.
We use analytic perturbation theory in order to control the location of the spectrum of h(k) when |kk 0 | is small (we assume without loss of generality that k 0 lies in the interior of Ω * ). By writing

h(k) = h(k 0 ) + (k -k 0 ) • (-i∇ + k 0 ) + (k -k 0 ) 2 /2 =: h(k 0 ) + W (k),
we see that we can find a constant C > 0 such that

W (k)(h(k 0 ) -i) -1 ≤ C|k -k 0 |, |k -k 0 | ≤ 1.
Take a circle γ with center at a N and radius r = (a N -λ)/2 ≤ δ/2. For any z ∈ γ, by virtue of the first resolvent equation:

(h(k 0 ) -z) -1 = (h(k 0 ) -i) -1 + (z -i)(h(k 0 ) -i) -1 (h(k 0 ) -z) -1
and by using the estimate (h(k 0 ) -z) -1 = 2/(a N -λ), we can find another constant C δ > 0 such that:

sup z∈γ W (k)(h(k 0 ) -z) -1 ≤ C δ |k -k 0 | (a N -λ) , |k -k 0 | ≤ 1. It turns out that if |k -k 0 |/(a N -λ) is smaller than some ǫ > 0, then sup z∈γ W (k)(h(k 0 ) -z) -1 ≤ ǫC δ whenever |k -k 0 | ≤ (a N -λ)ǫ.
Standard analytic perturbation theory insures that if ǫ is chosen small enough, h(k) will have exactly M eigenvalues inside γ. Thus for all k satisfying |k -

k 0 | ≤ ǫ(a N -λ), we have σ(h(k)) ∩ [a N -δ, a N ] ⊆ [(a N + λ)/2, a N ] ⊂ [λ, a N ].
In particular, λ < E N (k) ≤ a N for all such k's. But the ball in Ω * where |kk 0 | ≤ (a N -λ)ǫ has a volume which goes like (a N -λ) 3 , and the proof is finished.

Proof of Theorem 1.1

In this paragraph we prove the existence of the Fermi energy. We separately investigate the semiconducting case and the metallic case.

The semiconducting case (SC)

We here consider the same situation as in Lemma 2.3 in which there exists N ≥ 1 such that n ∞ (E) = ρ 0 for all E satisfying max E N ≤ E ≤ min E N +1 . We set a N := max E N and b N := min E N +1 . Let µ(β) := µ ∞ (β, ρ 0 ) be the unique solution of the equation ρ ∞ (β, e βµ , 0) = ρ 0 . We start with the following lemma: Lemma 2.4.

a N ≤ µ 1 := lim inf β→∞ µ(β) ≤ lim sup β→∞ µ(β) =: µ 2 ≤ b N .
(2.9)

Proof . We will only prove the inequality a N ≤ µ 1 , since the proof of the other one (µ 2 ≤ b N ) is similar. Assume the contrary: µ 1 < a N . Define ǫ := a N -µ 1 > 0. Then there exists a sequence {β n } n≥1 with β n → ∞ and an integer M ǫ ≥ 1 large enough such that:

lim n→∞ µ(β n ) = µ 1 and µ(β n ) ≤ a N -ǫ/2 < a N , ∀ n ≥ M ǫ .
Since ρ ∞ (β, e βµ , 0) is an increasing function of µ, we have:

ρ 0 = ρ ∞ (β n , e βnµ(βn) , 0) ≤ ρ ∞ (β n , e βn(aN -ǫ/2) , 0).
By letting n → ∞ in the above inequality, (2.6) implies:

ρ 0 ≤ n ∞ (a N -ǫ/2) < n ∞ (a N ) = ρ 0
where in the second inequality we used (2.7). We have arrived at a contradiction. Now if a N = b N , the proof of (1.13) is over. Thus we can assume that a N < b N , i.e. the gap is open. We have the following lemma:

Lemma 2.5. Define c N = (a N + b N )/2. For any 0 < ǫ < (b N -a N )/2, there exists β ǫ > 0 large enough such that µ(β) ∈ [c N -ǫ, c N + ǫ] whenever β > β ǫ .
Proof . We know that µ(β) exists and is unique, thus if we can construct such a solution in the given interval, it means that this is the one. We use (2.4) in which we introduce µ(β) and arrive at the following identities:

n ∞ (a N ) = ρ 0 = ∞ -∞ dλ ∂f F D ∂λ (β, µ(β); λ)n ∞ (λ) = - aN -∞ dλ ∂f F D ∂λ (β, µ(β); λ)n ∞ (λ) -n ∞ (b N )f F D (β, µ(β); b N ) + n ∞ (a N )f F D (β, µ(β); a N ) - ∞ bN dλ ∂f F D ∂λ (β, µ(β); λ)n ∞ (λ),
where in the last term we used the fact that n ∞ (•) is constant on the interval [a N , b N ], and this constant is nothing but ρ 0 . We can rewrite the above equation as:

aN -∞ dλ ∂f F D ∂λ (β, µ(β); λ){n ∞ (λ)-n ∞ (a N )} = ∞ bN dλ ∂f F D ∂λ (β, µ(β); λ){n ∞ (b N )-n ∞ (λ)} (2.10)
where we used the fact that

f F D (β, µ(β); -∞) = 1 and f F D (β, µ(β); λ) ≤ Ce -λβ for large λ.
In the left hand side of (2.10) we now introduce the explicit formula:

∂ λ f F D (β, µ(β); λ) = -β e β(λ-µ(β)) (e β(λ-µ(β)) + 1) 2 = -βe β(aN -µ(β)) e β(λ-aN ) (e β(λ-µ(β)) + 1) 2 ,
while in the right hand side of (2.10) we use another expression:

∂ λ f F D (β, µ(β); λ) = -β e -β(λ-µ(β)) (1 + e -β(λ-µ(β)) ) 2 = -βe -β(bN -µ(β)) e -β(λ-bN ) (1 + e -β(λ-µ(β)) ) 2 .
Then (2.10) can be rewritten as:

aN -∞ dλ e β(λ-aN ) (e β(λ-µ(β)) + 1) 2 {n ∞ (a N ) -n ∞ (λ)} = e β{2µ(β)-(aN +bN )} ∞ bN dλ e -β(λ-bN ) (1 + e -β(λ-µ(β)) ) 2 {n ∞ (λ) -n ∞ (b N )}, (2.11)
or by taking the logarithm:

µ(β) = c N + 1 2β ln aN -∞ dλ e β(λ-aN ) (e β(λ-µ(β)) + 1) 2 {n ∞ (a N ) -n ∞ (λ)} -ln ∞ bN dλ e -β(λ-bN ) (1 + e -β(λ-µ(β)) ) 2 {n ∞ (λ) -n ∞ (b N )} .
Let us define the smooth function f : [c N -ǫ, c N + ǫ] → R given by:

f (x) := c N + 1 2β ln aN -∞ dλ e β(λ-aN ) (e β(λ-x) + 1) 2 {n ∞ (a N ) -n ∞ (λ)} -ln ∞ bN dλ e -β(λ-bN ) (1 + e -β(λ-x) ) 2 {n ∞ (λ) -n ∞ (b N )} .
(2.12)

We will prove that if β is large enough, then f invariates the interval [c N -ǫ, c N + ǫ], which is already enough for the existence of a fixed point. This would also show that µ(β) must be in that interval. But in fact one can prove more: f is a contraction for large enough β.

The idea is to find some good upper and lower bounds when β is large for the integrals under the logarithms. We start by finding a lower bound in β for the first integral. Let δ > 0 sufficiently small. Using (2.7) in the left hand side of (2.11) we get:

aN -∞ dλ e β(λ-aN ) (e β(λ-x) + 1) 2 {n ∞ (a N ) -n ∞ (λ)} ≥ C 4 aN aN -δ e -β(aN -λ) (a N -λ) 3 (2.13)
where we used that x ≥ a N ≥ λ in order to get rid of the numerator. After a change of variables and using some basic estimates one arrives at another constant C > 0 such that for β sufficiently large:

aN -∞ dλ e β(λ-aN ) (e β(λ-x) + 1) 2 {n ∞ (a N ) -n ∞ (λ)} ≥ C β 5 .
(2.14)

By restricting the interval of integration to [b N , b N + δ] and by using (2.8), we obtain by the same method a similar lower bound for the second integral under the logarithm. Moreover, using the Weyl asymptotics which says that n ∞ (λ) ∼ λ 3 2 for large λ (see e.g. [START_REF] Kirsch | Random Schrödinger Operators[END_REF]), one can also get a power-like upper bound in β for our two integrals.

We deduce from these estimates that there exists a constant C ǫ > 0 such that:

sup x∈[cN -ǫ,cN +ǫ] |f (x) -c N | ≤ C ǫ ln(β) β , β > 1.
Thus if β is large enough, f invariates the interval. Being continuous, it must have a fixed point. Moreover, the derivative f ′ (x) decays exponentially with

β uniformly in x ∈ [c N -ǫ, c N + ǫ]. It implies that if β is large enough, then f ′ ∞ < 1, that is f is a contraction.

The metallic case (M)

Consider the situation in which there exists a unique solution E M of the equation n ∞ (E M ) = ρ 0 , and this solution lies in the interior of a Bloch band. In other words, there exists (a possibly not unique) integer N ≥ 1 such that min E N < E M < max E N . We will use in the following that the IDS n ∞ (•) is a strictly increasing function on the interval [min E N , max E N ].

Let µ(β) := µ ∞ (β, ρ 0 ) be the unique real solution of the equation ρ ∞ (β, e βµ(β) , 0) = ρ 0 . Let us show that :

E M ≤ lim inf β→∞ µ(β) ≤ lim sup β→∞ µ(β) ≤ E M , (2.15)
which would end the proof. We start with the first inequality.

Assume ad-absurdum that µ 1 := lim inf β→∞ µ(β) < E M . Then there exists ǫ > 0 and a sequence {β n } n≥1 satisfying β n → ∞ such that:

lim n→∞ µ(β n ) = µ 1 , µ(β n ) ≤ E M -ǫ, ∀n ≥ 1.
Since ρ ∞ (β, e βµ , 0) is increasing with µ, we have: (2.16) where in the last equality we used (2.6). But the inequality n

n ∞ (E M ) = ρ 0 = lim n→∞ ρ ∞ (β n , e βnµ(βn) , 0) ≤ lim n→∞ ρ ∞ (β n , e βn(EM -ǫ) , 0) = n ∞ (E M -ǫ),
∞ (E M ) ≤ n ∞ (E M -ǫ) is in contradiction with the fact that n ∞ (•) is a strictly increasing function near E M . Thus E M ≤ µ 1 .
Now assume ad-absurdum that µ 2 := lim sup β→∞ µ(β) > E M . Then there exists ǫ > 0 and a sequence {β n } n≥1 satisfying β n → ∞ such that:

lim n→∞ µ(β n ) = µ 2 , E M + ǫ ≤ µ(β n ), ∀n ≥ 1.
We again use that ρ ∞ (β, e βµ , 0) is increasing with µ and write: (2.17) where in the first equality we again used (2.6). But the inequality n

n ∞ (E M + ǫ) = lim n→∞ ρ ∞ (β n , e βn(EM +ǫ) , 0) ≤ lim n→∞ ρ ∞ (β n , e βnµ(βn) , 0) = ρ 0 = n ∞ (E M ),
∞ (E M + ǫ) ≤ n ∞ (E M ) is also in contradiction with the fact that n ∞ (•) is a strictly increasing function near E M . Therefore µ 2 ≤ E M .
3 The zero-field susceptibility at fixed density and positive temperature

In this section we prove a general formula for the zero-field grand-canonical susceptibility of a Bloch electrons gas at fixed density and positive temperature.

Here is the main result of this section:

Theorem 3.1. Let β > 0 and ρ 0 > 0 be fixed. Let µ ∞ = µ ∞ (β, ρ 0 ) ∈ R the unique solution of the equation ρ ∞ (β, e βµ , ω = 0) = ρ 0 . Then for each integer j 1 ≥ 1 there exists four families of functions c j1,l (• ), with l ∈ {0, 1, 2, 3}, defined on Ω * outside a set of Lebesgue measure zero, such that the integrand below can be extended by continuity to the whole of Ω * :

X (β, ρ 0 ) = - e c 2 1 2β 1 (2π) 3 ∞ j1=1 Ω * dk 3 l=0 ∂ l f ∂ξ l β, µ ∞ ; E j1 (k) c j1,l (k), (3.1)
with the convention

(∂ 0 ξ f)(β, µ ∞ ; E j1 (k)) = f(β, µ ∞ ; E j1 (k)) := ln(1 + e β(µ∞-Ej 1 (k)) ).
This formula is a necessary step in the proof of Theorem 1.2 (i) and (ii) (this is the aim of the following section) when we will take the limit of zero temperature.

The special feature of this formula lies in the fact that each function c j1,l (•) can be only expressed in terms of Bloch energy functions and their associated eigenfunctions. For each integer j 1 ≥ 1, the functions c j1,2 (•) and c j1,3 (•) are identified respectively in (3.28) and (3.27). As for the functions c j1,l (•) with l ∈ {0, 1}, they can also be written down but their explicit expression is not important for the proof of Theorem 1.2. Note as well that the above formula brings into play the Fermi-Dirac distribution and its partial derivatives up to the second order. This will turn out to be very important when we will take the limit β → ∞ in the following section.

Starting the proof: a general formula from the magnetic perturbation theory

We start by giving a useful formula for the thermodynamic limit of the grand-canonical susceptibility. Let β > 0 and z := e βµ ∈ (0, ∞) the fixed external parameters. Let Γ be the positively oriented contour defined in (2.1), going round the half-line [E 0 , ∞), and included in the analiticity domain of the map ξ → f(β, z; ξ) = ln(1 + ze -βξ ). Denote by R ∞ (ω, ξ) := (H ∞ (ω) -ξ) -1 for all ξ ∈ ρ(H ∞ (ω)) and ω ∈ R. Taking into account the periodic structure of our system, it is proved (see [START_REF] Briet | Diamagnetism of quantum gases with singular potentials[END_REF], Theorem 3.8) that the thermodynamic limit of the grand-canonical pressure of the Bloch electron gas at any intensity of the magnetic field B is given by:

P ∞ (β, z, ω) := 1 β|Ω| i 2π Tr L 2 (R 3 ) χ Ω Γ dξ f(β, z; ξ)R ∞ (ω, ξ) , (3.2)
where Ω is the unit cube centered at the origin of coordinates (χ Ω denotes its characteristic function). Although the integral kernel R ∞ (• , • ; ω, ξ) of the resolvent has a singularity on the diagonal, the integration with respect to ξ in (3.2) provides us with a jointly continuous kernel on R 3 ×R 3 . One can see this by performing an integration by parts in (3.2) and using the fact that the kernel of R 2 ∞ (ω, ξ) is jointly continuous. Moreover, one can prove [START_REF] Briet | Diamagnetic expansions for perfect quantum gases[END_REF][START_REF] Briet | Generalized susceptibilities for a perfect quantum gas. Markov Process[END_REF][START_REF] Briet | Diamagnetic expansions for perfect quantum gases II: uniform bounds[END_REF] that the thermodynamic limit of the grand-canonical pressure is jointly smooth on (z, ω) ∈ (-e βE0 , ∞) × R.

Let ω ∈ R and ξ ∈ ρ(H ∞ (ω)). Introduce the bounded operators T ∞,1 (ω, ξ) and T ∞,2 (ω, ξ) generated by the following integral kernels:

T ∞,1 (x, y; ω, ξ) := a(x -y) • (i∇ x + ωa(x))R ∞ (x, y; ω, ξ) (3.3)
T ∞,2 (x, y; ω, ξ) := 1 2 a 2 (x -y)R ∞ (x, y; ω, ξ), x = y, (3.4)
where a(•) stands for the usual symmetric gauge a(x) = 1 2 e 3 ∧ x = 1 2 (-x 2 , x 1 , 0). We introduce the following operators :

W ∞,1 (β, µ, ω) := i 2π Γ dξ f(β, µ; ξ)R ∞ (ω, ξ)T ∞,1 (ω, ξ)T ∞,1 (ω, ξ) (3.5) W ∞,2 (β, µ, ω) := i 2π Γ dξ f(β, µ; ξ)R ∞ (ω, ξ)T ∞,2 (ω, ξ) (3.6)
One can prove using the same techniques as in [START_REF] Cornean | The Faraday effect revisited: Thermodynamic limit[END_REF] that these operators are locally trace class and have a jointly continuous kernel on R 3 × R 3 . By a closely related method as the one in [START_REF] Briet | Diamagnetic expansions for perfect quantum gases[END_REF], [START_REF] Briet | Generalized susceptibilities for a perfect quantum gas. Markov Process[END_REF], it is proved in [START_REF] Savoie | [END_REF] that we can invert the thermodynamic limit with the partial derivatives w.r.t. ω of the grand-canonical pressure. Then the bulk orbital susceptibility reads as:

X GC ∞ (β, e βµ , ω) := e c 2 ∂ 2 P ∞ ∂ω 2 (β, e βµ , ω) = e c 2 2 β|Ω| Tr L 2 (R 3 ) χ Ω W ∞,1 (β, µ, ω) -Tr L 2 (R 3 ) χ Ω W ∞,2 (β, µ, ω)
We mention that the above formula is obtained using the so-called 'gauge invariant magnetic perturbation theory' applied to the resolvent integral kernel (see e.g. [START_REF] Cornean | The Faraday effect revisited: Thermodynamic limit[END_REF] for further details) which allows to control the linear growth of the magnetic vector potential.

The quantity which we are interested in is the orbital susceptibility at zero magnetic field and at fixed density of particles ρ 0 . Note that the pressure is an even function of ω, thus its first order derivative at ω = 0 is zero. This explains why the susceptibility is the relevant physical quantity for the weak magnetic field regime.

The orbital susceptibility at zero magnetic field and fixed density ρ 0 is given by (see also (1.8)):

X (β, ρ 0 ) : = X GC ∞ (β, e βµ∞(β,ρ0) , 0) = e c 2 2 β|Ω| Tr L 2 (R 3 ) χ Ω W ∞,1 (β, µ ∞ , 0) -Tr L 2 (R 3 ) χ Ω W ∞,2 (β, µ ∞ , 0) . (3.7)
The formula (3.7) constitutes the starting-point in obtaining (3.1). The next step consists in rewriting the local traces appearing in (3.7) in a more convenient way: Proposition 3.2. Let p α := -i∂ α with α ∈ {1, 2, 3} be the cartesian components of the momentum operator defined in L 2 (R 3 ). Then we have:

Tr L 2 (R 3 ) χ Ω W ∞,1 (β, µ ∞ , 0) = 1 4 i 2π Tr L 2 (R 3 ) χ Ω Γ dξ f(β, µ ∞ ; ξ) R ∞ (0, ξ)p 1 R ∞ (0, ξ)p 2 R ∞ (0, ξ) p 2 R ∞ (0, ξ)p 1 R ∞ (0, ξ) -p 1 R ∞ (0, ξ)p 2 R ∞ (0, ξ) + + R ∞ (0, ξ)p 2 R ∞ (0, ξ)p 1 R ∞ (0, ξ) p 1 R ∞ (0, ξ)p 2 R ∞ (0, ξ) -p 2 R ∞ (0, ξ)p 1 R ∞ (0, ξ) (3.8)
and

Tr L 2 (R 3 ) χ Ω W ∞,2 (β, µ ∞ , 0) = - 1 4 i 2π Tr L 2 (R 3 ) χ Ω Γ dξ f(β, µ ∞ ; ξ) R ∞ (0, ξ)R ∞ (0, ξ) p 2 R ∞ (0, ξ)p 2 R ∞ (0, ξ) + p 1 R ∞ (0, ξ)p 1 R ∞ (0, ξ) -R ∞ (0, ξ) . (3.9)
Proof . We begin with the justification of (3.9). By rewriting (3.4) as:

T ∞,2 (x, y; 0, ξ) = 1 8 {e 3 ∧ (x -y)} • {e 3 ∧ (x -y)}R ∞ (x, y; 0, ξ) = 1 8 (x 2 -y 2 ) 2 + (x 1 -y 1 ) 2 R ∞ (x, y; 0, ξ), from (3.6) it follows: W ∞,2 (x, x; β, µ, 0) (3.10) = 1 8 Γ dξ f(β, µ; ξ) R 3 dz R ∞ (x, z; 0, ξ) (z 2 -x 2 ) 2 + (z 1 -x 1 ) 2 R ∞ (z, x; 0, ξ), ∀ x ∈ R 3 .
Let l ∈ {1, 2}. Denote by X the multiplication operator with x. Then for all z = x we can write:

(z l -x l )R ∞ (z, x; 0, ξ) = X • e l , R ∞ (0, ξ) (z, x) = R ∞ (0, ξ) H ∞ (0), X • e l R ∞ (0, ξ) (z, x).
We know that [H ∞ (0), X • e l ] = -ip l . Thus:

(z l -x l )R ∞ (z, x; 0, ξ) = -i R ∞ (0, ξ)p l R ∞ (0, ξ) (z, x). (3.11)
Using standard commutation rules, we deduce from (3.11) that for l ∈ {1, 2} and for all z = x:

(z l -x l ) 2 R ∞ (z, x; 0, ξ) = -2R ∞ (0, ξ)p l R ∞ (0, ξ)p l R ∞ (0, ξ) -R ∞ (0, ξ)R ∞ (0, ξ) (z, x). (3.12)
It remains to put (3.12) in (3.10), and we get (3.9).

Let us now prove now (3.8). Since the divergence of a is zero, then for x = y we have:

T ∞,1 (x, y; 0, ξ) = i 2 ∇ x • {e 3 ∧ (x -y)}R ∞ (x, y; 0, ξ) = i∇ x • - (x 2 -y 2 ) 2 e 1 + (x 1 -y 1 ) 2 e 2 R ∞ (x, y; 0, ξ).
From (3.5) it follows that for all x ∈ R 3 :

W ∞,1 (x, x; β, µ, 0) = 1 4 Γ dξ f(β, µ; ξ) R 3 dz 1 R 3 dz 2 R ∞ (x, z 1 ; 0, ξ) (i∇ z1 • e 1 )[-(z 1,2 -z 2,2 )R ∞ (0, ξ)(z 1 , z 2 )] + (i∇ z1 • e 2 )[(z 1,1 -z 2,1 )R ∞ (0, ξ)(z 1 , z 2 )] • • (i∇ z2 • e 1 )[-(z 2,2 -x 2 )R ∞ (0, ξ)(z 2 , x)] + (i∇ z2 • e 2 )[(z 2,1 -x 1 )R ∞ (0, ξ)(z 2 , x)] .
Then by using (3.11), we get (3.8) from the following identity:

∀ x ∈ R 3 , W ∞,1 (x, x; β, µ, 0) = 1 4 Γ dξ f(β, µ; ξ) R 3 dz 1 R 3 dz 2 R ∞ (x, z 1 ; 0, ξ) ip 1 R ∞ (0, ξ)p 2 R ∞ (0, ξ) (z 1 , z 2 ) -ip 2 R ∞ (0, ξ)p 1 R ∞ (0, ξ) (z 1 , z 2 ) ip 1 R ∞ (0, ξ)p 2 R ∞ (0, ξ) (z 2 , x) -ip 2 R ∞ (0, ξ)p 1 R ∞ (0, ξ) (z 2 , x) .

Using the Bloch decomposition

We know that (see e.g. [START_REF] Berezin | The Schrödinger Equation[END_REF]) H ∞ (0) can be seen as a direct integral

⊕ Ω * dk h(k)
where the fiber Hamiltonians h(k) acting in L 2 (T 3 ) are given by :

h(k) = 1 2 (-i∇ + k) 2 + V. (3.13)
Recall that h(k) is essentially self-adjoint in C ∞ (T 3 ); the domain of its closure is the Sobolev space H 2 (T 3 ). For each k ∈ Ω * , h(k) has purely discrete spectrum. We have already denoted by {E j (k)} j≥1 the set of eigenvalues counting multiplicities and in increasing order. The corresponding eigenfunctions {u j (• ; k)} j≥1 form a complete orthonormal system in L 2 (T 3 ) and satisfy:

h(k)u j (•; k) = E j (k)u j (•; k).
The eigenfunctions u j 's are defined up to an arbitrary phase depending on k. These phases cannot be always chosen to be continuous at crossing points, and even less differentiable. For the following let us introduce another notation. Let α ∈ {1, 2, 3}, and let i, j ≥ 1 be any natural numbers. Then for all k ∈ Ω * we define: πi,j (α; k) :=

Ω dx u i (x; k)[(p α + k α )u j (x; k)] = u i (• ; k), (p α + k α )u j (• ; k) . (3.14)
Note that due to the phases presence in the eigenfunctions u j 's, we cannot be sure that the πi,j 's are continuous/differentiable at crossing points. But all these 'bad' phase factors will disappear when we take the traces (see (3.20) and (3.23) below).

We now can write the local traces of Proposition 3.2 in the following way: Proposition 3.3. Let β > 0 and ρ 0 > 0 be fixed. Let µ ∞ = µ ∞ (β, ρ 0 ) ∈ R be the unique solution of the equation ρ ∞ (β, e βµ , 0) = ρ 0 . Then both quantities (3.8) and (3.9) can be rewritten as:

Tr L 2 (R 3 ) χ Ω W ∞,1 (β, µ ∞ , 0) = - 1 4 1 |Ω * | ∞ j1,...,j4=1 Ω * dk C j1,j2,j3,j4 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 2 E j2 (k) -ξ E j3 (k) -ξ E j4 (k) -ξ , (3.15)
and

Tr L 2 (R 3 ) χ Ω W ∞,2 (β, µ ∞ , 0) = - 1 4 1 |Ω * | ∞ j1=1 Ω * dk 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 3 + - ∞ j1,j2=1 Ω * dk C j1,j2 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 3 E j2 (k) -ξ , (3.16)
where the functions Ω * ∋ k → C j1,j2,j3,j4 (k) and Ω * ∋ k → C j1,j2 (k) are defined by:

C j1,j2,j3,j4 (k) := πj1,j2 (1; k)π j2,j3 (2; k) -πj1,j2 (2; k)π j2,j3 (1; k) × πj3,j4 (2; k)π j4,j1 (1; k) -πj3,j4 (1; k)π j4,j1 (2; k) (3.17)
and

C j1,j2 (k) := πj1,j2 (1; k)π j2,j1 (1; k) + πj1,j2 (2; k)π j2,j1 (2; k) = πj1,j2 (1; k) 2 + πj1,j2 (2; k) 2 . (3.18)
We do not give more details since this result is just a straightforward application of the following rather non-trivial technical lemma (recently proved in [START_REF] Cornean | Faraday effect revisited: sum rules and convergence issues[END_REF]): Lemma 3.4. Let β > 0 and µ ∈ R be fixed. For n, m ∈ N with m, n ≥ 1, consider the local trace given by:

J (m) α1,...,αn := Tr L 2 (R 3 ) χ Ω Γ dξ f(β, µ; ξ)(H ∞ (0) -ξ) -m p α1 (H ∞ (0) -ξ) -1 • • • p αn (H ∞ (0) -ξ) -1
Then under the assumption that V ∈ C ∞ (T 3 ) we have:

J (m) α1,...,αn = 1 |Ω * | j1,...,jn≥1 Ω * dk πj1,j2 (α 1 ; k) . . . πjn,j1 (α n ; k) Γ dξ f(β, µ; ξ) E j1 (k) -ξ m+1 E j2 (k) -ξ • • • E jn (k) -ξ . (3.19)
where all the above series are absolutely convergent and πi,j (α; k) is defined by (3.14).

Applying the residue calculus

Consider the expression of the susceptibility at fixed density (3.7) in which the local traces are now given by (3.15) and (3.16). Remark that these quantities now are written in a convenient way in order to apply the residue theorem. Denote the integrands appearing in (3.15) and (3.16) by:

g j1,j2 (β, µ ∞ ; ξ) := f(β, µ ∞ ; ξ) E j1 (k) -ξ 3 E j2 (k) -ξ , j 1 , j 2 ∈ N * h j1,j2,j3,j4 (β, µ ∞ ; ξ) := f(β, µ ∞ ; ξ) E j1 (k) -ξ 2 E j2 (k) -ξ E j3 (k) -ξ E j4 (k) -ξ , j 1 , j 2 , j 3 , j 4 ∈ N * .
Note that g j1,j2 (β, µ ∞ ; • ) can have first order, third order, or even fourth order poles (in the case when j 1 = j 2 ). In the same way, h j1,j2,j3,j4 (β, µ ∞ ; • ) can have poles from the first order up to at most fifth order (in the case when j 1 = j 2 = j 3 = j 4 ). Hence we expect that the integrals of

h j1,j2,j3,j4 (β, µ ∞ ; • ) in (3.15) (resp. of g j1,j2 (β, µ ∞ ; • ) in (3. 16
)) to make appear partial derivatives of f(β, µ ∞ ; • ) with order at most 4 (resp. with order at most 3). But we will see below that the factor multiplying (∂ 4 ξ f)(β, µ ∞ ; • ) is identically zero. Getting back to the susceptibility formula in (3.7) and by virtue of the previous remarks, we expect to obtain an expansion of the orbital susceptibility of the type (3.1). The next two results identify the functions c j1,l (• ) coming from (3.15) and (3.16): Lemma 3.5. The quantity defined by (3.15) can be rewritten as:

Tr L 2 (R 3 ) χ Ω W ∞,1 (β, µ ∞ ) = - 1 4 1 |Ω * | ∞ j1=1 Ω * dk 3 l=0 ∂ l f ∂ξ l β, µ ∞ ; E j1 (k) a j1,l (k) (3.20)
where for all j 1 ∈ N * and k ∈ Ω * , the functions a j1,3 (•) and a j1,2 (•) are given by:

a j1,3 (k) := 1 3! πj1,j1 (1; k) 2 ∞ j2=1 j2 =j1 πj1,j2 (2; k) 2 E j2 (k) -E j1 (k) + πj1,j1 (2; k) 2 ∞ j2=1 j2 =j1 πj1,j2 (1; k) 2 E j2 (k) -E j1 (k) -πj1,j1 (1; k)π j1,j1 (2; k) ∞ j2=1 j2 =j1 2ℜ πj1,j2 (2; k)π j2,j1 (1; k) E j2 (k) -E j1 (k) (3.21)
and

a j1,2 (k) := - 1 2! ∞ j2=1 j2 =j1 ∞ j3=1 j3 =j1 C j1,j1,j2,j3 (k) + C j1,j2,j1,j3 (k) + C j1,j2,j3,j1 (k) E j2 (k) -E j1 (k) E j3 (k) -E j1 (k) + ∞ j2=1 j2 =j1 C j2,j1,j1,j1 (k) -C j1,j1,j2,j1 (k) E j2 (k) -E j1 (k) 2 . (3.22)
Note that it is possible to identify in (3.20) all the functions a j1,l (•) for j 1 ≥ 1 and l ∈ {1, 0} since such a result is based only on identities provided by the residue theorem. However, the number of terms is large and we will not need their explicit expressions in order to prove our theorem. Now we treat the next term.

Lemma 3.6. The quantity defined by (3.16) can be rewritten as:

Tr L 2 (R 3 ) χ Ω W ∞,2 (β, µ ∞ ) = 1 4 1 |Ω * | ∞ j1=1 Ω * dk 3 l=0 ∂ l f ∂ξ l β, µ ∞ ; E j1 (k) b j1,l (k) (3.23)
where for all integers j 1 ≥ 1 and all k ∈ Ω * we have:

b j1,3 (k) := 1 6 πj1,j1 (1; k) 2 + πj1,j1 (2; k) 2 , (3.24) b j1,2 (k) := - 1 2 ∞ j2=1 j2 =j1 πj1,j2 (1; k) 2 + πj1,j2 (2; k) 2 E j2 (k) -E j1 (k) + 1 2 , (3.25) b j1,s (k) := -(2 -s) ∞ j2=1 j2 =j1 πj1,j2 (1; k) 2 + πj1,j2 (2; k) 2 E j2 (k) -E j1 (k) 3-s , s ∈ {0, 1}.
Thus our Lemmas 3.5 and 3.6 provide an expansion of the type announced in (3.1), where the coefficients are given by:

c j1,l (k) := a j1,l (k) + b j1,l (k), l ∈ {0, 1, 2, 3}. (3.26)
In particular, for all integer j 1 ≥ 1 and for all k ∈ Ω * , the functions c j1,3 (• ) and c j1,2 (• ) are respectively given by:

c j1,3 (k) := 1 3! |π j1,j1 (1; k)| 2 1 + ∞ j2=1 j2 =j1 |π j1,j2 (2; k)| 2 E j2 (k) -E j1 (k) + + |π j1,j1 (2; k)| 2 1 + ∞ j2=1 j2 =j1 |π j1,j2 (1; k)| 2 E j2 (k) -E j1 (k) + -πj1,j1 (1; k)π j1,j1 (2; k) ∞ j2=1 j2 =j1 2ℜ πj1,j2 (2; k)π j2,j1 (1; k) E j2 (k) -E j1 (k) (3.27) and c j1,2 (k) := - 1 2 ∞ j2=1 j2 =j1 |π j1,j2 (1; k)| 2 + |π j1,j2 (2; k)| 2 E j2 (k) -E j1 (k) -1 + ∞ j2=1 j2 =j1 C j2,j1,j1,j1 (k) -C j1,j1,j2,j1 (k) (E j2 (k) -E j1 (k)) 2 + + ∞ j2=1 j2 =j1 ∞ j3=1 j3 =j1 C j1,j1,j2,j3 (k) + C j1,j2,j1,j3 (k) + C j1,j2,j3,j1 (k) (E j2 (k) -E j1 (k))(E j3 (k) -E j1 (k)) , (3.28)
where for all integers j 1 , j 2 , j 3 , j 4 ∈ N * , Ω * ∋ k → C j1,j2,j3,j4 (k) is defined in (3.17).

In order to conclude the proof of Theorem 3.1, it remains to use this last result (its proof is in the appendix of this section): Lemma 3.7. For all integers j 1 ≥ 1 and l ∈ {0, 1, 2, 3}, the maps Ω * ∋ k → a j1,l (k) and Ω * ∋ k → b j1,l (k) are bounded and continuous on any compact subset of Ω * where E j1 is isolated from the rest of the spectrum.

Thus for all integers j 1 ≥ 1 and k ∈ Ω * , the maps c j1,l (• ) appearing in (3.1) might be singular on a set with zero Lebesgue measure where E j1 can touch the neighboring bands. However, the whole integrand in (3.1) is bounded and continuous on the whole Ω * because it comes from two complex integrals ((3.15) and (3.16)) which do not have local singularities in k.

Appendix -Proofs of the intermediate results

Here we prove Lemmas 3.5, 3.6, and 3.7. Proof of Lemma 3.5. Let Ω * ∋ k → C j1,j2,j3,j4 (k) be the complex-valued function appearing in (3.15):

C j1,j2,j3,j4 (k) := πj1,j2 (1; k)π j2,j3 (2; k) -πj1,j2 (2; k)π j2,j3 (1; k) × πj3,j4 (2; k)π j4,j1 (1; k) -πj3,j4 (1; k)π j4,j1 (2; k) . (3.29)
Note that this function is identically zero for the following combinations of subscripts:

j 1 = j 2 = j 3 = j 4 , j 1 = j 2 = j 3 = j 4 , j 1 = j 3 = j 4 = j 2 . (3.30)
Therefore the expansion of (3.15) consists of partial derivatives of f(β, µ ∞ ; • ) of order at most equal to three. On the other hand, since the functions C j1,j1,j1,j4 (• ) and C j1,j2,j1,j1 (• ) are identically equal to zero (see (3.30)), the quadruple summation in (3.15) is reduced to :

∞ j1,...,j4=1 C j1,j2,j3,j4 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 2 E j2 (k) -ξ E j3 (k) -ξ E j4 (k) -ξ = ∞ j1=1 ∞ j3=1 j3 =j1 C j1,j1,j3,j1 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 4 E j3 (k) -ξ + ∞ j1=1 ∞ j2=1 j2 =j1 C j1,j2,j2,j2 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 2 E j2 (k) -ξ 3 + ∞ j1,...,j4=1 at most 2 equal subscripts C j1,j2,j3,j4 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 2 E j2 (k) -ξ E j3 (k) -ξ E j4 (k) -ξ . (3.31)
By applying the residue theorem in the first term of the right hand side of (3.31) we get:

∞ j3=1 j3 =j1 C j1,j1,j3,j1 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 4 E j3 (k) -ξ = ∞ j3=1 j3 =j1 C j1,j1,j3,j1 (k) 1 3! 1 E j3 (k) -E j1 (k) ∂ 3 f ∂ξ 3 β, µ ∞ ; E j1 (k) + 3 3! 1 E j3 (k) -E j1 (k) 2 ∂ 2 f ∂ξ 2 β, µ ∞ ; E j1 (k) + + others terms involving ∂ l f ∂ξ l (β, µ ∞ ; •), with l ≤ 1 . (3.32)
The function C j1,j1,j3,j1 (• ) appearing in front of

∂ 3 f ∂ξ 3 β, µ ∞ ; E j1 (k) in (3.32) corresponds to a j1,3 (• ) since: ∀ k ∈ Ω * , C j1,j1,j3,j1 (k) = πj1,j1 (1; k)π j1,j3 (2; k) -πj1,j1 (2; k)π j1,j3 (1; k) 2 .
Note that the function C j1,j1,j3,j1 (•) contributes to the term a j1,2 (•), too. By applying once again the residue theorem in the second term of the right hand side of (3.31) we obtain:

∞ j2=1 j2 =j1 C j1,j2,j2,j2 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 2 E j3 (k) -ξ 3 = ∞ j2=1 j2 =j1 C j1,j2,j2,j2 (k) - 1 2! 1 E j1 (k) -E j2 (k) 2 ∂ 2 f ∂ξ 2 β, µ ∞ ; E j2 (k) + others terms involving ∂ l f ∂ξ l (β, µ ∞ ; •), with l ≤ 1 . (3.33)
The function C j1,j2,j2,j2 (• ) appearing in front of ∂ 2 f ∂ξ 2 β, µ ∞ ; E j2 (k) contributes to a j1,2 (• ). It remains to isolate in (3.31) (where at most two subscripts are equal) all combinations which provide a second order derivative of f(β, µ ∞ ; • ). These combinations are:

j 1 = j 2 = j 3 , j 4 ; j 1 = j 3 = j 2 , j 4 ; j 1 = j 4 = j 2 , j 3 .
Finally, we once again apply the residue theorem and gathering all terms proportional with

∂ 2 f ∂ξ 2 (β, µ ∞ ; • ).
The proof is over. Proof of Lemma 3.6. By separating the cases j 1 = j 2 and j 1 = j 2 , the double summation in the right hand side of (3.16) reads as:

∞ j1=1 ∞ j2=1 C j1,j2 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 3 E j2 (k) -ξ = ∞ j1=1 C j1,j1 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 4 + ∞ j1=1 ∞ j2=1 j2 =j1 C j1,j2 (k) 1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k)(k) -ξ 3 E j2 (k) -ξ . (3.34)
By using the residue theorem in the first term of the r.h.s. of (3.34) :

1 2iπ Γ dξ f(β, µ ∞ ; ξ) E j1 (k) -ξ 4 = 1 3! ∂ 3 f ∂ξ 3 β, µ ∞ ; E j1 (k)
This is only the one term which provides a third-order partial derivative of f(β, µ ∞ ; • ). The rest of the proof is just plain computation using the residue theorem. We do not give further details.

Proof of Lemma 3.7. Let p α := -i∂ α be the α component of the momentum operator with periodic boundary conditions in L 2 (Ω), α ∈ {1, 2, 3}. Now assume that E j1 (k) is isolated and non-degenerate if k belongs to some compact K ⊂ Ω * . We have to investigate integrals of the type

Tr L 2 (Ω) Γ dξ f(β, µ ∞ ; ξ)(h(k) -ξ) -1 p α1 (h(k) -ξ) -1 • • • p α4 (h(k) -ξ) -1 . (3.35)
Let k 0 ∈ K, and let Γ 1 be a simple, positively oriented path surrounding E j1 (k 0 ) but no other eigenvalue of h(k 0 ). If |kk 0 | is small enough, Γ 1 will still only contain E j1 (k). The projection Π(k) corresponding to E j1 (k) is given by a Riesz integral. We have:

Π(k) = i 2π Γ1 dz (h(k) -z) -1 , (3.36)
and is continuous at k 0 in the trace norm topology. Moreover,

Π(k)(h(k) -ξ) -1 = 1 E j1 (k) -ξ Π(k), (1 -Π(k))(h(k) -ξ) -1 = 1 2πi Γ1 dz 1 z -ξ (h(k) -z) -1 . (3.37)
Clearly, Π(k)(h(k) -ξ) -1 is analytic in ξ in the exterior of Γ 1 . We can decompose the integral on Γ in (3.35) as a sum of three integrals, one of which being on a simple contour Γ 2 around E j1 (k 0 ), completely surrounded by Γ 1 . The other two integrals will never have E j1 (k) as a singularity, so they cannot contribute to the formula of a j1,l (k). On the other hand, in the integral on Γ 2 we can replace the resolvents with the decomposition in (3.37) and use the fact that (1-Π(k))(h(k)-ξ) -1 is analytic if ξ lies inside Γ 2 . Now one can apply the Cauchy residue formula. For example, we can compute the integral in which we have Π(k) at the extremities, and (1 -Π(k)) in the interior; in that case E j1 = E j1 (k) will be a double pole:

Tr L 2 (Ω) Γ2 dξ f(ξ)Π(k)(h(k) -ξ) -1 p α1 (h(k) -ξ) -1 (1 -Π(k)) • • • p α4 (h(k) -ξ) -1 Π(k) = 2πi (∂ ξ f)(E j1 (k))Tr L 2 (Ω) Π(k)p α1 (h(k) -E j1 ) -1 (1 -Π(k)) • • • p α4 Π(k) + (3.38) + f(E j1 (k)) d dξ Tr L 2 (Ω) Π(k)p α1 (h(k) -ξ) -1 (1 -Π(k)) • • • (h(k) -ξ) -1 (1 -Π(k))p α4 ξ=Ej 1 (k) .
Thus one contribution to a j1,1 (k) will be:

Tr L 2 (Ω) Π(k)p α1 (h(k) -E j1 ) -1 (1 -Π(k)) • • • p α4 Π(k) .
This expression does not use eigenvectors, only resolvents and projectors. Since E j1 is continuous at k 0 , the map

k → (1 -Π(k))(h(k) -E j1 (k)) -1 = 1 2πi Γ1 dz 1 z -E j1 (k) (h(k) -z) -1
is operator norm continuous at k 0 , and the map k → Π(k) is continuous in the trace norm. By using standard perturbation theory (see e.g. [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF]), the same holds for the maps:

k → (1 -Π(k))(h(k) -E j1 (k)) -1 p α l and k → p α l (1 -Π(k))(h(k) -E j1 (k)) -1 p α k .
Thus the trace defines a continuous function; all other coefficients can be treated in a similar way.

4 The zero-field susceptibility at fixed density and zero temperature

In this section, we separately investigate the semiconducting and metallic cases from the expansion (3.1). In particular, we prove Theorem 1.2 (i) and (ii).

4.1 The semiconducting case (SC)-Proof of Theorem 1.2 (i) By using that f F D (β, µ; ξ) = -β -1 ∂ ξ f(β, µ; ξ), (3.1) can be rewritten as:

X (β, ρ 0 ) = e c 2 1 2 1 (2π) 3 ∞ j1=1 Ω * dk 2 l=0 ∂ l f F D ∂ξ l β, µ ∞ ; E j1 (k) c j1,1+l (k) - 1 β f β, µ ∞ ; E j1 (k) c j1,0 (k) . (4.1)
From (4.1), the proof of Theorem 1.2 (i) is based on two main ingredients. The first one is that for any fixed µ ≥ E 0 we have the following pointwise convergences:

lim β→∞ 1 β f(β, µ; ξ) = (µ -ξ)χ [E0,µ] (ξ), lim β→∞ f F D (β, µ; ξ) = χ [E0,µ] (ξ), ∀ ξ ∈ [E 0 , ∞) \ {µ}, (4.2
) while in the distributional sense:

lim β→∞ ∂f F D ∂ξ (β, µ; ξ) = -δ(ξ -µ), lim β→∞ ∂ 2 f F D ∂ξ 2 (β, µ; ξ) = -∂ ξ δ(ξ -µ). (4.3)
The second ingredient is related to the decay of the derivatives of the Fermi-Dirac distribution: for all d > 0 and for all j ∈ N * , there exists a constant C j,d > 0 such that sup |ξ-µ|≥d>0

∂ j f F D ∂ξ j (β, µ; ξ) ≤ C j,d e -β|ξ-µ| 2 . ( 4.4)
Now assume that we are in the semiconducting case with a non-trivial gap, that is there exists N ∈ N * such that lim β→∞ µ ∞ (β, ρ 0 ) = (max E N + min E N +1 )/2 = E F (ρ 0 ) and max E N < min E N +1 . Since the Fermi energy lies inside a gap, all terms containing derivatives of the Fermi-Dirac distribution will converge to zero in the limit β → ∞. Here (4.4) plays a double important role: first, it makes the series in j 1 convergent, and second, it provides an exponential decay to zero. Then by taking into account (4.2), we immediately get (1.15) from (4.1).

The metallic case (M)-Proof of Theorem 1.2 (ii)

Now we are interested in the metallic case. The limit β → ∞ is not so simple as in the previous case, because the Fermi energy lies in the spectrum. The starting point is the same formula (3.1), but we have to modify it by getting rid of the third order partial derivatives of f in order to make appear a Landau-Peierls type contribution. However, this operation needs the already announced additional assumption of non-degeneracy (which will provide regularity in k) in a neighborhood of the Fermi surface: Assumption 4.1. We assume that there exists a unique N ∈ N * such that lim β→∞ µ ∞ (β, ρ 0 ) = E F (ρ 0 ) ∈ (min E N , max E N ), which means that the Fermi energy lies inside the N th Bloch band E N . We also assume that the Fermi surface defined by S F := {k ∈ Ω * : E N (k) = E F (ρ 0 )} is smooth and non-degenerate.

Recall that E N (k) is supposed to be non degenerate outside a (possibly empty) zero Lebesgue measure set of k-points. Our assumption leads to the following consequence:

dist E F (ρ 0 ), ∪ N -1 j=1 E j = d 1 > 0, dist E F (ρ 0 ), ∪ ∞ j=N +1 E j = d 2 > 0. (4.5)
Note that the minimum of the lowest Bloch band E 1 is always simple. If the density ρ 0 is small enough then Assumption 4.1 is automatically satisfied since the Bloch energy function k → E 1 (k) is non-degenerate in a neighborhood of k = 0 (see e.g. [START_REF] Kirsch | Comparison Theorems for The Gap of Schrödinger Operators[END_REF]).

In fact, the non-degeneracy assumption is indispensable for to use of the regular perturbation theory in order to express the functions defined by (3.24), (3.25) and (3.21) (only in the case where j 1 = N ) with the help of the partial derivatives of E N (•) with respect to the k i -variables, for k in a neighborhood of the Fermi surface:

∂E N (k) ∂k i = πN,N (i; k), i ∈ {1, 2, 3}, (4.6)
∂ 2 E N (k) ∂k 2 i = 1 + 2 ∞ j=1 j =N πj,N (i; k) 2 E N (k) -E j (k) , i ∈ {1, 2, 3}, (4.7)
∂ 2 E N (k) ∂k 1 ∂k 2 = ∞ j=1 j =N 2ℜ πj,N (1; k)π N,j (2; k) E N (k) -E j (k) = ∂ 2 E N (k) ∂k 2 ∂k 1 . (4.8)
Such identities have been studied in [START_REF] Cornean | Faraday effect revisited: sum rules and convergence issues[END_REF]. Note that the above series are absolutely convergent if the potential V is smooth enough ( [START_REF] Cornean | Faraday effect revisited: sum rules and convergence issues[END_REF]). Now using Assumption 4.1, we can group the coefficients corresponding to the third and second order derivatives of f appearing in (3.1). This operation allows us to isolate a Landau-Peierls type contribution (the proof can be found in the appendix of this section):

Ω * dk 3 l=2 ∂ l f ∂ξ l β, µ ∞ ; E N (k) c N,l (k) = Ω * dk ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) 1 3! ∂ 2 E N (k) ∂k 2 1 ∂ 2 E N (k) ∂k 2 2 - 1 3! ∂ 2 E N (k) ∂k 1 ∂k 2 2 + a N,2 (k) , (4.9)
where Ω * ∋ k → a j1,2 (k) are the functions defined in (3.22).

From (3.1) and Proposition 4.2 we get an expansion for the orbital susceptibility at fixed density ρ 0 > 0 and inverse of temperature β > 0: Proposition 4.3. Assume for simplicity that E N is a simple band. For every j 1 ∈ N * there exist four families of functions c j1,l (•) with l ∈ {0, 1, 2, 3}, defined on Ω * outside a set of Lebesgue measure zero, such that the second integrand below is bounded and continuous on Ω * :

X (β, ρ 0 ) = - e c 2 1 12β 1 (2π) 3 • Ω * dk ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) ∂ 2 E N (k) ∂k 2 1 ∂ 2 E N (k) ∂k 2 2 - ∂ 2 E N (k) ∂k 1 ∂k 2 2 -3F N (k) + 6 Ω * dk ∞ j1=1 j1 =N 3 l=2 ∂ l f ∂ξ l β, µ ∞ ; E j1 (k) c j1,l (k) + ∞ j1=1 1 l=0 ∂ l f ∂ξ l β, µ ∞ ; E j1 (k) c j1,l (k) , (4.10)
where by convention (∂ 0 ξ f)(β, µ ∞ ; •) := f(β, µ ∞ ; •) and:

F N (k) := -2a N,2 (k) = ∞ j2=1 j2 =N ∞ j3=1 j3 =N C N,N,j2,j3 (k) + C N,j2,N,j3 (k) + C N,j2,j3,N (k) E j2 (k) -E N (k) E j3 (k) -E N (k) + ∞ j2=1 j2 =N C j2,N,N,N (k) -C N,N,j2,N (k) E j2 (k) -E N (k) 2 . (4.11)
Note that we can use identities provided by the regular perturbation theory in order to express the functions c j1,l (as well as F N ) appearing in (4.10) in terms of derivatives of E j and u j w.r.t. the k-variable. But this formulation will only hold true outside a set of k-points of Lebesgue measure zero, while the formulation involving πi,j 's is more general, physically relevant, providing us with bounded and continuous coefficients on Ω * (see Lemma 3.7). Finally keep in mind that the main goal is the Landau-Peierls formula, and it will turn out that only the factor multiplying the second partial derivative of f will contribute to it.

In order to complete the proof of Theorem 1.2 (ii), it remains to take the limit when β → ∞ in (4.10). Since the Fermi energy lies inside the band E N and it is isolated from all other bands, then using (4.3) and (4.4) we have:

lim β→∞ 1 β Ω * dk ∞ j=1 j =N 3 l=2 ∂ l f ∂ξ l (β, µ ∞ (β, ρ 0 ); E j (k))c j,l (k) = 0 and lim β→∞ - 1 β Ω * dk ∂ 2 f ∂ξ 2 (β, µ ∞ ; E N (k)) ∂ 2 E N (k) ∂k 2 1 ∂ 2 E N (k) ∂k 2 2 - ∂ 2 E N (k) ∂k 1 ∂k 2 2 -3F N (k) = - SF dσ(k) |∇E N (k)| ∂ 2 E N (k) ∂k 2 1 ∂ 2 E N (k) ∂k 2 2 - ∂ 2 E N (k) ∂k 1 ∂k 2 2 -3F N (k)
where S F denotes the Fermi surface. Using these two identities together with (4.2) in (4.10), we obtain (1.16).

Appendix -Proof of Proposition 4.2

Using (3.26) we get:

Ω * dk 3 l=2 ∂ l f ∂ξ l β, µ ∞ ; E N (k) c N,l (k) = Ω * dk ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) a N,2 (k)
+

Ω * dk 3 l=2 ∂ l f ∂ξ l β, µ ∞ ; E N (k) b N,l (k) + ∂ 3 f ∂ξ 3 β, µ ∞ ; E N (k) a N,3 (k) .
Using (4.6) and (4.7), the functions b N,l (•), l ∈ {2, 3}, can be rewritten as:

b N,3 (k) = 1 3! ∂E N (k) ∂k 1 2 + ∂E N (k) ∂k 2 2 , b N,2 (k) = - 1 2! - 1 2
∂ 2 E N (k) ∂k 2 1 -1 - 1 2 ∂ 2 E N (k) ∂k 2 2 -1 -1 = 1 2! 1 2 ∂ 2 E N (k) ∂k 2 1 + ∂ 2 E N (k) ∂k 2 2 . Since E N (• ) ∈ C 2 (R 3 /(2πZ 3
)), a simple integration by parts gives us:

∀ i ∈ {1, 2}, π -π dk i ∂E N (k) ∂k i ∂f 3 ∂ξ 3 β, µ ∞ ; E N (k) ∂E N (k) ∂k i = - π -π dk i ∂f 2 ∂ξ 2 β, µ ∞ ; E N (k) ∂ 2 E N (k) ∂k 2 i (4.12)
whence:

Ω * dk ∂ 3 f ∂ξ 3 β, µ ∞ ; E N (k) b N,3 (k) = - 1 3! Ω * dk ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) ∂ 2 E N (k) ∂k 2 1 + ∂ 2 E N (k) ∂k 2 2
and:

Ω * dk 3 l=2 ∂ l f ∂ξ l β, µ ∞ ; E N (k) b N,l (k) = 1 3! 1 2 Ω * dk ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) ∂ 2 E N (k) ∂k 2 1 + ∂ 2 E N (k) ∂k 2
2 (4.13) On the other hand, using (4.6), (4.7) and (4.8), the function a N,3 (• ) can be rewritten as :

a N,3 (k) = 1 3! ∂E N (k) ∂k 1 2 1 2 1 - ∂ 2 E N (k) ∂k 2 2 + ∂E N (k) ∂k 2 2 1 2 1 - ∂ 2 E N (k) ∂k 2 1 - ∂E N (k) ∂k 1 ∂E N (k) ∂k 2 - ∂ 2 E N (k) ∂k 1 ∂k 2 . (4.14)
Note that by a simple integration by parts:

∀ i = j ∈ {1, 2}, π -π dk j ∂E N (k) ∂k j ∂ 3 f ∂ξ 3 β, µ ∞ ; E N (k) ∂E N (k) ∂k j ∂ 2 E N (k) ∂k 2 i (4.15) = - π -π dk j ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) ∂ ∂k j ∂E N (k) ∂k j ∂ 2 E N (k) ∂k 2 i = - π -π dk j ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) ∂ 2 E N (k) ∂k 2 j ∂ 2 E N (k) ∂k 2 i + ∂E N (k) ∂k j ∂ ∂k j ∂ 2 E N (k) ∂k 2 i = - π -π dk j ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) ∂ 2 E N (k) ∂k 2 j ∂ 2 E N (k) ∂k 2 i + ∂E N (k) ∂k j ∂ ∂k i ∂ 2 E N (k) ∂k j ∂k i .
By virtue of (4.14), using (4.15) and (4.12), we get:

Ω * dk ∂ 3 f ∂ξ 3 β, µ ∞ ; E N (k) a N,3 (k) = 1 3! 1 2 Ω * dk ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) 2 ∂ 2 E N (k) ∂k 2 1 ∂ 2 E N (k) ∂k 2 2 + ∂E N (k) ∂k 1 ∂ ∂k 2 ∂ 2 E N (k) ∂k 1 ∂k 2 + ∂E N (k) ∂k 2 ∂ ∂k 1 ∂ 2 E N (k) ∂k 2 ∂k 1 - ∂ 2 E N (k) ∂k 2 1 - ∂ 2 E N (k) ∂k 2 2 + 1 3! Ω * dk ∂ 3 f ∂ξ 3 β, µ ∞ ; E N (k) ∂E N (k) ∂k 1 ∂E N (k) ∂k 2 ∂ 2 E N (k) ∂k 1 ∂k 2 . (4.16)
Finally, by a last integration by parts:

∀ i = j ∈ {1, 2}, π -π dk j ∂E N (k) ∂k j ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) ∂ ∂k i ∂ 2 E N (k) ∂k j ∂k i = - π -π dk j ∂ ∂k i ∂E N (k) ∂k j ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) ∂ 2 E N (k) ∂k j ∂k i = - π -π dk j ∂ 2 E N (k) ∂k i ∂k j ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) + ∂E N (k) ∂k j ∂E N (k) ∂k i ∂ 3 f ∂ξ 3 β, µ ∞ ; E N (k) ∂ 2 E N (k) ∂k j ∂k i .
Then (4.16) is reduced to:

Ω * dk ∂ 3 f ∂ξ 3 β, µ ∞ ; E N (k) a N,3 (k) = 1 3! 1 2 Ω * dk ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) 2 ∂ 2 E N (k) ∂k 2 1 ∂ 2 E N (k) ∂k 2 2 -2 ∂ 2 E N (k) ∂k 1 ∂k 2 2 - ∂ 2 E N (k) ∂k 2 1 - ∂ 2 E N (k) ∂k 2 2
. (4.17)

By adding (4.13) to (4.17) we get :

Ω * dk 3 l=2 ∂ l f ∂ξ l β, µ ∞ ; E N (k) b N,l (k) + ∂ 3 f ∂ξ 3 β, µ ∞ ; E N (k) a N,3 (k) = 1 3! Ω * dk ∂ 2 f ∂ξ 2 β, µ ∞ ; E N (k) ∂ 2 E N (k) ∂k 2 1 ∂ 2 E N (k) ∂k 2 2 - ∂ 2 E N (k) ∂k 1 ∂k 2 2
and we are done. Note that the proof does not work if E N can touch other bands because we loose regularity. In that case the integration by parts have to be done across a tubular neighborhood of the Fermi surface S F , the price being the apparition of some extra terms. These terms will though disappear in the limit β → ∞ because they will decay exponentially with β.

The Landau-Peierls formula

The aim of this section is to establish an asymptotic expansion of (1.16) in the limit of small densities (ρ 0 → 0). Here we prove the expansion (1.17), of which (1.18) is a particular case which has allready been suggested by T. Kjeldaas and W. Kohn in 1957 [START_REF] Kjeldaas | Theory of the Diamagnetism of Bloch Electrons[END_REF].

Proof of Theorem 1.2 (iii)

Let us recall that E 0 = min k∈Ω * E 1 (k) = E 1 (0), and E 1 (k) is non degenerate near the origin with a positive definite Hessian matrix (see e.g. [START_REF] Kirsch | Comparison Theorems for The Gap of Schrödinger Operators[END_REF]). The same reference insures the existence of the following quadratic expansion of E 1 (k) for k → 0:

E 1 (k) = E 0 + 1 2! k T ∂ 2 E 1 ∂k i ∂k j (0) 1≤i,j≤3 k + O k 4 when k → 0
As the Hessian matrix is symmetric, then up to a change of coordinates this quadratic expansion can be rewritten as :

E 1 (k) = E 0 + 1 2 3 i=1 k 2 i m * i + O k 4 when k → 0 (5.1)
where 1 m * i 1≤i≤3 are the eigenvalues of the inverse effective-mass tensor. Consider the assumption of weak density ρ 0 ∈ (0, 1). In this case the Fermi energy defined by (1.12) lies in the interval (E 0 , max k∈Ω * E 1 (k)). When ρ 0 → 0 it follows that E F (ρ 0 ) converges to E 0 . The k-subset of Ω * where E 0 ≤ E 1 (k) ≤ E F (ρ 0 ) is therefore only localized near the origin.

From (5.1) we get the following asymptotic expansion of E F (ρ 0 ) -E 0 when ρ 0 → 0 (the proof is given in the appendix of this section): Proposition 5.1. When ρ 0 → 0, we have the following expansion: .

E F (ρ 0 ) -E 0 = sρ
(5.2)

In the particular case when m * i = m * > 0 for i ∈ {1, 2, 3} and by setting k F := (6π 2 ρ 0 )

1 3 : E F (ρ 0 ) -E 0 = 1 2m * k 2 F + O k 4 F . (5.3)
Before proving Theorem 1.2 (iii), we need one more technical result (its proof is also in the appendix of this section): where s is the coefficient defined in (5.2). Now we are ready to prove the Landau-Peierls formula in (1.17). For this, consider the formula (1.16). Remember that E 1 (•) is non-degenerate and analytic in a neighborhood of the origin. Let us concentrate ourselves on the first term appearing in (1.16):

- e c 2 1 12 1 (2π) 3 SF dσ(k) ∇E 1 (k) ∂ 2 E 1 (k) ∂k 2 1 ∂ 2 E 1 (k) ∂k 2 2 - ∂ 2 E 1 (k) ∂k 1 ∂k 2 2 -3F 1 (k) , (5.6)
since only this term will have a nonzero contribution to the leading term in (1.17). The other term will go to zero like ρ 0 ; this can be shown using (5.5), (5.1), and the fact that the coefficients c 1,1 and c 1,0 are continuous near 0 (see Lemma 3.7). Now consider the following function:

F (k) := ∂ 2 E 1 (k) ∂k 2 1 ∂ 2 E 1 (k) ∂k 2 2 - ∂ 2 E 1 (k) ∂k 1 ∂k 2 2 -3F 1 (k).
By taking into account that F 1 (• ) = -2a 1,2 (• ) (see (4.11)) and by virtue of Lemma 3.7, F (•) is continuous near the origin. According to (5.4), the only thing we need to do is to compute F (0). The determinant of the Hessian matrix gives after a short computation:

∂ 2 E 1 (k) ∂k 2 1 ∂ 2 E 1 (k) ∂k 2 2 - ∂ 2 E 1 (k) ∂k 1 ∂k 2 2 = 1 m * 1 m * 2 + O k 2 when k → 0. (5.7)
Thus we can write: The only thing we have left to do, is proving that F 1 (0) = 0. The definition of F 1 can be found in (4.11), while the coefficients entering in its definition are defined in (3.29).

X M (ρ 0 ) = - e c 2 1 24π 2 m * 1 m * 2 m * 3 1 3 1 m * 1 m * 2 -3F 1 (0) (6π)
Let us start by showing that for all integers j 2 , j 3 ≥ 2 we have:

C 1,1,j2,j3 (0) = C 1,j2,j3,1 (0) = C j2,1,1,1 (0) = C 1,1,j2,1 (0) = 0.

Indeed, in the expression of each of these functions it is possible to identify a factor of the type π1,1 (α; 0), α ∈ {1, 2} which are nothing but partial derivatives of E 1 at the origin, thus they must be zero. It follows that:

F 1 (0) = ∞ j2=2 ∞ j3=2
C 1,j2,1,j3 (0) E j2 (0) -E 1 (0) E j3 (0) -E 1 (0) .

(5.8)

Since:

C 1,j2,1,j3 (0) = π1,j2 (1; 0)π j2,1 (2; 0)π 1,j3 (2; 0)π j3,1 (1; 0) + π1,j2 (2; 0)π j2,1 (1; 0)π 1,j3 (1; 0)π j3,1 (2; 0) -π1,j2 (2; 0)π j2,1 (1; 0)π 1,j3 (2; 0)π j3,1 (1; 0) -π1,j2 (1; 0)π j2,1 (2; 0)π 1,j3 (1; 0)π j3,1 (2; 0), then (5.8) can be rewritten as :

F 1 (0) = 2 ∞ j=2
π1,j (2; 0)π j,1 (1; 0)

E j (0) -E 1 (0) 2 - ∞ j=2
π1,j (2; 0)π j,1 (1; 0)

E j (0) -E 1 (0) 2 - ∞ j=2
π1,j (2; 0)π j,1 (1; 0) E j (0) -E 1 (0) 2 .

(5.9) But for k = 0 we may choose all our eigenfunctions u l (• ; 0) to be real. It means that for all integers j ≥ 2 and α ∈ {1, 2}, the matrix elements π1,j (α; 0) are purely imaginary. As a result, the sums in (5.9) are real numbers and cancel each other, thus F 1 (0) = 0.

Appendix -Proofs of intermediate results

Here we prove Proposition 5.1 and Lemma 5.

Proposition 4 . 2 .

 42 Assume for simplicity that E N is a simple band. Let Ω * ∋ k → c N,2 (k) and Ω * ∋ k → c N,3 (k) the functions respectively defined by (3.28) and (3.27) with j 1 = N . Then:

Lemma 5 . 2 .

 52 Assume that E 1 (k) remains non-degenerate on the ball B ǫ0 (0) := {k ∈ Ω * : |k| ≤ ǫ 0 } with ǫ 0 > 0 small enough. Consider any continuous function F : B ǫ0 (0) → C. Then when ρ 0 → 0 we have the following asymptotic expansions:χ [E0,EF (ρ0)] E 1 (k) F (k) = Bρ 0 + o ρ 0 with B := m *

2 . 3 + 2 r 2 +dr r 2 . 3 .

 232223 Proof of Proposition 5.1. In (5.1) use the change of variables ki := ki √ m * i , with i ∈ {1, 2, 3}.This gives:Ẽ1 ( k) := E 1 ( m * i k) = E 0 O( k4 ). In spherical coordinates: Ẽ1 (r, θ, φ) = E 0 + 1 O(r 4 ) when r → 0.(5.10)We would like to express r as a function of Ẽ1 , θ and φ. Clearly, the equation Ẽ1 r(θ, φ), θ, φ = E 0 + ∆ has a unique solution r(θ, φ, ∆) if ∆ > 0 is small enough. This solution obeys a fixed point equation of the type:r(θ, φ, ∆) = √ 2∆ 1 + O r 2 (θ, φ, ∆) (5.11)which leads to the estimate:r(θ, φ, ∆) = √ 2∆ 1 + O ∆ when ∆ → 0.(5.12)We can finally determine ∆ (thus the Fermi energy) as a function ρ 0 . By setting Ω * := This is the equation we have to solve in order to find ∆ as a function of ρ 0 . Then by standard fixed point arguments we arrive at the estimate (5.2) and we are done.Proof of Lemma 5.2. We only prove (5.4), the other estimate being similar. As before, we prefer the new variables ki = ki √ m * i where i ∈ {1, 2, 3}. Denote by Ẽ1 ( k) = E 1 (k), by F ( k) = F (kThen we can formally write:SF dσ(k) ∇E 1 (k) F (k) = Ω * dk δ E F (ρ 0 ) -E 1 (k) F (k) = m * 1 m * 2 m * 3 Ω * d k δ E F (ρ 0 ) -Ẽ1 ( k) F ( k) = m * 1 m * 2 m *3 { k∈ Ω * s.t. Ẽ1( k)=EF (ρ0) } dσ( k) |∇ k Ẽ1 ( k)| F ( k). (5.13)

  .14) Remark 1 . In other words, a semiconductor/semimetal either has its Fermi energy in the middle of a non-trivial gap (this occurs if max E N < min E N +1 ), or where the two consecutive Bloch bands touch each other closing the gap (this occurs if max E N = min E N +1 ). As for a metal, its Fermi energy lies in the interior of a Bloch band. Remark 2 . According to the above result, E F is discontinuous at all values of ρ 0 for which the equation n ∞ (E) = ρ 0 does not have a unique solution. Each open gap gives such a discontinuity.

	Now here is our main result concerning the orbital susceptibility of a Bloch electrons gas at
	fixed density and zero temperature:
	Theorem 1.2. Denote by E 0

The quadratic expansion (5.1)

Using spherical coordinates, let us denote as before by r(θ, φ, ρ 0 ) the unique root of the equation Ẽ1 r(θ, φ, ρ 0 ), θ, φ = E F (ρ 0 ). Then (5.14) can be rewritten as:

Now by setting ∆ := E F (ρ 0 ) -E 0 and by using (5.12) when ∆ → 0:

Finally, we use (5.2) and the proof is over.

Acknowledgments

H.C. was partially supported by the Danish FNU grant Mathematical Physics. B.S. thanks Cyril Levy for fruitful discussions.


Maryvonne Chalony

Romain Pierrat

Dominique Delande

David Wilkowski

David Wilkowski Coherent

Coherent flash of light emitted by a cold atomic cloud

Keywords:

When sent on a turbid object, a coherent light beam is scattered by the heteoregeneities of the object. If the mean free path ℓ of light inside the scattering object is shorter than its thickness L, one enters into the multiple scattering regime where the photons follow a random walk. Nevertheless, a small part of the incoming photons cross the medium ballistically producing a coherent transmission equal to e -b , where b = L/ℓ is the optical thickness of the medium. When b is large, it seems thus hopeless to coherently transmit a laser beam through an optically thick medium. In this letter, we show experimentally that this is not true and that coherent transmission almost equal to 100% can be achieved, although only over a short temporal window [see e.g. fig. 1b].

At the microscopic level, the coherently transmitted field can be seen as the result of the destructive interference between the incoming laser beam and the field radiated in the forward direction by the dipoles induced in the medium. In the stationary regime, the depletion of the coherent beam is exactly compensated by the intensity scattered in other directions. The simple idea of our experiment is to abruptly switch off a monochromatic laser beam incoming on an optically thick medium. If the laser extinction is fast enough, the induced dipoles continue to radiate a coherent field, creating a coherent flash of light decaying over their lifetime.

The extremely fast response time (fs or ps) of standard dielectric scatterers make the observation of the coherent flash challenging. This issue can be in principle overcome using resonant scatterers with long dwell time. For this reason, the first reported observation of coherent emission after the source extinction -known as free induction decay (FID)-came from NMR few decades ago [1]. In the optical domain, observation of FID in an optically thin medium was reported on a molecular thermal gas [2,3]. More recently, optical FID in a cold cloud of rubidium gas, without inhomogeneous broadening, was also observed and compared with Maxwell-Bloch equations in the zero temperature limit [4][5][6][7]. In this letter, FID in an optically thick medium is studied, with emphasis on the effect of the optical thickness and temperature.

A cold atomic sample of Strontium with up to 10 7 atoms in a volume of about 0.01 mm 3 , is produced in a magneto-optical trap on the 1 S 0 → 3 P 1 intercombination line at λ = 689 nm with a excited lifetime of Γ -1 = 21µs [8]. The minimal temperature of the sample is 0.7 µK which corresponds to kv ≃ 1.6Γ, where v is the mean atomic velocity (along any direction). Thus, in contrast to standard laser cooling, Doppler effect still broadens the atomic transition. After loading, the magnetooptical trap is switched off. A bias B-field of 1 G is applied and a probe laser at exact resonance with the 1 S 0 → 3 P 1 , m = 0 is turned on. The waist of the probe laser is 0.5 mm, so it can be well approximated by a plane wave. The probe is kept on during 40 µs at an intensity lower than half the saturation intensity of the transition (3 µW/cm 2 ). On average, each atom scatters less than 0.3 probe photons, ensuring that mechanical effects are very small. The ignition and the extinction of the probe laser are achieved in less that 50 ns. The forward transmitted light is collected on an intensified CCD, placed on an image plane of the cloud, with a spatial resolution of 20 µm. The gate time is 400 ns. The imaging system is well adapted for spatio-temporal studies. Moreover, the collection solid angle, 5 × 10 -3 , is small enough for the incoherent fluorescence detected by the camera to be negligible. Some keys results are presented in fig. 1. In (a), the images show the shadow of the atomic cloud in the stationary regime (upper image), due to the depletion of the incoming laser beam by the atomic scatterers, and the flash of light emitted by the medium after the laser extinction (lower images). Note the ring structure, discussed below. flash immediately after extinction is almost equal to the incoming laser intensity while the coherent transmission in the stationary regime (just before switch off) was negligibly small. Fig. 1(c) shows the incoherent fluorescence signal, recorded at a angle (30 • ) from the forward direction -so that the coherent flash does not contribute -after the laser extinction. The decay is exponential with a characteristic time τ fluo = 57 µs = 2.7Γ -1 . As a single scattering event is characterized by a dwell time Γ -1 , this proves that multiple scattering and radiation trapping are at play [9,10]. The rather small value τ fluo suggests that multiple scattering remains limited. This is due to the "large" Doppler effect that can bring the scattered photon out of resonance, so that it can escape the medium before being rescattered. We performed Monte Carlo simulations [11] which include both Doppler and recoil effects which turn out to be in excellent agreement with the observed decay time.

In sharp contrast, typical transient times for the coherent flash, see fig. 1(b), are much shorter. It is a clear signature that it has a different physical origin than radiation trapping, and that dephasing phenomena occur at a time scale shorter than τ fluo . In the sequel of this paper, we analyze the physical processes responsible for the loss of phase coherence and show that they explain quan-titatively the experimental observations. A first process is the finite lifetime of the atomic excited state related to the finite width of the atomic resonance. A second effect is due to the residual atomic motion. Once the laser is turned off, the atomic motion leads to dephasing of the dipoles, with a time constant expected to be proportional to (kv) -1 , the time needed for an atom to travel one laser wavelength. A third identified dephasing mechanism occurs in a thick medium: deep inside the medium, the atomic dipoles are driven by the superposition of the incident field and the field radiated by other atoms. Thus, the phase coherence between the various layers is more fragile, producing a shorter coherent flash.

We now sketch the principle of the theoretical calculation. For the sake of simplicity, we assume that the laser beam is a plane wave with wavevector k perpendicular to a infinite slab of scattering medium of thickness L. Since the experiment is performed on an almost perfect twolevel atomic ensemble, the vectorial nature of the fields -i.e. the light polarization -is irrelevant. Note that we do not attempt here to describe the field (multiply) scattered by the medium in other modes; the linearity of the equations implies that it will not affect the coherent outgoing beam. In the stationary regime where the incoming laser is at fixed frequency ω, the outgoing electric field of the laser is proportional to the incoming field:

Ẽout (ω) = Ẽin (ω) exp i n(ω)ωL c (1)
where n(ω) is the complex index of refraction and c the light velocity. When the medium is diluted, the index of refraction is related to the atomic density ρ and the individual atomic polarizability α(ω): n(ω) = 1 + ρα(ω)/2 [START_REF] Hecht | Optics[END_REF]. The scattering mean free path ℓ(ω) is directly related to the imaginary part of the polarizability: ℓ = c/[ρωℑ(α(ω))], and the intensity transmission is simply e -b , with b = L/ℓ(ω) the optical thickness of the medium.

From now on, we consider the case of resonant scatterers -typically isolated atomic scatterers -which present a narrow resonance centered around ω 0 with a width Γ ≪ ω 0 , such that the atomic polarizability is proportional to 1/(ωω 0 + iΓ/2) [START_REF] Loudon | The quantum theory of light[END_REF]. If the laser beam is switched off much faster than Γ -1 , the atomic dipoles can be considered as frozen during the switching, implying that they will continue to radiate the same field immediately after the switching off, slowing decaying over a typical time scale Γ -1 , thus creating the coherent flash of light. The preceding expression for the atomic polarizability is valid for atoms at rest. For moving atoms, it is described by the convolution of the atomic polarizability shifted by Doppler effect 1/(ωkvω 0 + iΓ/2) (with v the atomic velocity along the laser axis) with the velocity distribution.

In order to analyze the transient phenomena, we make a Fourier decomposition of the incoming and outgoing In this example, the electric field Eout is not in phase with the ingoing field Ein. For a resonant laser, the theoretical prediction, essentially confirmed by the experimental observations, see fig. 4, is that the phase vanishes. electric fields. In particular we consider two complementary experiments: in the first one ("off" experiment), a resonant laser beam is sent on the medium during a long time, then abruptly switched off, while in the "on" experiment, the same beam is simply switched on abruptly and remains on forever [START_REF] Ducloy | [END_REF] (see a schematic description in fig. 2). The corresponding Fourier components of the incoming fields are simply:

Ẽoff in (ω) = iE 0 /2π ω -ω 0 + i0 + Ẽon in (ω) = iE 0 /2π ω -ω 0 -i0 + (2)
The temporal dependence of the outgoing electric field is readily obtained by injecting eqs. (2) in the stationary response, eq. ( 1), and Fourier transforming back from frequency to time. The superposition of the "on" and "off" incoming fields is a monochromatic field at frequency ω 0 , which produces a stationary coherent outgoing field E out at the same frequency. The coherently transmitted fields in the "on" and "off" experiments are thus such that:

E on out (t) + E off out (t) = E out (3)
When the laser is switched on abruptly, it takes some time for the atomic dipoles build up and radiate a coherent field antagonist to the incoming one, meaning that the medium is initially transparent, reaching its stationary transmission over a time of the order of Γ -1 . This phenomenon -known as optical nutation [15] -is also schematically described in fig. 2. For large optical thickness, E out is negligibly small compared to the incoming field, so that E off out (t = 0 + ) = E out -E on out (t = 0 + ) ≈ -E in : the coherent flash of light outgoing from the medium is then as intense as the incoming laser (with a π relative phase), although the latter is already switched off!

The detailed calculation of these various effects is possible, and will be presented in a forthcoming publication [16]. The imaginary part of the atomic polarizability is unchanged when the detuning changes sign while its real part changes sign. This implies that Ẽout (ω 0 + δ) ≈ Ẽ * out (ω 0 -δ) and that, for a laser frequency resonant with the atomic frequency, both E on out (t) and E off out (t) are real quantities, in phase (or in opposition) with the incoming field E in , as experimentally observed, see below. The shape of the coherent flash of light depends on the various parameters, but it turns out that, at resonance, the initial decay rate (at time t = 0 + ) has a simple expression:

τ = I off out (t = 0 + ) - dI off out dt (t = 0 + ) = 1 -exp(-b/2) b/2 g kv Γ Γ -1 (4)
where g is related to the complementary error function [START_REF] Abramowitz | Handbook of mathematical functions[END_REF] (assuming a Boltzman velocity distribution):

g(x) = π 8 1 x exp 1 8x 2 Erfc 1 √ 8x (5)
The first term in eq. ( 4) expresses the shortening of the coherent flash in optically thick media, while the g function (always smaller than or equal to unity) the effect of the atomic velocity. For large atomic velocity, g(x) ∝ 1/x, so that the decay time scales as 1/kv, as expected from the width of the atomic absorption profile. These predictions are confronted to the experimental observations in fig. 3: the agreement is good, both when the atomic velocity or the optical thickness are varied. Note that the characteristic decay time of the experimentally observed coherent flash is measured by fitting it with a decreasing exponential, although the shape decreases faster than an exponential. This is why expression (4) slightly overestimates the decay time. When the same fitting procedure is used on the theoretical curves, the characteristic time is typically 20% smaller and agrees perfectly well with the experimental observations. The experimental decay times, shown in fig. 3, are extracted assuming a relative phase φ = 0 between E on out (t) and E in as suggested by the theoretical calculation. However, using the complementarity of the "on" and "off" experiments, one can extract φ(t) using the expression (3), from the measurement of the intensity in both the "off" and the "on" experiment. We have checked, see fig. 4, that φ(t) remains small at short time (t < 2 µs), up to the typical decay time. At longer time, the phase shift is more important especially at large optical thickness. This phase shift, not predicted by the theory developed here, could be due to some experimental imperfections such as a systematic frequency offset of the probe laser with respect to the atomic resonance, or to some inherent limitations of the experiment like mechanical effect when the probe laser is on. It could also come from a breakdown of the independent scattering assumption since the medium is not extremely dilute (ρλ 3 ≃ 0.1). This interesting point is left for future investigations.

In the actual experiment, the atomic medium is not an infinite slab, but a spherical cloud with a Gaussian FIG. 3: (color online) Characteristic decay time of the coherent transmission at the laser ignition (blue symbols) and extinction (red symbols). The dashed line is the analytic prediction, eq. ( 4), based on the short time decay of the coherent pulse, while the solid curve corresponds to a fit of the full theoretical curve with a decreasing exponential. shaped atomic density. When we look at the coherent beam, only forward propagation is important so that the relevant parameter is the optical thickness along an optical ray. It is maximum when the ray passes through the center of the atomic cloud and decays like b(r) = b max exp(-r 2 /2R 2 ) as a function of the transverse distance r from the center (R being the radius of the cloud). The coherent flash is thus expected to decay faster at the cloud center -the larger the optical thickness, the shorter the decay time -than in the external layers. This is visible in the experimental records, producing a spatio-temporal dynamics with a characteristic ring shape visible in fig. 1(a) at time t = 2.4 µs. At large optical thickness, the broadening and the distortion -with respect to the natural Lorentzian shape -of the spectral absorption window [see eq. ( 1)] leads to an oscillatory temporal evolution of the flash, clearly observed in the experiment, see fig. 1(b), with a minimum of the intensity in the center at t = 2.4 µs, followed by a revival at t = 4.8 µs. This reenforces the ring shape. We have also performed ab initio calculations of the electromagnetic field coherently transmitted through a set of randomly placed atomic dipoles [START_REF] Lax | [END_REF] -using the full geometry of the atomic cloud -which confirm this interpretation.

To summarize, thanks to the slow response time of the Strontium intercombination line, we have studied the spatio-temporal dynamics of the coherent transmission of a resonant laser beam across a scattering medium. We have shown the existence of a strong coherent flash of light following the extinction of the laser, and measured its characteristic properties vs. temperature and optical thickness of the medium. Using the complementary information of the laser ignition, we have reconstructed the phase of the forward transmitted electric field.

M. Ducloy and C. Miniatura are acknowledged for fruitful discussions.

Fig. 1 (

 1 b) shows the temporal evolution of the coherent transmission at the cloud center after laser ignition and extinction. The intensity of the coherent

FIG

  FIG. online) (a): Images of the signal transmitted by a cold cloud with maximum optical thickness (at the center) b = 6.5 and temperature T = 3.8 µK. The images are taken downstream during the laser extinction. The first image (t = -0.4 µs) shows, in the stationary regime, the shadow of the atomic cloud which has scattered most of the incoming beam. After the incoming laser is switched off at time t = 0, a coherent flash of light -with a spatio-temporal structure -is emitted by the atomic cloud, decaying over few µs. (b): Coherent transmission at the laser ignition and extinction at the cloud center. (c): Temporal evolution of the incoherent fluorescence intensity (black line) at the laser extinction. The red curve is a fit by a decreasing exponential.

FIG. 2 :

 2 FIG. 2: Schematic view of the coherent transmission. (a) & (b): Two complementary "off" and "on" experiments, where a monochromatic laser is abruptly switched off or on. (c) Schematic temporal evolution of the outgoing electric field in the complex plane, for the two complementary experiments.In this example, the electric field Eout is not in phase with the ingoing field Ein. For a resonant laser, the theoretical prediction, essentially confirmed by the experimental observations, see fig.4, is that the phase vanishes.

  FIG.3:(color online) Characteristic decay time of the coherent transmission at the laser ignition (blue symbols) and extinction (red symbols). The dashed line is the analytic prediction, eq. (4), based on the short time decay of the coherent pulse, while the solid curve corresponds to a fit of the full theoretical curve with a decreasing exponential. (a): Data taken as a function of the average atomic velocity, at approximately constant optical thickness b = 1.2 -1.6. Theoretical curves are computed for b = 1.5. (b): Data taken as a function of the optical thickness b at two temperatures: T = 1.0(2) µK or kv = 1.8Γ (stars) and T = 3.8(4) µK or kv = 3.7Γ (triangles).