Eric Cancès

Clotilde Fermanian Kammerer

Antoine Levitt

Sami Siraj-Dine

Coherent electronic transport in periodic crystals

We consider independent electrons in a periodic crystal in their ground state, and turn on a uniform electric field at some prescribed time. We rigorously define the current per unit volume and study its properties using both linear response and adiabatic theory. Our results provide a unified framework for various phenomena such as the quantization of Hall conductivity of insulators with broken time-reversibility, the ballistic regime of electrons in metals, Bloch oscillations in the long-time response of metals, and the static conductivity of graphene. We identify explicitly the regime in which each holds.

Introduction

We consider a d-dimensional perfect crystal (d = 1, 2 or 3) with periodic lattice R and assume that its electronic structure can be described by an effective linear Hamiltonian H acting on some Hilbert space H. We focus here on the case of spinless continuous models, for which H = L 2 (R d ; C) and

H = 1 2 (-i∇ + A) 2 + V, (1)
where A ∈ L 4 per (R d ; R d ) and V ∈ L 2 per (R d ; R) are R-periodic functions. We adopt the Coulomb gauge choice, ∇ • A = 0 in the sense of distributions. At zero temperature, the ground-state density matrix is given by γ

(0) = 1(H ≤ µ F ), (2)
where µ F ∈ R is the Fermi level, chosen to have a prescribed number of electrons per unit cell.

Depending on the position of µ F in the spectrum σ(H) of H, this can model different types of physical systems. If µ F / ∈ σ(H), the system is an insulator. If µ F is an interior point of σ(H), the system is a metal, or a semi-metal, depending on the density of states of H at µ F . We refer to Section 2.4 for the precise hypotheses we use in each case.

The vector potential A is chosen to be periodic, which excludes the case of a uniform external magnetic field. Our analysis therefore does not directly cover the quantum Hall effect, but can be adapted to do so (see Remark 2.3). It is directly applicable to the quantum anomalous Hall effect [START_REF] Haldane | Model for a quantum Hall effect without landau levels: Condensedmatter realization of the" parity anomaly[END_REF]. We perform our analysis with this particular Hamiltonian, but it can easily be extended to spin-dependent continuous models, tight-binding models, or 2D materials such as graphene (for which the physical space is three-dimensional while the periodic lattice is two-dimensional); see Remark 2.3 for the exact structure needed.

The purpose of this article is to analyze mathematically the behavior of the electrical current appearing in the crystal when a uniform external electric field is turned on instantaneously at the initial time t = 0. In the case of a uniform stationary electric field of magnitude ε > 0 along a (not necessarily normalized) vector e β ∈ R d , the Hamiltonian of the system at time t > 0 is

H ε β = H + εx β , (3)
where x β = x • e β . This operator is self-adjoint on L 2 (R d ; C) (see Proposition 2.2 below), and therefore gives rise to a unitary group (e -itH ε β ) t∈R on L 2 (R d ; C). The electronic state of the system at time t ≥ 0 then is γ ε β (t) = e -itH ε β γ(0)e itH ε β .

The electrical current in the e α -direction at time t ≥ 0 is defined as

j ε α,β (t) = Tr(J α γ ε β (t)), ( 5
)
where Tr is the trace per unit volume (which will be precisely defined in Section 2.1) and J α the current operator along the vector e α ∈ R d (not necessarily normalized nor orthogonal to e β ), defined as

J α = -(-i∇ + A) • e α . (6)
Remark 1.1 (on units and sign convention). If a spinless particle with mass m and charge q is subjected to a electromagnetic field generated by a vector potential A and a scalar potential -εx β generated by a uniform electric field εe β , its Hamiltonian in atomic units is H = 1 2m (-i∇ -qA) 2 qεx β and the charge current operator is J = q(-i∇ -qA). In our definitions (1), ( 3) and ( 6), we have set m = 1 and q = -1 (atomic units) which are the physical values for the electron: this corresponds to applying a force in the direction -e β to the electrons, and measuring their velocity in the direction -e α .

In the limit of weak external fields (ε 1), the qualitative properties of the function t → j ε α,β (t) heavily depends on the physical nature of the material (insulator, metal, semi-metal), as well as on the regime (short, intermediate or long times). Our main results, stated in Theorems 2.7, 2.8 and 2.10, show that the behavior is as follows (see Figure 2 in Section 3)

• For insulators, the time-averaged conductivity

σ α,β = lim t→∞ 1 t ˆt 0 lim ε→0 j ε α,β (t ) ε dt (7)
has a finite value, which is zero in longitudinal directions, and, for 2D materials, is proportional to the Chern number in the transverse direction (quantum anomalous Hall effect).

• For metals, when t ε -1 , the electrons are in the ballistic regime, and the current increases linearly: j ε α,β (t) ≈ D α,β εt. Under some additional assumptions on the Bloch bands, the current displays Bloch oscillations of order 1 when ε -1 t ε -1 log(ε -ζ ) for some small enough ζ > 0.

• For time-reversible 2D semimetals such as graphene, the time-averaged conductivity σ α,β defined in [START_REF] Cornwell | Group Theory in Physics[END_REF] has a finite value equal to 1 16 e α • e β times the number of Dirac points in the Brillouin zone.

Although our formalism is different, our results for insulators and metals are formally consistent with those obtained using the semiclassical equations of motion ẋ = ∇λ n,k , k = -∇V + ẋ×(∇×A) and their higher-order refinements in the case when the n th band is isolated, where the λ n,k 's are the Bloch eigenvalues of H (see Section 2.1). We refer to [START_REF] Panati | Effective dynamics for Bloch electrons: Peierls substitution and beyond[END_REF][START_REF] Stiepan | Semiclassical approximations for hamiltonians with operator-valued symbols[END_REF] for a mathematical analysis of the insulating case.

Note that our results use an averaging in time, and we are unable to conclude anything about what would be the naive definition of the conductivity

lim t→∞ lim ε→0 j ε α,β (t) ε . (8)
A form of averaging of time fluctuations is always necessary to infer zero-frequency behavior from step responses in non-dissipative systems, even in the linear case. The easiest way to see this is by the very simple model for the relationship between an input I(t) and an output O(t):

i Ȯ(t) = ωO(t) + I(t).

This simplified model describes a forced oscillator with eigenfrequency ω, and arises from the linear response of the time-dependent Schrödinger equation of a two-level system. For a constant input I 0 , there is a steady state solution O 0 = R 0 I 0 , where R 0 = -1 ω is the zero-frequency transfer function of the system. However, since this system is oscillatory, this steady state may never be reached: if I is brutally switched on at time 0 with

I(t) = O(t) = 0, if t ≤ 0, I(t) = I 0 if t > 0, then O(t) = O 0 (1 -e -iωt ) = R 0 I 0 (1 -e -iωt
) and we cannot define R 0 as the limit of O(t)/I 0 when t goes to infinity. However, by averaging we obtain lim

t→∞ 1 t ˆt 0 O(t ) I 0 dt = R 0 .
Another common way of retrieving the value of R 0 is by an adiabatic switching of the electric field I(t) = I 0 e ηt for t ≤ 0, I(t) = I 0 for t > 0 [START_REF] Bouclet | Linear response theory for magnetic Schrödinger operators in disordered media[END_REF][START_REF] De | Linear Response Theory -An Analytic-Algebraic Approach[END_REF]. Another possibility is to represent the relationship between O and I by a convolution with a causal response function R(t): O(t) = (R * I)(t), and define the zero-frequency transfer function as lim η→0 + R(iη), as is often done implicitly in the physics literature. Yet another, more physical, possibility is to use a model with dissipation (in this case i Ȯη (t) + iηO η (t) = ωO η (t) + I(t)), compute the zero-frequency transfer function as the long-time limit of O η (t)/I 0 , and then let the dissipation η tend to zero. A particular variant of this scheme is known as the relaxation time approximation [START_REF] Bellissard | The noncommutative geometry of the quantum Hall effect[END_REF] (the relaxation time being proportional to 1/η). For simple systems, all these methods are equivalent.

Note that the problems in the toy model [START_REF] De | Linear Response Theory -An Analytic-Algebraic Approach[END_REF] are related to the presence of a resonance at ω in the model, i.e. a pole in the Fourier transform of the response function. For our perfect crystal model however, the oscillatory components of the response are integrated over the Brillouin zone of the periodic crystal, which induces an averaging. Therefore, these procedures might not be necessary. Indeed, we observe numerically in simple tight-binding models that the naive limit in (8) seems to be well-defined (see Section 3). Identifying precise conditions on the band structure so that this holds will be the subject of future work.

In the metallic case, the conductivity is either infinity or zero, depending on the definition adopted. Indeed, our results imply that lim

t→∞ 1 t ˆt 0 lim ε→0 j ε α,α (t ) ε dt = +∞ is infinite, because j ε α,α (t) ≈ D αα εt in the regime t ε -1 .
On the other hand, in tight-binding models, a simple argument [START_REF] Bellissard | The noncommutative geometry of the quantum Hall effect[END_REF]Proposition 4] shows that lim

ε→0 1 ε lim t→∞ 1 t ˆt 0 j ε α,α (t ) dt = 0.
These two limits correspond to different regimes. In the first one, the electrons undergo ballistic transport, being uniformly accelerated by the electric field. In the second one, the electrons undergo Bloch oscillations, a phenomenon whereby particles in a periodic potential accelerated by a constant force oscillate rather than propagate, as first noticed by Zener [START_REF] Zener | A theory of the electrical breakdown of solid dielectrics[END_REF].

Of course, our model is extremely simple. We assume that the electrons are at zero temperature and we ignore electron-electron interactions, the reaction of the lattice (electron-phonon interactions), and electron scattering by impurities in the crystal. These collision events play a relatively minor role in insulators at low temperatures, with the quantum Hall effect in particular being very robust to perturbations [START_REF] Bellissard | The noncommutative geometry of the quantum Hall effect[END_REF]. However, they impact significantly the properties of metals. In fact, in the linear response regime (ε 1, t ε -1 ), the current increases linearly, representing ballistic transport (see Theorem 2.8). This increase in the velocity of the electrons physically results in an increased collision rate, which acts as dissipation and eventually limits the current. This results in the finite conductivity observed experimentally in macroscopic physics (Ohm's law). The mathematical understanding of this effect in the mathematical framework considered here is left to future work.

The question of quantum transport in solids has attracted significant attention in the mathematical community, with one of the main drivers being the explanation of Anderson localization on the one hand, and the quantum Hall effect and its relation to topological properties on the other hand [START_REF] Teufel | Non-equilibrium almost-stationary states and linear response for gapped quantum systems[END_REF][START_REF] Monaco | Adiabatic currents for interacting electrons on a lattice[END_REF][START_REF] Fiorenza | Z 2 invariants of topological insulators as geometric obstructions[END_REF][START_REF] Panati | Triviality of Bloch and Bloch-Dirac bundles[END_REF][START_REF] Marcelli | Spin conductance and spin conductivity in topological insulators: Analysis of Kubo-like terms[END_REF]. Other topics of interest include the properties of graphene (see for instance [START_REF] Giuliani | Universality of conductivity in interacting graphene[END_REF]), and mesoscopic transport in the Landauer-Buttiker formalism. Comparatively few works have looked specifically at transport in metals. To the best of our knowledge, the present work is the first to present mathematically rigorous results on insulators, metals and semi-metals in a unified framework.

In our results in the linear response regime, we consider the quantity

j LR α,β (t) = lim ε→0 j ε α,β (t) ε .
and then compute the conductivity in the limit t → ∞. This order of limits is necessary to ensure that the electronic state never deviates significantly from its equilibrium. Considering the opposite limit, i.e. the infinite-time dynamics of j ε α,β (t) at ε fixed, is an extremely hard problem, as it is for any dynamics of non-dissipative systems.

The validity of this linear response approximation to derive conductivities, pioneered by Kubo in [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. I. general theory and simple applications to magnetic and conduction problems[END_REF], is by no means obvious. As was pointed out in [START_REF] Nicolaas | Case against linear response theory[END_REF], the most obvious way to derive it, "microscopic linearity" is not physically relevant: assuming a mean free path length d for the electrons, the validity condition that a free electron of (effective) mass m and charge e is not significantly affected by an electric field E in the time interval [0, t] is t 2 2 eE m d. Taking macroscopic times t and reasonable microscopic values for e and m, this limits fields to microscopic values ( [START_REF] Nicolaas | Case against linear response theory[END_REF] quotes 10 -18 Volt/cm), which is unrealistic in practice. The solution of this paradox is "macroscopic linearity": γ ε β (t) has to be understood not as the state of a single set of electrons, but rather as a thermodynamic ensemble. The effective evolution of γ ε β (t) then involves a dissipative term (coming from electron-phonon, electron-impurity or electron-electron interaction) that tends to restore the density matrix to its equilibrium state. The condition of validity of linear response is then that the driving force is negligible compared to the restoring force, which is usually satisfied in practice [START_REF] Vliet | On van kampen's objections against linear response theory[END_REF]. The mathematical justification of linear response is then to consider a more sophisticated model involving a dissipation strength η, and to perform the van Hove limit t → ∞, η → 0, t 2 η = cst, at ε fixed [START_REF] Vliet | On van kampen's objections against linear response theory[END_REF]. Then ε can be taken to zero, and the results of linear response are recovered. Alternatively, a Drude-type model like the relaxation-time approximation [START_REF] Bellissard | Coherent and dissipative transport in aperiodic solids: An overview[END_REF] can be used, with the similar effect of returning the density matrix to equilibrium. Yet another possibility is to never let the density matrix get out of its equilibrium state by switching on adiabatically the external field from negative infinity as e ηt εx β , and to consider the limit η → 0 first then ε → 0 [START_REF] Bachmann | The adiabatic theorem and linear response theory for extended quantum systems[END_REF]. In our naive model, we do not consider a dissipation term, and therefore simply assume the validity of linear response.

Our method of proof is based on the standard gauge change ψ(x, t) = e iεtx β ψ(x, t) that transforms the constant in time but non-spatially-periodic Hamiltonian

H ε β = H + εx β into the time- dependent Hamiltonian H ε β (t) = 1 2 (-i∇ + A -εe β t) 2 + V .
This Hamiltonian is spatially periodic, and the study of its dynamics can be reduced via Bloch-Floquet theory to that of its fibers

H ε β,k (t) = 1 2 (-i∇ + k + A -εe β t) 2 +
V acting on periodic functions (Section 4), for all values of the pseudo-momentum k ∈ R d . Fiber by fiber, this time-dependent Hamiltonian can then be treated using the tools of time-dependent perturbation theory (Section 5). Since time is scaled by ε, the Hamiltonian can be seen as either a small perturbation of the rest Hamiltonian H for small times (in which case we can use linear response to expand j ε α,β (t) to first order in ε for a fixed t, Proposition 5.7), or as a slow perturbation (in which case the adiabatic theorem allows us to access larger time scales t ≈ 1 ε , Proposition 5.3). For insulators and metals in the short-time regime, both tools are applicable and yield the same result. For metals in the Bloch oscillations regime, only the adiabatic theorem is applicable, and for semimetals, only linear response is applicable due to the gap closing at the Dirac points.

The techniques we use (linear response and adiabatic theory) are not new, nor are our results particularly surprising to experts in the field. Rather, we see the contribution of this paper as unifying in the same framework disparate studies on different systems, as well as providing insights on the current response without any specific regularization technique (such as adiabatic switching or dissipation). Our results on Bloch oscillations also appear to be new in the mathematical literature.

The structure of the paper is as follows. We describe our results in Section 2: we define the current in Proposition 2.2, and study its properties for insulators, metals and semi-metals in Theorems 2.7, 2.8 and 2.10. We illustrate numerically the different behaviors we obtain in each of the three settings in Section 3. We devote Section 4 to preliminaries about the regularity and Bloch decomposition of the current. Section 5 states and proves results in adiabatic and linear response perturbation theory. Sections 6, 7 and 8 are devoted to the proof of our results in the case of insulators, metals and semi-metals. Finally two short Appendices are devoted to technical issues.

2 Main results: electrical current in periodic materials

Notation

In this paper we fix

A ∈ L 4 per (R d ; R d ), V ∈ L 2 per (R d ; R) (see below
for the definition of these spaces), µ F ∈ R, and R is the lattice of the d-dimensional crystal. We fix a (non-necessarily orthonormal) basis (e α ) α=1,...,d of the momentum space R d , and set x α = x • e α , A α = A • e α for α = 1, . . . , d. We denote by R * the dual lattice of the periodic lattice R, by Ω an arbitrary unit cell in the physical space, and by B an arbitrary unit cell in the reciprocal space (which we will call by abuse of language the Brillouin zone). In the special case of a cubic crystal of lattice parameter a > 0, we have R = aZ d , R * = 2π a Z d , and we can take Ω = [0, a) d , B = [-π a , π a ) d . The R-periodic Lebesgue and Sobolev spaces are denoted by

L p per := u ∈ L p loc (R d ; C) | u R-periodic , H s per := u ∈ H s loc (R d ; C) | u R-periodic .
The space of bounded linear operators on a Hilbert space H is denoted by L(H), and the Schatten class of bounded operators A ∈ L(H) such that Tr(|A| p ) < ∞ by S p (H). For R ∈ R, we denote by τ R the translation operator formally defined by τ R φ = φ(•-R). Depending on the context, τ R will be seen as a unitary operator on L 2 (R d ; C), or as a linear operator on some R-translation invariant subspace of

D (R d ; C). A bounded operator on L 2 (R d ; C) is called R-periodic if it commutes with τ R for all R ∈ R. An unbounded self-adjoint operator on L 2 (R d ; C) is called R-periodic if its resolvent is R-periodic. A bounded R-periodic operator A ∈ L(L 2 (R d ; C)) is called locally trace-class if χAχ ∈ S 1 (L 2 (R d ; C)) for any compactly supported function χ ∈ L ∞ (R d ; C). For p ≥ 1,
∀χ ∈ C ∞ c (R d ; C), Tr(Aχ) = ˆRd ρ A χ.
The trace per unit volume of an operator A ∈ S 1,per is defined as

Tr(A) = 1 |Ω| Tr L 2 (R d ;C) (1 Ω A1 Ω ) = Ω ρ A ,
where 1 Ω is the characteristic function of the unit cell Ω, and ffl Ω is a shorthand notation for 1 |Ω| ´Ω. This formula is independent of the choice of the unit cell Ω.

Since we are dealing here with periodic materials, we will use the Bloch transform (also called Bloch-Floquet transform) [START_REF] Michael | Methods of Modern Mathematical Physics. IV. Analysis of operators[END_REF][START_REF] Kuchment | An overview of periodic elliptic operators[END_REF][START_REF] Wilcox | Theory of Bloch waves[END_REF]. For K ∈ R * , let T K be the unitary multiplication operator on L 2 per defined by

∀v ∈ L 2 per , (T K v)(x) = e -iK•x v(x) for a.a. x ∈ R d ,
and

L 2 qp (L 2 per ) := R d k → u k ∈ L 2 per | ˆB u k 2 L 2 per dk < ∞, u k+K = T K u k for all K ∈ R * and a.a. k ∈ R d ,
the Hilbert space of R * -quasi-periodic L 2 per -valued functions on R d endowed with the inner product

u, v L 2 qp (L 2 per ) = B u k , v k L 2 per dk.
Here and below, the subscript qp refers to the quasi-periodicity property. The Bloch transform then is the unitary map from

L 2 (R d ; C) to L 2 qp (L 2 per ) defined for u ∈ C ∞ c (R d ; C) by ∀k ∈ R d , ∀x ∈ R d , u k (x) = R∈R u(x + R) e -ik•(x+R) . ( 10
)
Its inverse is given by

u(x) = B e ik•x u k (x) dk, for a.a. x ∈ R d . (11)
Any R-periodic operator A ∈ L(L 2 (R d ; C)) is decomposed by the Bloch transform in the sense that there exists a function

k → A k in L ∞ qp (L(L 2 per )) such that for any u ∈ L 2 (R d ; C) and almost all k ∈ R d , (Au) k = A k u k , and
A k+K = T K A k T * K , for all K ∈ R * and a.a. k ∈ R d . (12)
The A k 's are called the fibers of the operator A. The Bloch decomposition theorem can be extended to unbounded R-periodic self-adjoint operators using the resolvent [START_REF] Michael | Methods of Modern Mathematical Physics. IV. Analysis of operators[END_REF].

If A ∈ S 1,per , then the function k → A k is in L 1 qp (S 1 (L 2 per )), the function k → Tr(A k ) is in L 1 loc (R d ),
In the case of the periodic Hamiltonian operator H given by (1), we have

H k = 1 2 (-i∇ + k + A) 2 + V. (13)
For each k ∈ R d , H k is a bounded below self-adjoint operator on L 2 per with domain H 2 per and compact resolvent. Let (λ n,k ) n∈N * be the non-decreasing sequence of eigenvalues of

H k counting multiplicities λ 1,k ≤ λ 2,k ≤ λ 3,k ≤ • • • , lim n→∞ λ n,k = +∞,
and we use the convention λ 0,k = -∞. We denote by (u n,k ) n∈N * ∈ (H 2 per ) N * an L 2 per -orthonormal basis of associated eigenfunctions:

H k u n,k = λ n,k u n,k , u m,k , u n,k L 2 per = δ m,n .
For N ∈ N * and k ∈ R d , we will denote by

P N,k = 1(H k ≤ λ N,k ). ( 14
)
Whenever λ N,k < λ N +1,k , P N,k is the spectral projector on the eigenspace associated with the lowest N eigenvalues of H k (counting multiplicities):

P N,k = N n=1 |u n,k u n,k |. ( 15
)
Since H k is quasi-periodic, so is P N,k , and the eigenvalues λ n,k are R * -periodic functions of k. By a min-max argument (see e.g. [START_REF] Michael | Methods of Modern Mathematical Physics. IV. Analysis of operators[END_REF][START_REF] Cancès | Numerical quadrature in the Brillouin zone for periodic Schrödinger operators[END_REF]), there exists C 1 , C 1 ∈ R, and C 2 , C 2 > 0 such that

C 1 + C 2 n 2/d ≤ λ n,k ≤ C 1 + C 2 n 2/d . (16)
Denoting by N k the number of eigenvalues below the Fermi level µ F at k

N k = λ n,k ≤ µ F , n ∈ N * , (17)
we see that N k is bounded uniformly in k.

Let us now consider the ground-state density matrix γ(0) = 1(H ≤ µ F ) defined in [START_REF] Bellissard | Coherent and dissipative transport in aperiodic solids: An overview[END_REF]. Its Bloch fibers are

γ k (0) = 1(H k ≤ µ F ) = P N k ,k . (18)
The current operator J α = -(-i∇ + A) • e α defined in ( 6) is also R-periodic, with fibers

J α,k = -(-i∇ + k + A) • e α = -∇ k H k • e α =: -∂ α H k .
Note that the notation ∂ α denotes a derivative along the (not necessarily normalized) vector e α .

Lastly, for each q ∈ R d , we denote that G q the unitary multiplication operator on

L 2 (R d ; C) defined by ∀u ∈ L 2 (R d ; C), (G q u)(x) = e iq•x u(x) for a.a. x ∈ R d . ( 19
)
The operator G q is not R-periodic, except when q ∈ R * (in which case G q is fibered, with G q,k = T -q for all k). However, for any R-periodic operator A ∈ L(L 2 (R d ; C)) and any q ∈ R d , the operator G q AG * q is R-periodic and its Bloch decomposition is given by

(G q AG * q ) k = A k-q , for a.a. k ∈ R d . (20)

The Bloch theorem

Before attacking the well-posedness of the current j ε α,β (t) = Tr(J α γ ε β (t)) for ε, t = 0, we first study an easier special case. Proposition 2.1 (Bloch theorem). The current satisfies j 0 α,β (t) = 0, ∀t ≥ 0 (no current in the absence of external field), j ε α,β (0) = 0, ∀ε ≥ 0 (continuity of the current at t = 0).

Proof. This is a classical statement going back to Bloch, valid in a more general context. We adapt here the proof in [START_REF] Bohm | Note on a theorem of bloch concerning possible causes of superconductivity[END_REF]. We have

J := j 0 α,β (t) = j ε α,β (0) = Tr(J α γ(0)) = -(2π) -d ˆB Tr(∂ α H k γ k (0)) dk.
Assume that this quantity is non-zero. Construct for δ ∈ R a trial state

γ δ = e -iδeαx γ(0)e iδeαx ,
a periodic operator with fibers

γ δ k = γ k+δeα ,
and compute

Tr(Hγ δ ) = (2π) -d ˆB Tr(H k γ k+δeα (0)) dk = (2π) -d ˆB Tr(H k-δeα γ k (0)) dk = Tr(Hγ(0)) -δ(2π) -d ˆB Tr(∂ α H k γ k (0)) dk + O(δ 2 ).
Since γ(0) is the ground state, Tr(Hγ δ ) ≤ Tr(Hγ(0)) for all δ, and therefore J = 0.

Definition of the current

For ε > 0, the operator

H ε β = H + εx β = 1 2 (-i∇ + A) 2 + V + εx • e β
already introduced in ( 3) is not R-periodic, and we would naively expect that the density matrix

γ ε β (t) = e -itH ε β γ(0)e itH ε β
at time t > 0 (already introduced in (4)) is not either. Yet, this operator is in fact R-periodic.

Physically, this is due to the fact that although the potential V el (x) := εx • e β is not periodic, the field E = -∇V el = -εe β to which the electrons are subjected is constant, hence periodic. The proof of this result relies on the standard gauge transform

ψ(x, t) = (G εte β ψ(•, t) (x) = e iεtx β ψ(x, t), (21)
where the operator G q has been defined in [START_REF] Hagedorn | Molecular propagation through electron energy level crossings[END_REF], and the introduction of the gauge-transformed operators

U ε β (t, t ) := G εte β e -i(t-t )H ε β G * εt e β , (22)
and

γ ε β (t) := G εte β γ ε β (t)G * εte β = U ε β (t)γ(0) U ε β (t) * , (23)
where

U ε β (t) is a short-hand notation for U ε β (t) := U ε β (t, 0) = G εte β e -itH ε β .
Through the change of gauge ( 21), the dynamics induced by the time-independent but non-periodic Hamiltonian H ε β is equivalent to the dynamics induced by the time-dependent periodic Hamiltonian

H ε β (t) = G εte β H ε β G * εte β = 1 2 (-i∇ + A -εe β t) 2 + V. (24)
This change of gauge is standard in both the mathematical and physical literature, as it turns the spatially inhomogeneous electric potential V el = -εx β into a homogeneous (but time-dependent) magnetic potential A el = -εe β t, more convenient to deal with here because it does not break periodicity. Physically, this is a manifestation of the gauge invariance of the Schrödinger equation, where an electric field E = -∇V el -∂A el ∂t can be realized either through a scalar or vector potential. The Bloch fibers of H ε β (t) are

H ε β,k (t) = 1 2 (-i∇ + k + A -εe β t) 2 + V = H k-εe β t . (25)
We sum up these arguments in the proposition below, together with elements that we shall use for defining the current. The reader can refer to the articles [START_REF] Bouclet | Linear response theory for magnetic Schrödinger operators in disordered media[END_REF][START_REF] Leinfelder | Schrödinger operators with singular magnetic vector potentials[END_REF] where part of the results of that Proposition are proved.

Proposition 2.2. Let A ∈ L 4 per (R d ; R d ) such that ∇ • A = 0, and V ∈ L 2 per (R d ; R). 1. For all ε ∈ R, the operator H ε β defined in (3) is essentially self-adjoint on C ∞ c (R d ; C)
, and therefore admits a unitary propagator (e

-itH β ε ) t∈R in L 2 (R d ; C). 2. For all t ∈ R, and ε ∈ R, the operator H ε β (t) defined in (24) is self-adjoint on L 2 (R d ) with do- main H 2 (R d ; C), and R-periodic. The strongly continuous unitary propagator ( U ε β (t, t )) (t,t )∈R×R on L 2 (R d ; C) defined in (22) is R-periodic for all t, t ∈ R, with fibers U ε β,k (t, t ) solving i∂ t U ε β,k (t, t ) = H ε β,k (t) U ε β,k (t, t ), U ε β,k (t , t ) = Id L 2 per . (26)
3. For all t ≥ 0, and

ε ∈ R, J α γ ε β (t) ∈ S 1,per . The current j ε α,β (t) = Tr(J α γ ε β (t)
) is well-defined and

j ε α,β (t) = -(2π) -d ˆB Tr ∂ α H ε β,k (t) γ ε β,k (t) dk (27) = -(2π) -d ˆB Tr ∂ α H k-εe β t U ε β,k (t)γ k (0) U ε β,k (t) * dk (28) γ(0) γ ε β (t) γ(0) γ ε β (t) e -itH ε β Id Gεte β U ε β (t) {γ k (0)} k∈R d = {γ ε β,k (0)} k∈R d {γ ε β,k (t)} k∈R d { γ ε β,k (0)} k∈R d = {γ ε β,k (0)} k∈R d { γ ε β,k (t)} k∈R d = {γ ε β,k-εte β (t)} k∈R d {Id} k∈R d translation in k space by-εte β { U ε β,k (t)} k∈R d
B k = U k A k U * k .
The results of Proposition 2.2 are not new (some are classical) but are nevertheless proved in Section 4 for the sake of completeness. The situation can be summed up in the commutative diagrams of Figure 1.

Remark 2.3. This proposition reduces the study of j ε α,β (t) to that of the dynamics of the timedependent Hamiltonian H k-εe β t . In particular, although we have focused on the specific Hamiltonian H given by ( 1), all computations beyond the proof of this proposition will be based on the use of the three formulae: for all k ∈ R d , t ∈ R + ,

γ k (0) = 1(H k ≤ µ F ), i∂ t U ε β,k (t) = H k-εe β t U ε β,k (t), U ε β,k (0) = Id H f , j ε α,β (t) = -(2π) -d ˆB Tr ∂ α H k-εe β t U ε β,k (t)γ k (0) U ε β,k (t) * dk, (29) (30) (31)
where the fiber H f is equal to L 2 per in our setting. Our results in the following sections can therefore be extended to other Hamiltonians where (H k ) k∈R d is a family of bounded below selfadjoint operators on a Hilbert space H f with compact resolvent satisfying the quasi-periodicity conditions [START_REF] Kammerer | Single switch surface hopping for molecular dynamics with transitions[END_REF]), and the boundedness conditions in Section 5. This includes in particular spin-dependent continuous models, tight-binding lattice models (for which H f = C M ), and 2D materials. It also contains the case of systems with a constant magnetic field where the flux per unit cell satisfies an adequate commensurability condition (see [START_REF] Freund | Peierls substitution for magnetic Bloch bands[END_REF] or [START_REF] Monaco | Optimal decay of Wannier functions in Chern and quantum Hall insulators[END_REF], where ideas from [START_REF] Zak | Magnetic translaion group[END_REF] are implemented).

H k+K = T K H k T * K , ∀K ∈ R * , k ∈ R d , where (T K ) K∈R * is a unitary representation of the group R * on H f (see

Insulators, non-degenerate metals, semimetals

As we said before, the position of the Fermi level in the band diagram (λ n,k ) n∈N * , k∈B is key to determining the electronic properties of the medium. We define the Fermi surface sheets

S n = {k ∈ B | λ n,k = µ F }, n ∈ N *
and the Fermi surface

S = n∈N * S n = {k ∈ B | ∃n ∈ N * s.t. λ n,k = µ F }. ( 32
)
We will be interested here in three types of systems that we now describe in three mutually exclusive assumptions.

Assumption 2.4 (insulator). The Fermi surface S is empty, and there exists

N ins ∈ N * such that N k = N ins for all k ∈ B, i.e. ∀k ∈ B, λ Nins,k < µ F < λ Nins+1,k , or equivalently µ F / ∈ σ(H).
In the case of insulators, we have for all k ∈ R d γ k (0) = P Nins,k , and γ k (0) is a real-analytic R * -quasi-periodic function.

Assumption 2.5 (non-degenerate metal). The Fermi surface S is non-empty and the following conditions are satisfied: for all n ∈ N * ,

• S n ∩ S n+1 = ∅ (no crossing at the Fermi level);

• for all k ∈ S n , ∇λ n,k = 0 (no flat bands at the Fermi level).

Note that this assumption was used in [START_REF] Cancès | Numerical quadrature in the Brillouin zone for periodic Schrödinger operators[END_REF]. It ensures a smooth density of states at the Fermi level. In this case, the Fermi surface consists of a finite union of disjoint smooth closed surfaces S n . Letting

B n = {k ∈ B | λ n,k < µ F < λ n+1,k }, we obtain a partitioning B = S n∈N * B n .
Both N k and the fibers γ k (0) = P N k ,k of the density matrix γ(0) are smooth on each B n , and have discontinuities on the sheets S n . Assumption 2.6 (semimetal). The dimension d is equal to 2, there is N sm such that λ Nsm,k ≤ µ F for all k ∈ B, and the Fermi surface S consists of a finite number of isolated points (k i ) i∈I ("Dirac points"). All these points are conical crossings: for all i ∈ I,

λ Nsm-1,ki < λ Nsm,ki = µ F = λ Nsm+1,ki < λ Nsm+2,ki , (33)
λ Nsm,k = µ F -v F,i |k -k i | + O(|k -k i | 2 ), (34)
λ Nsm+1,k = µ F + v F,i |k -k i | + O(|k -k i | 2 ), (35)
for some v F,i ∈ R. Furthermore, in this case we assume that A = 0, so that the system has the time-reversal symmetry

H -k = H k .
Note that we assumed in Assumption 2.6 that A = 0 to ensure time-reversal symmetry. We require more regularity on V than in the previous assumptions to be able to prove a Dyson expansion for the propagator (see Proposition 5.7). For the sake of clarity, we consider a model of 2D semimetals set in R 2 , but our arguments can be adapted to the more physical case of a model set in R 3 (see also Remark 2.3).

Assumption 2.6 is generic in the case of potentials possessing the symmetry of honeycomb lattices, such as graphene [START_REF] Fefferman | Honeycomb lattice potentials and Dirac points[END_REF]. In this case, there are two non-equivalent Dirac points in the Brillouin zone (|I| = 2), usually denoted by K and K , and we have

K = -K and v F,1 = v F,2 . The constant v F = v F,1 = v F,2
is known as the Fermi velocity. More generally, Dirac points generate specific dynamical behaviors that have been studied in [START_REF] Kammerer | Mesures semi-classiques et croisements de mode[END_REF][START_REF] Kammerer | A kinetic model for the transport of electrons in a graphene layer[END_REF] in the context of the Dirac operator. Such phenomena also appear in molecular dynamics (see [START_REF] Hagedorn | Molecular propagation through electron energy level crossings[END_REF][START_REF] Kammerer | An Egorov theorem for avoided crossings of eigenvalue surfaces[END_REF][START_REF] Kammerer | Single switch surface hopping for molecular dynamics with transitions[END_REF]).

In the semimetal case, N k = N sm for almost every k ∈ R 2 , and γ k (0) is singular at each k i ∈ S.

Main results: the current

In the following results, we use the notation O(f (ε, t)) to denote a quantity bounded by Cf (ε, t) where C is a constant that might depend on the material through V , A and µ F , but not on t and ε. Theorem 2.7 (insulators). Assume the system is an insulator (Assumption 2.4). Then there exists η > 0 such that for all ε, t ∈ R + ,

1 t ˆt 0 j ε α,β (t ) ε dt = -i(2π) -d ˆB Tr (γ k (0)[∂ α γ k (0), ∂ β γ k (0)]) dk + O 1 t + ε(1 + t) e ηεt .
Note that this implies in particular that

σ α,β = lim t→∞ 1 t ˆt 0 lim ε→0 j ε α,β (t ) ε dt = e T α σ ⊥ e β ,
where σ ⊥ is a real antisymmetric matrix with components

σ ⊥ ij := (2π) -d ˆB -i Tr γ k (0) ∂γ k ∂k i (0), ∂γ k ∂k j (0) dk. ( 36
)
The integrand in ( 36) is related to the well-known Berry curvature associated to the first N ins bands, that is to the 2-form

1≤i<j≤d

Ω ij (k) dk i ∧ dk j where Ω ij := -i Tr γ k (0) ∂γ k ∂k i (0), ∂γ k ∂k j (0) For d = 2, we have σ ⊥ 12 = (2π) -1 Ch 1 (γ • (0)), where Ch 1 (γ • (0)) ∈ Z
is the first Chern of the fiber bundle defined by the quasi-periodic function k → γ k (0) [START_REF] David J Thouless | Quantized Hall conductance in a two-dimensional periodic potential[END_REF][START_REF] Simon | Holonomy, the quantum adiabatic theorem, and Berry's phase[END_REF]. This relationship between the transverse bulk transport properties and the Chern number, characteristic of the integer quantum Hall effect, is known as the TKNN formula.

If A = 0, then the system has the time-reversal symmetry H -k = H k . As is classical, the Berry curvature is then odd, and the transverse conductivity matrix σ ⊥ equal to zero [START_REF] David J Thouless | Quantized Hall conductance in a two-dimensional periodic potential[END_REF].

Theorem 2.8 (conductivity in non-degenerate metals). Assume the system is a non-degenerate metal (Assumption 2.5).

1. Let θ > 0. For all ε > 0 small enough and 0 ≤ t ≤ 1 ε ε θ , we have

j ε α,β (t) = D α,β εt + O(ε + ε 2 t 2 ) ( 37
)
where

D α,β := (2π) -d n∈N * ˆSn ∂ α λ n,k (ds • e β ) (38)
2. If furthermore there exists

N met ∈ N * such that λ Nmet-1,k < µ F < λ Nmet+1,k for all k ∈ B
and there are uniform gaps between λ Nmet-1,k and λ Nmet,k on the one hand, and λ Nmet,k and λ Nmet+1,k on the other hand, then there exists η > 0 such that, for all ε, t ∈ R + ,

j ε α,β (t) = -(2π) -d ˆB1(λ Nmet,k ≤ µ F )∂ α λ Nmet,k+εe β t dk + O((ε + ε 2 t)e ηεt ). ( 39
)
Note that under the assumptions of the case 2 above, the lowest N -1 bands are completely filled, the N th band is partially filled, and the other bands are empty. Still in the setup of case 2, it follows from (37) and ( 39) that four different regimes can be observed for ε 1

1. For very short times t 1, quantum fluctuations of order O(ε) dominate the current:

j ε α,β (t) = O(ε); 2. For 1 t 1 ε
, the electrons undergo ballistic transport:

j ε α,β (t) ≈ D α,β εt,
where D α,β is defined in (38);

3. For

1 ε t 1 ε log ε -ζ with ζ ∈ (0, η -1
), we observe Bloch oscillations

j ε α,β (t) ≈ -(2π) -d ˆB 1(λ Nmet,k ≤ µ F )∂ α λ Nmet,k+εe β t dk.
In particular, when e β is commensurate with the reciprocal lattice R * , the current is well approximated in this regime by a periodic function of time with zero mean;

4. for times t 1 ε log ε -ζ , our estimates do not allow us to conclude. The proofs show that the factor e ηεt is due to the unboundedness of the operator H defined in [START_REF] Bachmann | The adiabatic theorem and linear response theory for extended quantum systems[END_REF]. For tight-binding models, this factor e ηεt is not present, and we would observe Bloch oscillations up to times t 1 ε 2 . The behavior for larger times is open. Note that some periodic metallic systems have a more complex crossing structure than that assumed in the second case of Theorem 2.8. This is the case in particular for the free electron gas (V = 0, A = 0, seen as a periodic system with an arbitrary periodic lattice), which does not display Bloch oscillations.

Remark 2.9. The coherent electronic transport model considered here neglects all sources of dissipation (phonons, impurities, electron-electron interactions). In the Drude approximation, these phenomena give rise to an effective timescale τ such that 1 τ 1/ε (larger than the coherence timescale of the electrons, but smaller than the Bloch oscillations timescale), yielding a finite DC conductivity σ α,β ∼ D α,β τ . In usual metals at room temperature, dissipation is dominated by phonon scattering, and the relaxation time τ is of the order of tens of femtoseconds [START_REF] Gall | Electron mean free path in elemental metals[END_REF]. By contrast, the timescale of Bloch oscillations in most experiments is much larger. Only in structures such as semiconductor superlattices or cold atoms have Bloch oscillations been observed experimentally [START_REF] Leo | Observation of Bloch oscillations in a semiconductor superlattice[END_REF].

Theorem 2.10 (conductivity in semi-metals). Assume that the system is a semimetal (Assumption 2.6). Assume furthermore that V ∈ H 1 per . Then,

σ α,β = lim t→∞ 1 t ˆt 0 lim ε→0 j ε α,β (t ) ε dt = |I| 16 e α • e β .
Semimetals are intermediate between insulators and metals, possessing a finite longitudinal conductivity in the linear response regime. This is due to the peculiar properties of the Dirac points. Note that the value of the conductivity is universal, not depending on the characteristics of the Hamiltonian but only on the number of conical crossings. More precisely, the conductivity tensor is isotropic and each conical intersection contributes as 1 16 to the total conductivity. Note that this result is consistent with formula (1.17a) in [START_REF] Das Sarma | Electronic transport in two-dimensional graphene[END_REF].

Numerics

Before turning to the proofs, we illustrate our results with numerical simulations. As mentioned in Remark 2.3, our results also apply to tight-binding models, and only depend on the form of H k . We test on a very simple model of H k , adapted from the Haldane model [START_REF] Haldane | Model for a quantum Hall effect without landau levels: Condensedmatter realization of the" parity anomaly[END_REF] (itself based on a tight-binding model of graphene), that can support many phases depending on the values of its parameters. The graphene lattice R is spanned by the vectors

a 1 = √ 3 2 , 1 2 , a 2 = √ 3 2 , - 1 2 ,
and R * by the vectors

b 1 = 2π 1 √ 3 , 1 , b 2 = 2π 1 √ 3 , -1 .
The Hamiltonian fibers are

H k = m(k) f (k) f (k) -m(k) , with m(k) = g + 2t 2 sin(k • a 1 ) + sin(k • a 2 ) + sin(k • (a 1 -a 2 )) , f (k) = 3 i=1 e ik•δi , δ 1 = 1 √ 3 , 0 , δ 2 = - 1 2 √ 3 , 1 2 , δ 3 = - 1 2 √ 3 , - 1 2 .
The eigenvalues of

H k are λ ± = ± m(k) 2 + |f (k)| 2 .
With g = 0, t 2 = 0, this is the standard model of graphene: two bands touching at level 0 at two inequivalent points in the Brillouin zone, where f (k) vanishes. The parameter g = 0 opens a gap of size 2g. The parameter t 2 models an internal magnetic field, and can turn the system into a Chern insulator (in particular, with g = 1, t 2 = -1, the system is a Chern insulator with Chern number +1). Therefore, varying the parameters g, t 2 and µ F , we can obtain a normal insulator, a Chern insulator, a semimetal or a metal. For a given set of parameters, we compute the current by using formulae ( 29)- [START_REF] Rackauckas | DifferentialEquations.jl-a performant and feature-rich ecosystem for solving differential equations in Julia[END_REF]. We sample the Brillouin zone using a uniform grid with N grid = 300 points per direction, and solve the ordinary differential equation

i du dt (t) = H k-εe β t u(t), u(0) = u n,k ,
for various n and k using the DifferentialEquations.jl Julia package [START_REF] Rackauckas | DifferentialEquations.jl-a performant and feature-rich ecosystem for solving differential equations in Julia[END_REF] with the default Tsitouras method of order 5.

Our parameter values are collected in Table 1.

Panel

g µ F t 2 Phase (a) 1 0 0 Normal insulator (b) 1 0 -1 Chern insulator (c) 1 -2 0 Metal (d) 0 0 0 Semimetal Table 1:
Parameter values for the experiments in Figure 2 Our results are presented in the linear response regime (ε = 10 -6 , t 1 ε ) in Figure 2. These results are consistent with our theoretical results, including the limit values of the conductivity in cases (b) and (d), where we obtain 4π/ √ 3 ≈ 7.26 and |b 1 | 2 /8 = 2π 2 /3 ≈ 6.58 respectively. However, there is an additional phenomenon worth of note: in the case of insulators and graphene, the linear response instantaneous conductivity j α,β (t) = lim ε→0 j ε α,β (t) ε seems to possess a finite limit as t → +∞. This is not captured by our results, where we used an averaging process to suppress the oscillations. Note that for a finite N grid , the linear response oscillates with frequencies λ n ,k -λ n,k for λ n,k < µ F < λ n ,k , and k in the discrete Brillouin zone. Only in the limit N grid → ∞ do these resonances merge together to yield a finite limit for the current. This is linked to the absence of resonances (parallel bands) in our model. A deeper investigation of this effect would be interesting future work.

We also investigate the Bloch oscillations regime ε 1, 1 ε t in Figure 3, where we use the same parameters as in case (c) above. The result is consistent with our theoretical result: periodic or quasi-periodic oscillations, depending on whether e β is commensurate with the reciprocal lattice or not.

4 Bloch decomposition of γ ε β (t) and regularity of the current

In this section, we prove Proposition 2.2. We first point out an alternative interpretation that helps shedding some light on the gauge change G εte β . Formally, γ ε β (t) satisfies the equation

i∂ t γ ε β = [H, γ ε β ] + ε[x β , γ ε β ].
The operator [x β , γ ε β ] can easily be seen to be R-periodic, with fibers i∂ β γ ε β,k (where ∂ β = e β • ∇ k ), and therefore, γ ε β (t) is R-periodic and its fibers γ ε β,k (t) satisfy the equation where

i∂ t γ ε β,k -iε∂ β γ ε β,k = [H k , γ ε β,k ] = L H k γ ε β,k ,
L H k := [H k , •]
is the Liouvillian associated with the operator H k (see Section 5.1). The left-hand side of this equation is a linear advection equation, which suggests the use of the method of characteristics: setting

γ ε β,k (t) = γ ε β,k-εe β t (t), (40)
we obtain

i∂ t γ ε β,k (t) = [H k-εe β t , γ ε β,k (t)] = L H k-εe β t γ ε β,k (t),
which is equivalent to [START_REF] Panati | Effective dynamics for Bloch electrons: Peierls substitution and beyond[END_REF]. The use of the gauge transform operator G εte β , equivalent to the change of variable [START_REF] Vliet | On van kampen's objections against linear response theory[END_REF], makes these remarks rigorous.

We now prove Proposition 2.2. As outlined above, the results of this proposition are wellknown; they can in fact be extended to the more general setting of ergodic magnetic Schrödinger operators (see [START_REF] Bouclet | Linear response theory for magnetic Schrödinger operators in disordered media[END_REF]). We provide here an elementary proof specific to the periodic case, and take this opportunity to introduce notations and tools which will be useful in the sequel.

Proof of the first assertion. The essential self-adjointness of H ε β follows from an extension of the Faris-Lavine theorem [START_REF] Michael | Methods of Modern Mathematical Physics. II. Fourier analysis, self-adjointness[END_REF]Theorem X.38]. Let C = C ∞ c (R d ; C) be the set of infinitely differentiable, compactly supported functions. Lemma 4.1 (Faris-Lavine theorem with periodic vector potentials). Let V and W be real-valued measurable functions on

R d , W ∈ L 2 loc (R d ; R) and A ∈ L 4 per (R d ; R d ) such that ∇ • A = 0 in the sense of distributions. Suppose that 1. there exists c, f ∈ R + such that W (x) ≥ -c|x| 2 -f , for a.a. x ∈ R d ; 2. 1 2 (-i∇ + A) 2 + V + W + 2c|x| 2 is essentially self-adjoint on C; 3. for some a < 1, a 2 (-i∇ + A) 2 + V is bounded below on C. Then 1 2 (-i∇ + A) 2 + V + W is essentially self-adjoint on C. The proof of the above lemma is postponed until Appendix A.1. We apply Lemma 4.1 with V ∈ L 2 per (R d ; R), W = εx β . The operator 1 2 (-i∇ + A) 2 + V + εx β + 2|x| 2 is essentially self-adjoint on the core C in view of [24, Theorem 3] (note that εx • e β ≥ -|x| 2 - ε 2 |e β | 2 4
). Moreover, since V is L 2 per (R d ; R), there exists 0 < a < 1, such that a 2 (-i∇+A) 2 +V is bounded below. This can be seen directly, or as a consequence of [START_REF] Leinfelder | Schrödinger operators with singular magnetic vector potentials[END_REF]Theorem 3]. Then, Lemma 4.1 gives that H ε β is essentially self-adjoint on C and therefore admits a unique self-adjoint extension on L 2 (R d ; C). Hence, the propagator of the associated Schrödinger equation is welldefined, and explicitly given by (e -itH ε β ) t∈R .

Proof of the second assertion. The self-adjointness of the time-dependent Hamiltonian

H ε β (t) = 1 2 (-i∇ + A -εe β t) 2 + V,
is a consequence of Lemma 4.1, by replacing A with (A -εe β t). To show the well-posedness of the dynamics, since H ε β (t) is R-periodic, it suffices to study its fibers. Hence we consider the dynamics of a Schrödinger equation with Hamiltonian

H ε β,k (t) = 1 2 (-i∇ + A + k -εe β t) 2 + V
and we use the following lemma on the dynamics generated by time-dependent perturbations of the free-particle Hamiltonian on L 2 per . Lemma 4.2. Let H 0 := -1 2 ∆ be the free-particle Hamiltonian on L 2 per , and a map

[0, T ] t → H 1 (t)
taking its values in the set of H 0 -bounded self-adjoint operators on L 2 per with relative bound lower than 1, that is: there exist 0 < a < 1 and b > 0 such that

∀t ∈ [0, T ], ∀φ ∈ H 2 per , H 1 (t)φ L 2 per ≤ a H 0 φ L 2 per + b φ L 2 per . (41)
Then, for all t ∈ [0, T ], the operator defined by

H(t) = H 0 +H 1 (t)
φ(t) = H(t)φ(t), φ(0) = φ 0 .
The proof of the above lemma is postponed to Appendix A.2.

For t ∈ [0, T ], k ∈ R d , we have H ε β,k (t) = H 0 + H 1 (t), with H 1 (t) = 1 2 (-i∇) • (A + k -εe β t) + (A + k -εe β t) • (-i∇) + (A + k -εe β t) 2 + V.
Using the Sobolev embeddings H 2 per ⊂ L ∞ per , H 1 per ⊂ L 6 per (recall that we assume d ≤ 3), the Coulomb gauge choice ∇ • A = 0 and the fact that A ∈ L 4 per (R d ; R d ) and V ∈ L 2 per (R d ; R), it is standard that H 1 satisfies the conditions of Lemma 4.2, and the result follows.

Proof of the third assertion. We first compute the fibers of the R-periodic operator γ ε β (t). Using [START_REF] Haldane | Model for a quantum Hall effect without landau levels: Condensedmatter realization of the" parity anomaly[END_REF], we have

γ ε β,k (t) = G * εte β U ε β (t)γ(0) U ε β (t) * G εte β k (42) = U ε β (t)γ(0) U ε β (t) * k+εte β = U ε β,k+εte β (t)γ k+εte β (0) U ε β,k+εte β (t) * (43)
=

N k+εte β n=1 | U ε β,k+εte β (t)u n,k+εte β U ε β,k+εte β (t)u n,k+εte β |.
Since the u n,k are in H 2 per , we deduce that

(J α γ ε β (t)) k = -(-i∇ + k + A) • e α γ ε β,k (t) = -∂ α H k γ ε β,k (t)
is trace-class (and finite-rank) uniformly in k ∈ B and therefore that the current

j ε α,β (t) = Tr J α γ ε β (t) is well-defined. As the function k → Tr ∂ α H k γ ε β,k (t) is R * -periodic, we also have j ε α,β (t) = -(2π) -d ˆB Tr ∂ α H k γ ε β,k (t) dk = -(2π) -d ˆB Tr ∂ α H k-εe β t γ ε β,k-εe β t (t) dk.

Perturbation theory for time-dependent Hamiltonians

In this section we consider the dynamics generated by a Hamiltonian H(s) = H(εt), and in particular its action on eigenspaces of H(0). We begin with some elementary properties of the Liouvillian in Section 5.1, then use it to study subspace perturbation theory in Section 5.2. We establish an adiabatic theorem in Section 5.3, and use it to study the time-dependent Hamiltonian H k-εe β t in Section 5.4. Finally, we prove a result in linear response with a remainder independent of the gap in Section 5.5.

The Liouvillian and its partial inverse

In order to formulate and interpret our results, it is convenient to make use of the formalism of the Liouvillian and its partial inverse, a classical tool in adiabatic theory and eigenvalue perturbation theory [START_REF] Teufel | Adiabatic Perturbation Theory in Quantum Dynamics[END_REF][START_REF] Kato | Perturbation Theory for Linear Operators[END_REF], although sometimes used implicitly. This formalism was for instance used in the context of transport properties in [START_REF] Bachmann | The adiabatic theorem and linear response theory for extended quantum systems[END_REF][START_REF] Teufel | Non-equilibrium almost-stationary states and linear response for gapped quantum systems[END_REF][START_REF] Monaco | Adiabatic currents for interacting electrons on a lattice[END_REF]. Recall that if h is a bounded self-adjoint operator on a Hilbert space H f , the Liouvillian L h associated with h is the bounded linear operator on L(H f ) (such a mathematical object is sometimes called a superoperator in the physics literature) defined by

∀A ∈ L(H f ), L h A = [h, A]. (44) Note that if A is self-adjoint, then L h A is anti-self-adjoint (iL h A = i[h, A] is self-adjoint). The restriction of L h to the space S 2 (H f ) of Hilbert-Schmidt operators on H f is self-adjoint: for all A, B ∈ S 2 (H f ), (L h A, B) S2 = Tr ([h, A] * B) = Tr ((A * h -hA * )B) = Tr (A * (hB -Bh)) = (A, L h B) S2 . (45)
The operator L h is to density matrices what the Hamiltonian h is to pure states: it is the infinitesimal generator of the norm-continuous unitary group (U h (t)) t∈R on L(H f ) defined by

∀A ∈ L(H f ), U h (t)A = e -ith Ae ith . ( 46
)
In the case when h is an unbounded self-adjoint operator, (44) does not make sense for all A ∈ L(H f ), but it is still possible to define the Liouvillian L h as the infinitesimal generator of the strongly-continuous unitary group (U h (t)) t∈R on L(H f ) defined by (46). It is then an unbounded operator on L(H f ), self-adjoint on S 2 (H f ).

If H f is of finite-dimension N f , the action of L h is easily understood in an orthonormal eigenbasis (e n ) 1≤n≤N f of h with associated eigenvalues

λ 1 ≤ • • • ≤ λ N f . Then, L h |e n e m | = (λ n -λ m )|e n e m |.
The operator L h is not invertible (for instance, L h |e n e n | = 0). However, it is invertible when restricted to the subspace of block off-diagonal matrices, i.e. matrices A such that A nn = A mm = 0 for n, n ≤ N < m, m for a given N such that λ N +1 > λ N . Its partial inverse L + h,N is given by

L + h,N 1≤n≤N, N <m≤N f A mn |e m e n | + A nm |e n e m | = 1≤n≤N, N <m≤N f A mn |e m e n | -A nm |e n e m | λ m -λ n (47)
and L + h,N is bounded in operator norm by 1 λ N +1 -λ N . More generally, if h is an unbounded self-adjoint operator, let I be a closed bounded interval of R, and assume that

g := min (1, dist (I, σ(h) \ (σ(h) ∩ I))) > 0.
The associated spectral projector is

P I,h := 1 I (h) = 1 2πi ˛C(z -h) -1 dz, (48)
where C is a Cauchy contour in the complex plane such that σ(h)∩I is inside C and σ(h)\(σ(h)∩I) is outside C. Generalizing the terminology of the finite-dimensional case, we call off-diagonal operators (with respect to the splitting of H f induced by P I,h ) the elements of the closed subspace

L OD h,I := {A ∈ L(H f ) | P h,I AP h,I = (1 -P h,I )A(1 -P h,I ) = 0}
of L(H f ). This defines a S 2 -orthogonal splitting of operators into their diagonal and off-diagonal parts. It is easily seen that L OD h,I is L h -stable, and that L h is invertible on L OD h,I with a bounded inverse. We denote its partial inverse by L + h,I , extended to all of L(H f ) by imposing that it vanishes on diagonal operators. We then have

L h,I L + h,I A = L + h,I L h,I A = P h,I A(1 -P h,I ) + (1 -P h,I )AP h,I for all A ∈ L(H f ).
It is easy to check that L + h,I has an explicit contour integral representation:

L + h,I A = 1 2πi ˛C(z -h) -1 [P h,I , A](z -h) -1 dz, ∀A ∈ L(H f ), ( 49
)
where C is a contour as above. From (49), we see that, when P h,I is of finite rank N , L + h,I A is of rank of most 2N .

Subspace perturbation theory

The Liouvillian is a powerful tool to write concisely the results of subspace perturbation theory, which studies the time dependence of a gapped subspace of a time-dependent Hamiltonian. We consider T > 0 and (H(s)) s∈[0,T ) a family of self-adjoint operators on a Hilbert space H f sharing the same domain D ⊂ H f and satisfying the following assumptions: H1 H(s) ≥ 1 for all s ∈ [0, T ); H2 for each φ ∈ D, the map s → H(s)φ is in C n ([0, T ), H f ) for some n ≥ 1. For all 1 ≤ l ≤ n, the operator H (l) (s) is self-adjoint on H f for all s ∈ [0, T ), and

α l := sup s∈[0,T ) H (l) (s)H(s) -1 L(H f ) < ∞; (50
H(s)P (s) ≤ M, (52)
H(s)P (l) (s) ≤ C l g(s) l+1 , (53)
H(s)(L + ) (l) (s)A ≤ C l g(s) l+3 A . ( 54
)
In addition, P (l) (s) has rank at most (l + 1)N , and (L + ) (l) (s)A has rank at most c l N where c l is a constant that only depends on l (in particular, c 0 = 2 and c 1 = 10).

Remark 5.2. The powers of the gap in the bounds ( 53) and ( 54) are too pessimistic, as could be shown by a more detailed analysis. For instance, in the case l = 0, L + (s) can be seen from the arguments at the beginning of this section to be bounded by a constant times 1 g(s) . Similarly, the operator Ṗ is bounded by a constant times 1 g(s) , using (51). Nevertheless, the above bounds are more straightforward to establish and will suffice for our purposes.

Proof. Differentiating L H(s) P (s) = 0, we get L H(s) Ṗ (s) = [P (s), Ḣ(s)].

Since both [P (s), Ḣ(s)] and Ṗ (s) are off-diagonal operators (the first by direct calculation, the second by differentiating the relationship P (s) 2 = P (s)), we deduce (51). By the functional calculus, H(s)P (s) = H(s)1 I(s) (H(s)) ≤ a + (s) ≤ M , whence (52).

In the following we take for C(s) the rectangular contour centered at the center of I(s), of length |I(s)| + g(s) and height g(s), so that

|C(s)| ≤ 2M + 4 and for all z ∈ C(s), 1 z -H(s) ≤ 2 g(s) . ( 55
)
We use the integral representation (48):

P (s) = 1 2πi ˛C(s) 1 z -H(s) dz. ( 56
)
Using for all z ∈ C(s) the bound

H(s) z -H(s) = sup λ∈σ(H(s)) λ z -λ ≤ 1 + sup λ∈σ(H(s)) z z -λ ≤ 1 + 2(M + g(s)) g(s) ≤ 2M + 3 g(s) (57)
establishes (53) for l = 0.

The contour C(s) in (56) above can be kept fixed equal to C(s 0 ) for s in a neighborhood of any s 0 ∈ [0, T ). Using

d ds 1 z -H(s) = 1 z -H(s) Ḣ(s) 1 z -H(s) (58) it follows that P ∈ C 1 ([0, T ), L(H f )) and Ṗ (s) = 1 2πi ˛C(s) 1 z -H(s) Ḣ(s) 1 z -H(s) dz.
Using the bounds (50), ( 55) and (57), it follows that

H(s) Ṗ (s) ≤ (2M + 3)(2M + 4)α 1 πg(s) 2
which proves (53) for l = 1. The general case for l > 1 follows from repeated application of the chain rule to (56) and ( 58), and the bounds (50), ( 55) and (57).

The differentiability and bounds on the inverse Liouvillian are treated using the same arguments on the representation

L + (s)A = 1 2πi ˛C(s) 1 z -H(s) [P (s), A] 1 z -H(s) dz.
Let (u 0 n ) n=1,...,N be an orthonormal basis of P (0). Then the solutions to the parallel transport equation un (s) = Ṗ (s)u n (s) with u n (0) = u 0 n are easily checked to be a C n orthogonal basis of RanP (s). It follows that one has

P (l) (s) = N n=1 l m=0 l m |u (m) n (s) u (l-m) n (s)|.
Therefore, P (l) (s) is of rank at most (l + 1)N . From the integral representation of L + (s) (see ( 49)), it follows that, for any bounded operator A, L + (s)A is of rank at most 2N . Its derivatives are sums of terms which all contain as a factor P (s) or one of its derivative, and the result follows with c l = 2 k1+k2+k3=l, kj ∈N (k 2 + 1).

Adiabatic theory

The following proposition is an adaptation in our context of the classical adiabatic theorem that the Schrödinger evolution with a slowly evolving Hamiltonian H(εt) approximately preserves gapped eigenspaces [START_REF] Teufel | Adiabatic Perturbation Theory in Quantum Dynamics[END_REF]. We explicitly compute the corrections to first order in ε.

Proposition 5.3. Assume the same hypotheses as in Proposition 5.1. Let (U ε (t, t )) 0≤t ≤t<ε -1 T be the propagator associated with the family of time-scaled Hamiltonians (H(εt)

) t∈[0,ε -1 T ) , i.e. i ∂U ε ∂t (t, t ) = H(εt)U ε (t, t ), t ∈ [t , ε -1 T ) U ε (t , t ) = Id, ( 59
)
and U ε (t) = U ε (t, 0). For all ε ≥ 0 and t ∈ [0, ε -1 T ), it holds

U ε (t)P (0)U ε (t) * = P (εt) + iεL + (εt) Ṗ (εt) -iεU ε (t) L + (0) Ṗ (0) U ε (t) * + R ε (t), ( 60
) with R ε (t) = -iε 2 ˆt 0 U ε (t, t ) d ds L(s) -1 Ṗ (s) s=εt U ε (t, t ) * dt . ( 61
)
In addition, we have the following estimates:

∀0 ≤ t ≤ t < ε -1 T, H(εt)U ε (t, t )H(εt ) -1 L(H f ) ≤ e α1ε(t-t ) , (62)
H(εt) 1/2 U ε (t, t )H(εt ) -1/2 L(H f ) ≤ e α1ε(t-t )/2 . ( 63
)
Proof. The existence and uniqueness of the strongly-continuous unitary propagator (U ε (t, t )) satisfying (59) can be obtained using (50) for l = 1, and Theorem X.70 and the arguments in the proof of Theorem X.71 in [START_REF] Michael | Methods of Modern Mathematical Physics. II. Fourier analysis, self-adjointness[END_REF]. We pass to the interaction picture defined by H(εt) and compute the evolution of a C 1 time-dependent Hilbert-Schmidt observable A ε (t) in that picture:

d dt (U ε (t) * A ε (t)U ε (t)) = U ε (t) * Ȧε (t) + i[H(εt), A ε (t)] U ε (t). ( 64
)
We first apply (64) to A ε (t) = P (εt) and obtain

d dt (U ε (t) * P (εt)U ε (t)) = ε U ε (t) * Ṗ (εt)U ε (t). ( 65
)
Estimating this to be of size ε is not enough because we look at long time scales. What allows us to proceed further is that this quantity is oscillating on a timescale of order O(1). Indeed, applying (64) to A ε (t) = L + (εt) Ṗ (εt), for which [H(εt), A ε (t)] = Ṗ (εt), we obtain

U ε (t) * Ṗ (εt)U ε (t) = -i d dt U ε (t) * (L + (εt) Ṗ (εt))U ε (t) + iU ε (t) * d dt L + (εt) Ṗ (εt) U ε (t).
Integrating (65) over [0, t] and using the above equality leads to

U ε (t) * P (εt)U ε (t) = P (0) + ε ˆt 0 U ε (t ) * Ṗ (εt )U ε (t ) dt = P (0) -iεU ε (t) * L + (εt) Ṗ (εt) U ε (t) + iεL + (0) Ṗ (0) + r ε (t) with r ε (t) = iε ˆt 0 U ε (t ) * d dt L + (εt ) Ṗ (εt ) U ε (t ) dt = iε 2 ˆt 0 U ε (t ) * d ds L + (s) Ṗ (s) s=εt U ε (t ) dt .
This implies

U ε (t)P (0)U ε (t) * = P (εt) + iεL + (εt) Ṗ (εt) -iεU ε (t) L + (0) Ṗ (0) U ε (t) * + R ε (t), ( 66
) with R ε (t) = -iε 2 ˆt 0 U ε (t, t ) d ds L + (s) Ṗ (s) s=εt U ε (t, t ) * dt ,
which establishes (60).

Let us now prove (62). Let ψ ∈ D. For all t ∈ [t , ε -1 T ), we set ψ ε (t) = U ε (t, t )ψ and φ ε (t) = H(εt)ψ ε (t). We have

i dφ ε dt (t) = i d dt (H(εt)ψ ε (t)) = H(εt)φ ε (t) + iε Ḣ(εt)H(εt) -1 φ ε (t),
from which we obtain

φ ε (t) = U ε (t, t )H(εt )ψ + iε ˆt t Ḣ(εt )H(εt ) -1 φ ε (t ) dt ,
and finally

φ ε (t) H f ≤ H(εt )ψ H f + α 1 ε ˆt t φ ε (t ) H f dt .
By the Grönwall lemma,

H(εt)U ε (t, t )ψ H f = φ ε (t) H f ≤ H(εt )ψ H f e α1ε(t-t ) .
Applying this inequality to ψ = H(εt ) -1 φ for all φ ∈ H f gives (62). We obtain (63) by interpolation (see e.g. [32, Section IX.4, Proposition 9]).

The third term

-U ε (t) iL + (0) Ṗ (0) U ε (t) *
of the right-hand side of (60) is oscillatory, and can be written as the derivative of a bounded function up to higher order terms. Its time-average therefore becomes negligible in the considered regimes. Let us introduce the space

L OD (s) := {A ∈ L(H f ) | P (s)AP (s) = (1 -P (s))A(1 -P (s)) = 0}
of bounded off-diagonal operators relatively to the decomposition H f = Ran(P (s)) ⊕ Ker(P (s)).

Lemma 5.4. Under the assumptions of Propositions 5.1 and 5.3, we have for any self-adjoint operator A ∈ L OD (0),

U ε (t)AU ε (t) * = d dt iL + (εt) (U ε (t)AU ε (t) * ) + R ε A (t), ( 67
)
where

R ε A (t) = 2U ε (t)r ε (t)Ar ε (t)U ε (t) * -U ε (t)(1-2P (0))Ar ε (t)U ε (t) * +h.c. +εi dL + ds (εt) U ε (t)AU ε (t) * and r ε (t) = -iεU ε (t) * L + (εt) Ṗ (εt) U ε (t) + iεL + (0) Ṗ (0) + U ε (t) * R ε (t)U ε (t).
Proof. We have

d dt iL + (εt) (U ε (t)AU ε (t) * ) = εi dL + ds (εt) (U ε (t)AU ε (t) * ) + L + (εt)L(εt) (U ε (t)AU ε (t) * ) = εi dL + ds (εt)U ε (t)AU ε (t) * + P (εt)U ε (t)AU ε (t) * (1 -P (εt)) + h.c.,
and we deduce from (60) that P (εt)U ε (t) = U ε (t) (P (0) + r ε (t)). We therefore have

P (εt)U ε (t)AU ε (t) * (1 -P (εt)) + h.c. = U ε (t) (P (0) + r ε (t)) A(1 -P (0) -r ε (t))U ε (t) * + h.c. = U ε (t)AU ε (t) * + (U ε (t) (1 -2P (0)) Ar ε (t)U ε (t) * + h.c.) + 2U ε (t)r ε (t)Ar ε (t)U ε (t) * ,
where we have used that A = P (0)A(1 -P (0)) + (1 -P (0))AP (0).

Application to coherent transport in Bloch representation

Let H be the periodic magnetic Hamiltonian defined in (1), J the current operator whose components are defined in [START_REF] Cancès | Numerical quadrature in the Brillouin zone for periodic Schrödinger operators[END_REF], µ F the Fermi level,

µ := 1 + min σ(H) and η = max |e|≤|eα|,|e β | (J • e)(H + µ) -1 L(L 2 (R d ;C)) < ∞. Let k ∈ R d . Assume that λ N k +1,k -λ N k ,k > 0 and set s k = inf{s > 0 | g k (s) = 0} where g k (s) := min(1, λ N k +1,k-se β -λ N k ,k-se β ).
We consider the family of Hamiltonians

H(s) := H k-e β s + µ. (68)
We have

Ḣ(s) = -∂ β H k-se β = -((-i∇ + k + A -se β ) • e β ) = J β,k-se β , (69)
Ḧ(s) = |e β | 2 Id L 2 per , (70)
and so hypotheses H1-H3 of Proposition 5.1 are satisfied with

H f = L 2 per , D = H 2 per , T = s k , n arbitrarily large, α 1 ≤ η, α 2 = |e β | 2 , α l = 0 for l ≥ 3, a -(s) = min σ(H)+µ, a + (s) = λ N k ,k-se β +µ, M = max k ∈B λ N k +1,k + µ, g(s) = g k (s), and N = N k . Corollary 5.5. Let k ∈ R d such that λ N k +1,k -λ N k ,k > 0. Then, for all ε > 0 and t ∈ [0, ε -1 s k ), the operator ∂ α H k-εe β t γ ε β,k-εe β t (t) is in S 1,
per , and we have

Tr(∂ α H k-εe β t U ε β,k (t)γ k (0) U ε β,k (t) * ) =∂ α Tr H k-εe β t P N k ,k-εe β t +iε Tr ∂ α H k-εe β t L + N k ,k-εe β t ∂ β P N k ,k-εe β t -iε Tr ∂ α H k-εe β t U ε β,k (t) L + N k ,k ∂ β P N k ,k U ε β,k (t) * + R ε k (t), (71)
where each term of the right-hand side is a well-defined real number and L + N,k is a shorthand notation for the inverse Liouvillian L

+ H k ,[λ 1,k ,λ N,k ] .
In addition, we have the following bounds

|Tr(∂ α H k-εe β t U ε β,k (t)γ k (0) U ε β,k (t) * )| ≤ Ce ηεt , (72)
|∂ α Tr H k-εe β t P N k ,k-εe β t | ≤ C, |ε Tr ∂ α H k-εe β t L + N k ,k-εe β t ∂ β P N k ,k-εe β t | ≤ C ε g k (εt) 4 , |ε Tr ∂ α H k-εe β t U ε β,k (t) L + N k ,k ∂ β P N k ,k U ε β,k (t) * | ≤ C εe ηεt g k (0) 4 , (73)
|R ε k (t)| ≤ Cε 2 te ηεt min s∈[0,εt] g k (s) 6 , (74)
for a constant C ∈ R + independent of k, ε and t.

Proof. Applying the second assertion in Proposition 5.3, we get

U ε β,k (t)γ k (0) U ε β,k (t) * =P N k ,k-εe β t + iεL + N k ,k-εe β t ∂ β P N k ,k-εe β t -iε U ε β,k (t) L + N k ,k ∂ β P N k ,k U ε β,k (t) * + R ε k (t). (75)
Each term A in (75) being a finite-rank self-adjoint operator, it holds

∂ α H k-εe β t A S1 ≤ Rank(A) ∂ α H k-εe β t A ≤ ηRank(A) (H k-εe β t + µ)A ,
and again by Proposition 5.3 we get

Tr(∂ α H k-εe β t U ε β,k (t)γ k (0) U ε β,k (t) * ) = Tr ∂ α H k-εe β t P N k ,k-εe β t +iε Tr ∂ α H k-εe β t L + N k ,k-εe β t ∂ β P N k ,k-εe β t -iε Tr ∂ α H k-εe β t U ε β,k (t) L + N k ,k ∂ β P N k ,k U ε β,k (t) * + Tr ∂ α H k-εe β t R ε k (t) with R ε k (t) = iε 2 ˆt 0 U ε β,k (t, t )∂ k β L + N k ,k-εe β t ∂ k β P N k ,k-εe β t U ε β,k (t, t ) * dt . (76)
It results from the Hellmann-Feynman formula that

Tr ∂ α H k-εe β t P N,k-εe β t = ∂ α Tr H k-εe β t P N,k-εe β t .
Finally, using Propositions 5.1 and 5.3, we obtain the bounds (72)-(74). In particular, 6 , where C ∈ R + is independent of k, ε and t.

|R ε k (t)| = Tr ∂ α H k-εe β t R ε k (t) ≤ 4N k ε 2 ηt sup t ∈[0,t) e ηε(t-t ) H k-εe β t ∂ k β L + N k ,k-εe β t ∂ k β P N k ,k-εe β t ≤ C ηε 2 te ηεt inf s∈[0,εt) g k (s)
Remark 5.6. The decomposition (71) will be key to computing the current in insulators, nondegenerate metals and semimetals. The first three terms in the right-hand side of (71) have different physical meanings. The first term is the adiabatic term: electrons simply are transported adiabatically across the Brillouin zone. This term will be responsible for the ballistic transport of electrons in metals. The second is the first-order static response, and will be the cause of the Hall conductivity in insulators. The third is oscillatory, and is related to the AC response of solids (not treated here). This decomposition only makes sense for a non-zero gap; in particular, it cannot be used to compute the contribution to the current for k points close to Dirac points for semimetals.

Linear response

We now aim at obtaining an expansion of the current to first order in ε for a given t, based on a Dyson expansion instead of the adiabatic theorem. This is a classical computation in response theory, sometimes known as the Kubo formula [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. I. general theory and simple applications to magnetic and conduction problems[END_REF]. In contrast to the previous result, this gives a remainder that does not depend on a gap, and will therefore be useful for the study of semimetals near Dirac points. Proposition 5.7. Let H be the periodic magnetic Hamiltonian defined in [START_REF] Bachmann | The adiabatic theorem and linear response theory for extended quantum systems[END_REF]. Under the additional assumptions that V ∈ H 1 per and A ∈ (H 2 per ) d , there exists a constant

C ∈ R + such that for all k ∈ R d such that λ N k +1,k -λ N k ,k > 0, we have for all ε, t ∈ R + , Tr ∂ α H k-εe β t U ε β,k (t) γ k (0) U ε β,k (t) * =∂ α Tr (H k γ k (0)) -εt∂ α ∂ β (Tr(H k γ k (0))) + iε Tr ∂ α H k (e -itL k -1)L + k ∂ β γ k (0) + ρ ε k (t), (77)
with, when εt ≤ 1,

|ρ ε k (t)| ≤ Cε 2 t 3 (1 + t 3 ). ( 78
) Proof. Let k be such that λ N k +1,k -λ N k ,k > 0. Since k → Tr (H k γ k (0)
) is real-analytic in a neighborhood of k, we have by Hellmann-Feynman theorem

∂ α Tr (H k γ k (0)) = Tr (∂ α H k γ k (0)) and ∂ α ∂ β Tr (H k γ k (0)) = Tr (∂ α ∂ β H k γ k (0)) + Tr (∂ α H k ∂ β γ k (0)) .
We also have

∂ α H k-εe β t = ∂ α H k -εte α • e β . It follows that ρ ε k (t) =Tr ∂ α H k U ε β,k (t) γ k (0) U ε β,k (t) * -Tr (∂ α H k γ k (0)) + εtTr (∂ α H k ∂ β γ k (0)) -iε Tr ∂ α H k (e -itL k -1)L + k ∂ β γ k (0) .
We now expand the first term in the right-hand side of this equation. We set µ := 1 + min σ(H),

H(s) := H k-se β + µ, A = ∂ α H k , I 0 = [0, 1 2 (λ N k ,k + λ N k +1,k ) + µ], P (s) = 1 I0 (H(s)).
It holds

H(s) = h 0 + sh 1 + s 2 |e β | 2 2 with h 0 = H k + µ and h 1 = J β,k = -∂ β H k .
The operators h 0 , h 1 and A are self-adjoint on L 2 per and we have h 0 ≥ 1 and h 1 h

-1/2 0
and Ah -1/2 0 bounded. Besides, P (s) = γ k-se β (0), so that Ṗ (0) = -∂ β γ k (0). Let (U ε (t, t )) t,t ∈R be the propagator associated with the family (H(εt)) t∈R and U ε (t) := U ε (t, 0). We have U ε (t, t ) = e -iµ(t-t ) U ε k (t, t ) and U ε (t) = e -iµt U ε k (t). With these notations, we have

ρ ε k (t) = Tr (AU ε (t) P (0)U ε (t) * ) -Tr (AP (0)) -εtTr A Ṗ (0) + iε Tr A(e -itL0 -1)L + 0 Ṗ (0) ,
where L 0 = L h0,I0 and L + 0 = L + h0,I0 , and we focus on expanding the operator U ε (t) P (0)U ε (t) * close to t = 0. Lemma 5.8. We have

U ε (t) P (0)U ε (t) * = P (0) + ε t Ṗ (0) -i e -itL0 -1 (L + 0 Ṗ (0)) + Π ε 2 (t) (79) with ρ ε k (t) = Tr (AΠ ε 2 (t))
. Moreover, we have the bound (78)

|ρ ε k (t)| ≤ Cε 2 t 3 (1 + t 3 ).
Lemma 5.8 closes the proof of Proposition 5.7.

Proof of Lemma 5.8. We deduce from the Dyson expansion that

U ε (t) = U 0 (t) + V ε (t) + W ε (t),
where U 0 (t) = e -ith0 and

V ε (t) = -iε ˆt 0 U 0 (t -t )t h 1 U 0 (t ) dt , W ε (t) = ε 2 -i t 3 6 U 0 (t) + ˆt 0 ˆt 0 U ε (t, t )t (h 1 + εt /2)U 0 (t -t )t (h 1 + εt /2)U 0 (t ) dt dt .
This induces U ε (t)

P (0)U ε (t) * = P (0) + Π ε 1 (t) + Π ε 2 (t) where Π ε 1 (t) = V ε (t)P (0)U 0 (t) * + h.c. = -iε ˆt 0 t U 0 (t -t )[h 1 , P (0)]U 0 (t -t ) * dt , Π ε 2 (t) = V ε (t)P (0)V ε (t) * + W ε (t)P (0)(U 0 (t) + V ε (t)) * + h.c. + W ε (t)P (0)W ε (t) * .
We first analyze Π ε 1 (t) by computing

U 0 (t -t )[h 1 , P (0)]U 0 (t -t ) * = -e -i(t-t )L0 L 0 Ṗ (0) = i d dt e -i(t-t )L0 Ṗ (0) = d 2 dt 2 e -i(t-t )L0 L + 0 Ṗ (0),
where we have used Ṗ (0) = L + 0 [P 0 , h 1 ] and Ṗ (0) = L + 0 L 0 Ṗ (0). Using integration by parts, we obtain

Π ε 1 (t) = ε t Ṗ (0) -i e -itL0 -1 (L + 0 Ṗ (0))
and (79) follows.

We now work on the bound (78). For that purpose, we introduce the following quantities, which are independent of k, ε and t:

ν 0 = max |e|≤|eα|,|e β | (J • e)(H + µ) -1/2 L(L 2 (R d ;C)) , ν 1 = max |e|≤|eα|,|e β | (H + µ) 1/2 (J • e)(H + µ) -1 L(L 2 (R d ;C)) , ν 2 = max |e|≤|eα|,|e β | (H + µ)(J • e)(H + µ) -2 L(L 2 (R d ;C)) , λ = max k∈R d ,|k-k |≤|e β | (H k + µ) 1/2 (H k + µ) -1/2 L(L 2 per ) .
Note that the assumptions A ∈ (L 4 per ) d , ∇ • A = 0, and V ∈ L 2 per are sufficient to ensure that the quantities ν 0 , ν 1 and λ are finite. Besides, since

h 0 h 1 h -2 0 ≤ (H + µ)J β (H + µ) -2 and (H + µ)J β (H + µ) -2 = J β (H + µ) -1 -2i d α=1 (∂ α A β -∂ β A α )J α (H + µ) -2 -(∆A β )(H + µ) -2 + i∂ β V (H + µ) -2 , we deduce from the assumptions A ∈ (H 2 per ) d and V ∈ H 1 per that (H + µ)J β (H + µ) -2 < ∞, hence that ν 2 < ∞.
We now aim at controlling ρ k (t) thanks to ν 0 , ν 1 , ν 2 and λ. Using the relations P (0) = P (0) 2 and P (0) = h -m 0 h m 0 P (0) with

h m 0 P (0) ≤ (µ F + µ) m , h -1 0 ≤ h -1/2 0 ≤ 1 and Rank(P (0)) = N k ≤ N := max k N k , we deduce |ρ ε k (t)| ≤ N (µ F + µ) 2 AV ε (t)h -1 0 V ε (t)h -1 0 + 2(µ F + µ) 3 AW ε (t)h -2 0 ( h -1 0 + V ε (t)h -1 0 ) + (µ F + µ) 3 AW ε (t)h -2 0 W ε (t)h -1 0 ≤ N ν 0 (µ F + µ) 2 h 1/2 0 V ε (t)h -1 0 2 + (µ F + µ) 3 h 1/2 0 W ε (t)h -2 0 (2 + 2 V ε (t)h -1 0 + W ε (t)h -1 0
) .

Next, we get

V ε (t)h -1 0 ≤ εt 2 2 h 1 h -1 0 ≤ εt 2 2 ν 0 , h 1/2 0 V ε (t)h -1 0 ≤ εt 2 2 ν 1 , W ε (t)h -1 0 ≤ ε 2 t 3 1 6 + ν 0 ν 1 t 8 + (εt)t ν 0 30 + ν 1 20 + (εt) 2 t 1 72 , h 1/2 0 W ε (t)h -2 0 ≤ ε 2 t

Insulators

In this section and the following ones, we use the notation O(f (ε, t, t , δ)) to denote a term that is bounded in absolute value by Cf (ε, t, t , δ), where C is a constant that can depend on the system under consideration (through A, V , µ F , e α and e β ), but not on the parameters ε, t, t , δ. We will use the notation γ 0 k for γ k (0). We now prove Theorem 2.7. For insulators, N k = N ins for all k, and λ N k +1,k -λ N k ,k , hence g k , is uniformly bounded away from zero and s k = +∞. We use the notation L + k for L + Nins,k . We apply Corollary 5.5 and obtain by integrating over the Brillouin zone

j ε α,β (t) =-(2π) -d ˆB ∂ α Tr H k-εe β t γ 0 k-εe β t dk-iε(2π) -d ˆB Tr ∂ α H k-εe β t L + k-e β t ∂ β γ 0 k-εe β t dk + iε(2π) -d ˆB Tr ∂ α H k-εe β t U ε β,k (t) (L + k ∂ β γ 0 k ) U ε β,k (t) * dk + O ε 2 te ηεt .
As mentioned in Remark 5.6, these three terms are adiabatic, static and oscillatory respectively.

• The first term of the right-hand side vanishes for all t, as the integral of the derivative of the smooth periodic function k → Tr H k γ 0 k on a unit cell. • The second term is dealt with using the relation

L + k ((∂ α H k ) OD ) = [γ 0 k , ∂ α γ 0 k ],
where

(∂ α H k ) OD = γ 0 k (∂ α H k )(1 -γ 0 k ) + (1 -γ 0 k )(∂ α H k )γ 0 k . By periodicity, we have ˆB Tr ∂ α H k-εe β t L + k-εe β t ∂ β γ 0 k-εe β t dk = ˆB Tr ∂ α H k L + k ∂ β γ 0 k dk,
and we observe that

Tr ∂ α H k L + k ∂ β γ 0 k = Tr (∂ α H k ) OD L + k ∂ β γ 0 k = Tr L + k ((∂ α H k ) OD )∂ β γ 0 k = Tr [γ 0 k , ∂ α γ 0 k ]∂ β γ 0 k , (80)
so that ˆB Tr ∂ α H k-εe β t L + k-εe β t ∂ β γ 0 k-εe β t dk = ˆB Tr γ 0 k [∂ α γ 0 k , ∂ β γ 0 k ] dk.
• We now focus on the time-average of the oscillating term

ω ε (t) := 1 t ˆt 0 dt ˆB Tr ∂ α H k-εe β t U ε β,k (t ) (iL + k ∂ β γ 0 k ) U ε β,k (t ) * dk.
In order to bound this term, we apply Lemma 5.4 to A = iL + k ∂ β γ 0 k , which is a self-adjoint off-diagonal operator for the decomposition L 2 per = Ran(γ 0 k ) ⊕ Ker(γ 0 k ). We thus get

U ε β,k (t) iL + k ∂ β γ 0 k U ε β,k (t) * = d dt iL + k-εe β t U ε β,k (t)(iL + k ∂ β γ 0 k ) U ε β,k (t) * + R ε k (t),
where

R ε k (t) = 2 U ε β,k (t)r ε k (t)(iL + k ∂ β γ 0 k )r ε k (t) U ε β,k (t) * -U ε β,k (t)(1 -2γ 0 k )(iL + k ∂ β γ 0 k )r ε k (t) U ε β,k (t) * + h.c. -εi∂ β L + k-εte β U ε β,k (t)(iL + k ∂ β γ 0 k ) U ε β,k (t)) * and r ε k (t) = iε U ε β,k (t) * L + k-εte β ∂ β γ 0 k-εte β U ε β,k (t) -iεL + k ∂ β γ 0 k + U ε β,k (t) * R ε k (t) U ε β,k (t),
where R ε k (t) is defined in (76). Therefore,

Tr ∂ α H k-εe β t U ε β,k (t) (iL + k ∂ β γ 0 k ) U ε β,k (t) * = Tr ∂ α H k-εe β t d dt U ε β,k (t) (iL + k ∂ β γ 0 k ) U ε β,k (t) * + Tr ∂ α H k-εe β t R ε k (t) = d dt Tr ∂ α H k-εe β t U ε β,k (t) (iL + k ∂ β γ 0 k ) U ε β,k (t) * + Tr ∂ α H k-εe β t R ε k (t) since Tr d dt (∂ α H k-εe β t ) U ε β,k (t) (iL + k ∂ β γ 0 k ) U ε β,k (t) * = -ε Tr ∂ kαk β H k-εe β t U ε β,k (t) (iL + k ∂ β γ 0 k ) U ε β,k (t) * = -εe α • e β Tr U ε β,k (t) (iL + k ∂ β γ 0 k ) U ε β,k (t) * = -εe α • e β Tr iL + k ∂ β γ 0 k = 0,
where we have used the fact that ∂ kαk β H k = -e α • e β and the off-diagonal character of iL + k ∂ β γ 0 k . Hence, using the bounds from Proposition 5.1, we obtain

ω ε (t) = 1 t ˆB Tr ∂ α H k-εe β t U ε β,k (t) (iL + k ∂ β γ 0 k ) U ε β,k (t) * dk + 1 t ˆt 0 dt ˆB Tr ∂ α H k-εe β t R ε k (t ) dk = O 1 t + ε e ηεt .
The result follows.

Metals

We prove the two assertions of Theorem 2.8 in sequence.

Linear response

We prove the first assertion of Theorem 2.8: We first note that, for ε > 0 small enough and t ≤ 1 ε ε θ , the function k → λ N k +1,k-εe β t -λ N k ,k-εe β t is bounded away from zero, and therefore so is g k (εt). We can therefore apply Corollary 5.5 on each B N to obtain

j ε α,β (t) = (2π) -d N ∈N - ˆBN Tr ∂ α H k-εe β t P N,k-εe β t dk (81)
-iε

ˆBN Tr ∂ α H k-εe β t L + k-εe β t ∂ β P N,k-εe β t dk +iε ˆBN Tr ∂ α H k-εe β t U ε β,k (t) L + k ∂ β P N,k U ε β,k (t) * dk + O ε 2 te ηεt
and so

j ε α,β (t) = -(2π) -d N ∈N ˆBN Tr ∂ α H k-εe β t P N,k-εe β t dk + O(ε) (82)
when t ≤ 1 ε ε θ . In contrast to the case of insulators however, the adiabatic term

-(2π) -d N ∈N ˆBN ∂ α Tr H k-εe β t P N,k-εe β t dk = εt(2π) -d N ∈N ˆBN ∂ α ∂ β (Tr(H k P N,k )) dk + O(ε 2 t 2 )
has a non-zero first-order contribution (the zeroth-order term vanishes by Proposition 2.1). The proportionality factor is computed by the Stokes formula as

N ∈N * ˆBN ∂ α ∂ β Tr H k γ 0 k dk = N ∈N * ˆBN ∂ α ∂ β N n=1 λ n,k = N ∈N * ˆSN - ˆSN-1 ∂ α N n=1 λ n,k (ds • e β ) = N ∈N * ˆSN ∂ α λ N,k (ds • e β ) = (2π) d D α,β (83)
and the result follows.

Bloch oscillations

Under the assumptions of the second assertion, N k is either N met or N met -1, and in both cases

λ N k +1,k-εe β t -λ N k ,k-εe β t
is bounded away from zero uniformly in k,t. We can therefore apply Corollary 5.5 and obtain

j ε α,β (t) = -(2π) -d N ∈N ˆBN ∂ α Tr H k-εe β t P N k ,k-εe β t dk + O((ε + ε 2 t)e ηεt ).
From the decomposition

P N k ,k-εe β t = P Nmet-1,k-εe β t + 1(λ Nmet,k ≤ µ F )|u Nmet,k-εe β t u Nmet,k-εe β t |
and since k → P Nmet-1,k-εe β t is smooth and R * -periodic, we have

j ε α,β (t) = -(2π) -d ˆB1(λ Nmet,k ≤ µ F ) u Nmet,k-εe β t |∂ α H k-εe β t |u Nmet,k-εe β t dk + O((ε + ε 2 t)e ηεt ) = -(2π) -d ˆB1(λ Nmet,k ≤ µ F )∂ α λ Nmet,k-εe β t dk + O((ε + ε 2 t)e ηεt ),
which concludes the proof.

Semi-metals

We prove here Theorem 2.10. We decompose the integral defining j ε α,β (t) into several parts depending whether one integrates far from the Dirac points or not.

We introduce a small parameter δ > 0 controlling the size of the neighborhood of the Dirac points, which is independent of t, ε. We decompose B as the disjoint union

B = B δ out ∪ (∪ i∈I B δ i ) with B δ i = B(k i , δ), where δ > 0 is small enough so that B δ out ⊂ {k ∈ B, λ Nsm,k ≤ µ F -cδ}
for some constant c > 0. Note that this decomposition is time-reversal symmetric in the sense that

-B δ out = B δ out and -(∪ i∈I B δ i ) = (∪ i∈I B δ i ).
We work in the regime εt δ 1, ε δ 1.

In the following analysis, we first treat the regions B δ out , where we will use adiabatic theory with a non-zero gap larger than a constant times δ. In the sets B δ i , where the gap closes, we study the structure of the Taylor expansion of the Hamiltonian H k close to the Dirac points and construct two-band reduced Hamiltonians H R i,k . Then, we use the linear response Proposition 5.7, reducing successively from the Hamiltonian H to the reduced Hamiltonian H R i,k , and finally to the Dirac Hamiltonian

H D k = 0 k 1 -ik 2 k 1 + ik 2 0
for which we can explicitly compute the current. Adding the contributions, we will obtain

σ α,β := lim t→∞ lim ε→0 1 εt ˆt 0 j ε α,β (t )dt = |I| 16 e α • e β + O(δ)
Finally, we will pass to the limit δ → 0.

Far from the Dirac points

We set j ε,out α,β (t ) := -

1 4π 2 ˆBδ out Tr ∂ α H k U ε β,k (t ) γ 0 k U ε β,k (t ) * dk.
Let k ∈ B δ out . In the regime we consider, γ 0 k-εe β t = P Nsm,k-εe β t is gapped with a gap larger than a constant times δ. Applying the analysis of the previous sections, we obtain that

-Tr ∂ α H k U ε β,k (t) γ 0 k U ε β,k (t) * = -∂ α Tr H k-εe β t γ 0 k-εe β t -iε Tr ∂ α H k-εe β t L + k-εe β t ∂ β γ 0 k-εe β t + iε Tr ∂ α H k-εe β t U ε β,k (t) L + k ∂ β γ 0 k U ε β,k (t) * + O ε 2 tδ -6 .
We treat each term separately.

• For the first (adiabatic) term, we have

-∂ α Tr H k-εe β t γ 0 k-εe β t = -∂ α Tr H k γ 0 k + εt∂ α ∂ β Tr H k γ 0 k + O ε 2 t 2 δ -4 .
By time-reversal symmetry, the first term vanishes when integrated on B δ out . Using Stokes formula for the second term as in the metallic case, we get ˆBδ

out -∂ α Tr H k-εe β t γ 0 k-εe β t dk = εt n≤Nsm ˆ∂B δ in ∂ α λ n,k (ds • e β ) + O ε 2 t 2 δ -4 (84)
• For the second (static) term we similarly expand in ε -iε

ˆBδ out Tr ∂ α H k-εe β t L + k-εe β t ∂ β γ 0 k-εe β t dk = -iε ˆBδ out Tr γ 0 k [∂ α γ 0 k , ∂ β γ 0 k ] dk + O(ε 2 tδ -6 ) = O(ε 2 tδ -6 ),
where we used the fact that the function k

→ Tr γ 0 k [∂ α γ 0 k , ∂ β γ 0 k ] is odd.
• For the third (oscillatory) term, we use the same arguments as in the insulating case, and obtain

1 t ˆt 0 dt ˆBδ out iε Tr ∂ α H k-εe β t U ε β,k (t ) L + k ∂ β γ 0 k U ε β,k (t ) * dk = O εδ -6 1 t + ε .
We are left with

1 t ˆt 0 j ε,out α,β (t ) dt = ε 4π 2 t ˆt 0 t dt n≤Nsm ˆ∂B δ in ∂ α λ n,k (ds • e β ) + O ε δ 6 1 t + ε(1 + t 2 ) = εt 2π 2 n≤Nsm ˆ∂B δ in ∂ α λ n,k (ds • e β ) + O ε δ 6 1 t + ε(1 + t 2 ) . ( 85
)
8.2 Close to Dirac points: reduction to the 2-band case

We set j ε,in α,β (t) := -

1 4π 2 i∈I ˆBδ in Tr ∂ α H k U ε β,k (t) γ 0 k U ε β,k (t) * dk. (86)
Using the linear response Proposition 5.7, we have for almost all k ∈ B δ in and all ε, t ≥ 0,

-Tr ∂ α H k-εe β t U ε β,k (t) γ 0 k U ε β,k (t) * = -∂ α Tr H k γ 0 k + εt∂ α ∂ β Tr H k γ 0 k -iε Tr ∂ α H k (e -itL k -1)L + k ∂ β γ 0 k + O(ε 2 t 3 (1 + t 3 )).
Using (72), the equality Tr H k γ 0 k = Nsm n=1 λ n,k and Assumption 2.6, it is easily seen that the left-hand side, as well as the first and fourth terms of the right-hand side of that equation, are bounded uniformly in k and therefore integrable on B δ i . Besides, for k ∈ B δ in , the second term is bounded by a constant multiple of (1/|k -k i |) as k → k i , and is therefore integrable. It follows that the third term is also integrable on B δ i . We treat the three leading terms of the right-hand side in sequence.

• The first term vanishes when integrated on the time-reversal symmetric set ∪ i∈I B δ i . • For the second, arguing as in the metallic case, we get ˆBδ

i εt∂ αβ Tr H k γ 0 k dk = -εt n≤Nsm ˆ∂B δ i ∂ α λ n,k (ds • e β ),
so that the corresponding term in (86) cancels the contribution (85) from B δ out . • For the third term, we use

Tr ∂ α H k (e -itL k -1)L + k ∂ β γ 0 k = Tr ∂ α H k (e -itL k -1)(L + k ) 2 [γ 0 k , ∂ β H k ] = n≤Nsm m>Nsm (e -it(λ n,k -λ m,k ) -1) u n,k , ∂ β H k u m,k u m,k , ∂ α H k u n,k (λ m,k -λ n,k ) 2 -c.c.
with the sum converging from the asymptotics [START_REF] Freund | Peierls substitution for magnetic Bloch bands[END_REF].

When n = N sm or m = N sm + 1, the denominators in that equation are bounded from below independently of δ. The constant term vanishes when integrated over the time-reversal symmetric set ∪ i∈I B δ i , and the oscillatory term can be treated using the formula

1 t ˆt 0 e -iωt dt = e -iωt -1 -iωt
with ω = λ n,k -λ m,k bounded away from zero independently of δ.

Putting all the results of the previous two sections together, we get that

1 t ˆt 0 j ε α,β (t )dt = 1 t ˆt 0 j ε,out α,β (t )dt + 1 t ˆt 0 j ε,in α,β (t )dt = - iε 4π 2 t i∈I ˆt 0 dt ˆBδ i (e -it (λ Nsm ,k -λ N sm+1 ,k ) -1) × u Nsm,k , ∂ β H k u Nsm+1,k u Nsm+1,k , ∂ α H k u Nsm,k (λ Nsm+1,k -λ Nsm,k ) 2 dk -c.c. + O ε 2 δ -6 (1 + t 2 ) + εδ -6 t -1 + ε 2 t 3 (1 + t 3 ) .
At this stage of the proof, only two modes are involved in the formula giving the current, namely the two modes that cross at the Fermi level. Everything happens as for a two-band model that we now study. We write for short

1 t ˆt 0 j ε α,β (t ) ε dt = - i 4π 2 t i∈I ˆt 0 I R,i α,β (t ) dt + O εδ -6 (1 + t 2 ) + δ -6 t -1 + εt 3 (1 + t 3 ) , with I R,i α,β (δ, t) = ˆBδ i (e -it(λ Nsm,k -λ N sm+1 ,k ) -1) u Nsm,k , ∂ β H k u Nsm+1,k u Nsm+1,k , ∂ α H k u Nsm,k (λ Nsm+1,k -λ Nsm,k ) 2 dk -c.c..

Close to the Dirac points: the local model

We now are interested in the computation of I R,i α,β (t). In the following, we drop the index i and assume without loss of generality that k i = 0.

Hypothesis [START_REF] Michael | Methods of Modern Mathematical Physics. IV. Analysis of operators[END_REF] implies that, for k small enough, the Bloch Hamiltonian H k has exactly two eigenvalues (counting multiplicities) close to µ F . Consider an arbitrary orthonormal basis (v 0 , w 0 ) of Ran(P Nsm+1,0 -P Nsm-1,0 ). For all k small enough, we can construct an orthonormal basis (v k , w k ) of Ran(P Nsm+1,k -P Nsm-1,k ) by Löwdin orthonormalization of ((P Nsm+1,k -P Nsm-1,k )v 0 , (P Nsm+1,k -P Nsm-1,k )w 0 ), and set

H R k = [v k |w k ] * H k [v k |w k ] = v k , H k v k v k , H k w k w k , H k v k w k , H k w k .
It follows that the reduced Hamiltonian

H k Ran(P Nsm+1,k -P N sm-1 ,k )
is equivalent through a unitary transform that depends analytically on k to the reduced 2 × 2 Hamiltonian

H R k = 3 p=0 b p (k)σ p , (87)
where where the (q p ) p=1,2,3 are the rows of a 3 × 2 matrix Q with orthogonal columns, so that

σ 0 = Id C 2 , σ 1 = 0 1 1 0 , σ 2 = 0 -i i 0 , σ 3 = 1
H R k = µ F + v F (Qk) • σ + O(|k| 2 ). ( 88
)
Let R ∈ SO(3) be a rotation matrix that maps Ran(Q) to Span(e 0 1 , e 0 2 ), where (e 0 1 , e 0 2 , e 0 3 ) is the canonical basis of R 3 . Let U be one of its associated 2 × 2 unitary matrices through the two-to-one SU(2) → SO(3) mapping, so that [START_REF] Cornwell | Group Theory in Physics[END_REF] R pq = 1 2 Tr(σ p U σ q U * ).

It follows that

Tr Up to a unitary transform, we can therefore assume Q to be a 2 × 2 matrix in (88). for some phases θ ± (k) ∈ R. We have

σ 3 U H R k U * = v F 3
∂ α H R k = [v k |w k ] * ∂ α H k [v k |w k ] + ∂ α [v k |w k ] * [v k |w k ]H R k + H R k [v k |w k ] * ∂ α [v k |w k ] = [v k |w k ] * ∂ α H k [v k |w k ] + O(|k|)
where we have used for the first line that Since (λ R -,k -λ R +,k ) is bounded from below by a constant multiple of |k|, it follows that

I R α,β (δ, t) = ˆB(0,δ) (e -it(λ R -,k -λ R +,k ) -1) u R -,k , ∂ β H R k u R +,k u R +,k , ∂ α H R k u R -,k (λ R +,k -λ R -,k ) 2 dk -c.c. + O(δ). ( 89
)
Reduction to H D k By standard results of perturbation theory [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF] applied to H R k = µ F + v F (Qk) • σ + O(|k| 2 ) with gap greater than a constant multiple of |k|,

λ R ±,k = λ Q ±,k + O(|k| 2 ) and u R ±,k = u Q ±,k + O(|k|)
where the superscript Q refers to eigenvalues and appropriately chosen orthonormal eigenvectors of the Hamiltonian

H Q k = v F (Qk) • σ.
It follows that

I R α,β (δ, t) = I Q α,β (δ, t) + O(δ),
where I R α,β (δ, t) is defined similarly to (89) as

I Q α,β (δ, t) = ˆB(0,δ) (e -it(λ Q -,k -λ Q +,k ) -1) u Q -,k , ∂ β H Q k u Q +,k u Q +,k , ∂ α H Q k u Q -,k (λ Q +,k -λ Q -,k ) 2
dk -c.c..

We perform the change of variable k = Qk (recall that Q is orthogonal) and obtain

I Q α,

A Proofs of two technical lemmata

A.1 Proof of Lemma 4.1

Proof. We replicate the proof of the Faris-Lavine Theorem given in [START_REF] Michael | Methods of Modern Mathematical Physics. II. Fourier analysis, self-adjointness[END_REF], replacing the Laplacian by Hypotheses 1 and 3 guarantee that A + b + c|x| 2 is bounded below. Hence, increasing the value of b if necessary to make this operator positive, we have

(A + b)φ 2 L 2 ≤ N φ 2 L 2 + 4cd φ 2 L 2 ,
which proves (90). A.2 Proof of Lemma 4.2

Proof. By the Kato-Rellich theorem, for any 0 ≤ t ≤ T , H(t) is self-adjoint on L 2 per with domain H 2 per , and bounded below. We will show that there exists µ > 0 so that the graph norm of (H(t) + µ) for any 0 ≤ t ≤ T is equivalent to the H 2 per -norm. This will prove Lemma 4.2 by Proposition 2.1 in [START_REF] Teufel | Adiabatic Perturbation Theory in Quantum Dynamics[END_REF] (see also Theorem X.70 in [START_REF] Michael | Methods of Modern Mathematical Physics. II. Fourier analysis, self-adjointness[END_REF]).

We have for any µ > 0, 0 ≤ t ≤ T and φ ∈ H and so the graph norm is controlled by the H 2 per -norm. For the other inequality, we relate the resolvent of H(t) to that of H 0 by a bounded operator, with bounded inverse. Notice that, for any µ > 0, since H 0 is positive, ∀ 0 ≤ t ≤ T, (H(t) + µ) = (1 + H 1 (t)(H 0 + µ) -1 )(H 0 + µ).

Furthermore,

∀ 0 ≤ t ≤ T, H 1 (t)(H 0 + µ) -1 ≤ a H 0 (H 0 + µ) -1 + b (H 0 + µ) -1 ≤ a + b µ .
and so, for µ > b 1-a , the operator 1 + H 1 (t)(H 0 + µ)

  R * -periodic, and we have Tr(A) = (2π) -d ˆB Tr(A k ) dk.

Figure 1 :

 1 Figure 1: Commutative diagrams of the relationships between density matrices γ ε β and γ ε β (top) and the fibers γ ε β,k of γ ε β which decompose both γ ε β and γ ε β (bottom). In the top diagram, A U → B means that B = U AU * . In

Figure 2 : 1 ε

 21 Figure 2: Instantaneous conductivity j ε α,β (t) ε (solid line) and running average 1 t ´t 0 j ε α,β (t ) ε dt (dotted line) for several phases, in the linear response regime (ε = 10 -4 , t 1 ε ). In all cases e β = b 1 , and e α = e β , except in panel (b) where e α = b 2 .

Figure 3 : 5 2 b 2 .

 352 Figure 3: Instantaneous conductivity j ε α,β (t) ε in the Bloch oscillations regime (ε = 10 -2 , 1 ε t). We take e α = b . The left figure is with e β = b 1 + b 2 , the right one with e β = b 1 + 1+ √ 5 2 b 2 .

)

  H3 there exist M ∈ R + and bounded continuous functionsa ± : [0, T ] → R with 0 ≤ a -≤ a + ≤ M defining bounded closed intervals I(s) = [a -(s), a + (s)] ⊂ R such that, for all s ∈ [0, T ), g(s) := min (1, dist(I(s), σ(H(s)) \ (σ(H(s)) ∩ I(s)))) > 0, P (s) := 1 I(s) (H(s)) has a finite (constant) rank N ∈ N * ,Under these assumptions, we set L + (s) := L + H(s),I(s) . Proposition 5.1. Assume H1, H2 and H3. Then,P ∈ C n ([0, T ), L(H f )), L + ∈ C n ([0, T ), L(L(H f ))),andṖ (s) = L(s) + [P (s), Ḣ(s)].(51)Furthermore, there exist constantsC 1 , • • • , C n ∈ R + depending only on α 1 , • • • , αn and M such that the following bounds hold for all 0 ≤ l ≤ n, s ∈ [0, T ) and A ∈ L(H f ):

0 0 - 1 3 p=1 b p (k) 2 .

 132 are the Pauli matrices, and (b p ) p=0,1,2,3 are real-valued analytic functions of k in a neighborhood of 0.The matrix H R k has eigenvaluesλ ± (k) = b 0 (k) ± It follows that b 0 (k) = µ F + O(|k| 2 ), b p (k) = v F q p , k + O(|k| 2 ), p = 1, 2, 3,

q=1( 3 ,

 3 Qk) q Tr(σ 3 U σ q U * ) + O(|k| 2 ) RQk + O(|k| 2 ) = O(|k| 2 ).

8. 4

 4 The two-band case: reduction to the Dirac Hamiltonian Reduction to H R k For k = 0, let λ R ±,k be the larger and smaller eigenvalues of H R k respectively, and u R ±,k associated orthonormal eigenvectors in C 2 . We have λ R -,k = λ Nsm,k , λ R +,k = λ Nsm+1,k , and[v k |w k ]u R +,k = e iθ+(k) u Nsm+1,k , [v k |w k ]u R -,k = e iθ-(k) u Nsm,k

  H k commutes with P Nsm+1,k -P Nsm,k = [v k |w k ][v k |w k ] * ,and for the second thatH R k = µ F Id 2 + O(|k|) and ∂ α [v k |w k ] * [v k |w k ] = ∂ α Id 2 = 0. We therefore obtain u R +,k , ∂ α H R k u R -,k = e -i(θ+(k)-θ-(k)) u Nsm+1,k , ∂ α H k u Nsm,k + O(|k|) u R -,k , ∂ β H R k u R +,k = e +i(θ+(k)-θ-(k)) u Nsm,k , ∂ β H k u Nsm+1,k + O(|k|).

2 -e -iθ 1 .

 21 β (δ, t) = e T β I D (δ, t)e α ,where the coefficients I D ij (δ, t) of the 2 × 2 matrix I D (δ, t) are given byI D ij (δ, t) = ˆB(0,δ) (e -it(λ D -,k -λ D +,k ) -1) u D -,k , ∂ kj H D k u D +,k u D +,k , ∂ ki H D k u D -,k (λ D +,k -λ D -,k ) 2 dk -c.c.and the superscript D refers to the Dirac HamiltonianH D k = v F k • σ.The Dirac Hamiltonian H D kWe finish by computing I D (δ, t) explicitly. Let k = r(cos θ, sin θ). We haveλ D ±,k = ±v F r, u D +,By an explicit calculation, we obtainˆ2π 0 (e -it(λ D -,k -λ D +,k ) -1) u D -,k , ∂ ki H D k u D +,k u D +,k , ∂ kj H D k u D -,k (λ D +,k -λ D -,k ) 2dθ -c.c. = iπ 1 2r

1 2 ( 1 2 φ 2 .

 212 -i∇ + A)2 . It consists in verifying the following two hypotheses of[START_REF] Michael | Methods of Modern Mathematical Physics. II. Fourier analysis, self-adjointness[END_REF] Theorem X.37]. Let A = 1 2 (-i∇ + A) 2 + W + V and N = A + 2c|x| 2 + b,where b ∈ R will be specified below: there exists h, such that for any φ ∈ C, Aφ ≤ h N φ ; (90) for some , for any φ ∈ C, |(Aφ, N φ) -(N φ, Aφ)| ≤ N (91) By hypothesis 3 in Lemma 4.1 and the conditions on W , it is possible to choose b so that N ≥ 1. As quadratic forms on C, N 2 = (A + b) 2 + 4c d j=1 x j (A + b + c|x| 2 )x j -2cd.

For ( 91 )

 91 , we observe that±i[A, N ] = ±2c(x • (-i∇ + A) + (-i∇ + A) • x) ≤ 2c (-i∇ + A) 2 + |x| 2 ≤ N,where we have used(-i∇ + A) 2 + |x| 2 ± (x • (-i∇ + A) + (-i∇ + A) • x) = (-i∇ + A ± x) + A) 2 + V + (W + c|x| 2 ) + 1 -a 2 (-i∇ + A) 2 + c|x| 2 + b ≥ e((-i∇ + A) 2 + |x| 2 ),where e = min(c, 1-a 2 ) > 0 and where b is chosen so thatb -f + min σ a 2 (-i∇ + A) 2 + V ≥ 0.This proves (91). Hence A is essentially self-adjoint on C.

  t) + µ)φ L 2 per ≤ (1 + a) H 0 φ L 2 per + (b + µ) φ L 2 per ≤ (1 + a + b + µ) φ H 2 per ,

  It follows that there exists a constant C depending only on V , A and µ F , such that|ρ ε k (t)| ≤ Cε 2 t 3 t + e ηεt/2 1 + t(1 + (εt) 2 ) + (εt)t 2 (1 + (εt) 3 ) + (εt) 4 t 3 (1 + (εt) 2 ) ,which leads to (78) when εt ≤ 1.

	3 λe ηεt/2 1 6	+ ν 1 ν 2	t 8	+ (εt)t	ν 1 30	+	ν 2 20	+ (εt) 2 t	1 72	.

  2 sin(2v F rt)δ ij .

	It follows that
	1 t ˆt 0	I D ij (δ, t ) dt = =	iπ 2 iπ 4	δ ij δ ij	1 t ˆt 0 ˆδ 0 ˆδv F t 0 1 -cos(2r ) sin(2v F rt ) r (r ) 2 dr = dr dt = iπ 2 4	iπ 4t δ ij + O (δt) -1 . δ ij ˆδ 0 1 -cos(2v F rt) v F r 2	dr
	We finally get by summing all the estimates
	1 t ˆt 0	j ε α,β (t ) ε	dt =	|I| 16	e α • e β + O δ + εt 3 (1 + t 3 ) +	1 δ 6	1 t	+ ε(1 + t) 2	,
	hence the result.

  -1 is bounded and invertible with bounded inverse in L 2 per . Therefore (H(t) + µ) -1 is bounded from L 2 per to H 2 per , which means there exists C > 0 such that, for any φ ∈ H 2 per and 0 ≤ t ≤ T , φ H 2 per = (H(t) + µ) -1 (H(t) + µ)φ H 2 per ≤ C (H(t) + µ)φ L 2 per , which concludes the proof.

Acknowledgements. We are grateful to Caroline Lasser for stimulating discussions all along the preparation of this article. We thank the two reviewers for their constructive suggestions, and in particular for pointing out to us a more elegant proof of Proposition 2.1 and the extension to uniform magnetic fields discussed in Remark 2.3. This project has been supported by Labex Bezout and has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 810367). The second author thanks the mathematics department of the Technische Universität München for hosting her during the final writing of this article.