Adiabatic Evolution of Low-Temperature Many-Body Systems

We consider finite-range, many-body fermionic lattice models and we study the evolution of their thermal equilibrium state after introducing a weak and slowly varying time-dependent perturbation. Under suitable assumptions on the external driving, we derive a representation for the average of the evolution of local observables via a convergent expansion in the perturbation, for small enough temperatures. Convergence holds for a range of parameters that is uniform in the size of the system. Under a spectral gap assumption on the unperturbed Hamiltonian, convergence is also uniform in temperature. As an application, our expansion allows us to prove closeness of the time-evolved state to the instantaneous Gibbs state of the perturbed system, in the sense of expectation of local observables, at zero and at small temperatures. As a corollary, we also establish the validity of linear response. Our strategy is based on a rigorous version of the Wick rotation, which allows us to represent the Duhamel expansion for the real-time dynamics in terms of Euclidean correlation functions, for which precise decay estimates are proved using fermionic cluster expansion.


Introduction
Adiabatic quantum dynamics. The adiabatic theorem is a fundamental result in quantum mechanics, dating back to the work of Born and Fock [13] and Kato [34]. Let us review its basic statement. Let H (s) be a family of time-dependent Hamiltonians, depending smoothly on the time parameter s for s ∈ [−1, 0]. We shall suppose that H (s) has a unique ground state ϕ s , and that the energy of the ground state is separated from the rest of the spectrum by a positive gap. We are interested in the adiabatic regime, defined as follows. Let η > 0 and consider the time-dependent Schrödinger equation: The constant C in Eq. (1.2) depends on the details of the model, in particular on the regularity of H (s). The way in which the regularity of H (s) is quantified is typically via an estimate for Ḣ (s) . This quantity is badly behaved in situations in which the Hamiltonian describes a many-body system, say an interacting Fermi gas on a lattice Λ L = [0, L] d ∩ Z d , due to the fact that the norm of the Hamiltonian and of its derivatives grow linearly with the size of the system. Thus, the standard adiabatic theorem fails in describing the evolution of many-body quantum systems for η small uniformly in L.
This is not a technical point. In fact, it turns out that the notion of convergence in Eq. (1.2) is not the natural one for many-body systems: one cannot expect norm convergence for extensive many-body quantum systems uniformly in their size, see for instance the discussion in [7]. Instead, a much more natural notion of convergence involves the expectation value of local observables, which only probe a finite region in space. In this setting, a many-body adiabatic theorem for quantum spin systems has been recently proved in [6]. The result has then been extended to Fermi gases in [40]. Specifically, let H (s) be a timedependent Hamiltonian for a quantum spin systems, or for lattice fermions, on a large but finite lattice Λ L ⊂ Z d . Suppose that H (s) has a spectral gap for all times s ∈ [−1, 0], and let Π L (s) be the projector associated with the ground state of H (s) on Λ L . Let P L (t ) be the solution of the evolution equation: Consider the expectation value of a local operator on the time-dependent state, Tr O X P L (t ). Then, under reasonable regularity and locality assumptions on the Hamiltonian, the many-body adiabatic theorem states that [6,40]: Tr O X P L (t ) − Tr O X Π L (ηt ) ≤ C η for all t ∈ [−1/η, 0], (1.4) where the constant C depends on the observable O X , but it is independent of L. An important application of this result is the proof of validity of linear response for extended, many-body quantum systems. To introduce the notion of linear response, let us further assume that the many-body Hamiltonian has the form: H (ηt ) = H + εg (ηt )P (1.5) where H and P are given by sums of local potentials, and g (t ) is a switch function, that is a bounded function that decays fast enough at negative times. A standard choice is the exponential switch function, g (t ) = e t . Consider the dynamics generated by (1.5) for t ∈ (−∞, 0]. Let P L (t ) be the solution of (1.3) with initial datum P L (−∞) = Π L (−∞). Then, [6,40] proved that, see also the reviews [7,29]: where χ O ,P agrees with the well-known Kubo formula for linear response. The statement (1.6) holds provided the thermodynamic limit of P L (−∞) exists. In general, a similar result holds for ε, η, L fixed, where ε, η are small uniformly in L and where ε is small uniformly in η. The Kubo formula is equivalent to the first order term in the Duhamel expansion for the non-autonomous evolution: These results have interesting applications in condensed matter physics. In particular, combined with [28,24,9], they allow to prove the quantization of the Hall conductivity for gapped interacting fermions starting from the fundamental many-body Schrödinger equation. Among other extensions that have been obtained in the last few years we mention: the construction of non-equilibrium almost-stationary states and the application to the proof of validity of linear response for a class of perturbations that might close the spectral gap [46]; the proof of exactness of linear response for the quantum Hall effect [8]; the extension of the many-body adiabatic theorem to infinite systems with a bulk gap [30]. Despite all this progress, an important limitation of the existing approaches is that they do not allow to study many-body quantum systems at positive temperature. In particular, the zero temperature limit is taken before the thermodynamic limit. It is of obvious physical relevance to consider the situation in which the thermodynamic limit is taken at fixed positive temperature, to make contact with experimental settings in which the temperature is possibly small but necessarily non-zero. In what follows, we will focus on interacting lattice fermionic models, which we shall describe in the grand-canonical Fock space formalism. We are interested in the following evolution equation: with o(1) a quantity that vanishes as η → 0 + , uniformly in L.
Our result. In this work, we introduce a different approach to study manybody quantum dynamics in the adiabatic regime, that also applies to many-body systems at positive temperature. We consider finite-range, time-dependent Hamiltonians of the form (1.5), under suitable assumptions discussed below. In our main result, Theorem 3.7, we derive a representation of Tr O X ρ(t ) via a convergent expansion in ε, uniformly in η and in L, for small temperatures. "Small" means that: Under a gap condition on the many-body Hamiltonian, the range of allowed ε for which convergence holds is also uniform in β. Our method allows to prove the validity of an adiabatic theorem for local observables, in the form of Eq. (1.9), for small temperatures in the sense of (1.10). This allows, in particular, to recover the usual zero-temperature many-body adiabatic theorem, by taking the infinite volume limit before the zero temperature limit. Furthermore, the method also allows to prove the validity of linear response, and more generally to compute all higher-order response coefficients in terms of equilibrium correlations, see Corollary 3.9. The proof is based on a rigorous Wick rotation, that allows to rewrite the Duhamel expansion for the quantum evolution of the system in terms of timeordered, connected Euclidean (or imaginary-time) correlation functions. Previously, this idea has been used to study rigorously the linear response and the quadratic response in a number of interacting gapped or gapless systems [24,4,27,39]. Here, we extend this strategy at all orders in the Duhamel expansion for the time-evolution of the state, and we use it to prove convergence of the Duhamel series for the real-time dynamics.
Our method applies to a class of switch functions g (ηt ) that can be approximated, for β large, by functions g β,η (t ) such that g β,η (t ) = g β,η (t −i β). This identity plays a key role in the proof of the Wick rotation. This requirement of course restricts the class of switch functions that we are able to consider; however, let us anticipate that this approximation holds for the standard exponential switch function, and more generally for the Laplace transform of suitable integrable functions.
Theorem 3.7 holds under a decay assumption on the Euclidean correlations of the Gibbs state of H , which is satisfied by finite-range Hamiltonians of the form H = H 0 + λV (1.11) with H 0 the second-quantization of a gapped Hamiltonian, V the many-body interaction and |λ| small. This is precisely the setting considered in [24], for the universality of the Hall conductivity, or in [19], where the stability of the spectral gap for the many-body Hamiltonian has been proved. This class of weakly interacting models can be analyzed via fermionic cluster expansion techniques, which allow to prove essentially optimal estimates for the decay of the Euclidean correlation functions. The assumptions on the Euclidean correlations actually hold even without a gap condition, however in this case one is forced to consider a range of ε that shrinks as T → 0 (but it is still uniform in L).
Besides the actual result, we believe that a relevant contribution of the present work is to import methods developed for interacting fermionic models at equi-librium to the study of real-time quantum dynamics. In perspective, if combined with rigorous renormalization group techniques, see [11,44,37] for reviews, we think that the approach of this paper could be extended to study the evolution of the Gibbs state of metallic or semimetallic systems, where the Fermi energy of the initial datum lies in the spectrum of the Hamiltonian. There, one does not expect an adiabatic theorem to hold; however, one might still have a convergent series expansion for the expectation of local observables in terms of Euclidean correlations, in a physically relevant range of parameters. This would be useful to establish the validity of linear response for gapless systems, widely used in applications.
Specifically, the combination of cluster expansion with rigorous renormalization group recently allowed to study the low temperature properties of a wide class of interacting gapless systems, and in particular to access their transport coefficients defined in the framework of linear response. Among the recent works, we mention the construction of the ground state of the two-dimensional Hubbard model on the honeycomb lattice [22] and the proof of universality of the longitudinal conductivity of graphene [23]; the construction of the topological phase diagram of the Haldane-Hubbard model [25,26]; the proof of the nonrenormalization for the chiral anomaly of Weyl semimetals [27]; the proof of Luttinger liquid behavior for interacting edge modes of two-dimensional topological insulators and the proof of universality of edge conductance [4,38,39]. It would be very interesting to prove the validity of linear response in the setting considered in these works, starting from many-body quantum dynamics.
Furthermore, it would be interesting to extend the methods presented in this work in the direction of studying spin transport, and prove the validity of linear response for spin-noncommuting many-body Hamiltonians. For noninteracting models, recent progress has been obtained in [36,35].
The adiabatic evolution of positive temperature quantum systems has been studied in the last years, e.g. in [1,2,31], see [32,33] for more recent developments. The setting considered in these works is however different from the one of the present paper. In [1,2,31] the Authors consider a small system coupled to reservoirs, and study the dynamics of the small system when the coupling with the reservoirs is introduced slowly in time. The key technical tool introduced in [1,2,31] is an isothermal adiabatic theorem, that proves norm-convergence of the evolved equilibrium state to the instantaneous equilibrium state of the perturbed system, in the adiabatic limit. The result holds under a suitable ergodicity assumption, which as far as we know, has not been proved for the class of extended, interacting Fermi systems considered here.
Ideas of the proof. Let us give a few more details about the method introduced in this paper. The proof starts by approximating the real-time dynamics generated by H (ηt ) by a suitable auxiliary dynamics, obtained from H (ηt ) by replacing the switch function g (ηt ) by a function g β,η (t ) such that lim β→∞ g β,η (t ) = g (ηt ) and g β,η (t ) = g β,η (t − i β). This approximation of course introduces an error, whose influence on the expectation values of local observables is estimated via Lieb-Robinson bounds for non-autonomous quantum dynamics [12,15]. This error is responsible for the limitation in the range of temperatures that we are able to consider. The big advantage of replacing g (ηt ) with g β,η (t ) is that the Duhamel series in ε for the auxiliary evolution can be written exactly in terms of Euclidean correlations, using the Wick rotation. This is ultimately possible thanks to the validity of the Kubo-Martin-Schwinger (KMS) identity for the thermal expectation of observables of the form g β,η (t )e i H t O X e −i H t . Then, once the Duhamel series is represented in terms of Euclidean correlations, the convergence of the series follows from the good decay properties of Euclidean correlations, which can be proved via cluster expansion techniques. Finally, the connection with the instantaneous Gibbs state of H (ηt ) is obtained by noticing that, for η small, the Wick-rotated Duhamel series agrees with the equilibrium perturbation theory in ε for the Gibbs state of the Hamiltonian H + εg (ηt )P .
An important ingredient of our proof is the complex deformation argument of Propositions 4.4, 4.5, which allows to prove the Wick rotation at all orders in the Duhamel series. Propositions 4.4, 4.5 are the adaptation of Propositions 5.4.12, 5.4.13 of [14] to our adiabatic setting. The main difference with respect to [14] is that in our case the observables involved in the correlations are "damped" in time by g β,η (t ): this allows to rule out the presence of spurious boundary terms at infinity in the complex deformation argument. In [14], these boundary terms are controlled by a suitable clustering assumption on the real-time correlation functions of the equilibrium state, which are very hard to prove for interacting models in the infinite volume limit. We are not aware of any result in this direction for the class of many-body lattice models considered in the present work.
Structure of the paper. The paper is organized as follows. In Section 2 we introduce the class of models considered in this work, we introduce their Gibbs state and Euclidean correlation functions, and we define the quantum dynamics. In Section 3 we state our main result, Theorem 3.7, which provides a representation for the average of the real-time evolution of local observables via a convergent expansion in ε. As an application, this representation establishes a many-body adiabatic theorem for the evolution of thermal states at low temperature. A relevant consequence of the proof of our main result is the validity of linear response, Corollary 3.9. The proof of the main result will be given in Section 4. In Appendix A we further discuss the class of switch functions considered in the present work; in Appendix B we discuss some well-known properties about time-ordered Euclidean correlations, which we include for completeness; and in Appendix C we review the proof of our Assumption 3.1, which is known to hold for many-body systems at zero temperature and with a spectral gap, or at positive temperature. This is done using fermionic cluster expansion, whose convergence is guaranteed by the celebrated Brydges-Battle-Federbush-Kennedy formula for cumulants.

The model
In this section we define the class of models we shall consider in this paper. We will focus on lattice fermionic systems, with finite-range interactions. We will then define the time-evolution of such systems, after introducing a timedependent perturbation.
Remark 2.1. Unless otherwise specified, the constants C , K etc. appearing in the bounds do not depend on β, L, η, ε and on time. Their values might change from line to line. Also, it will be understood that the natural numbers N include zero.

Lattice fermions
Let Γ be a d -dimensional lattice, namely where a 1 , . . . , a d are d linearly independent vectors in R d . Let L ∈ N, L > 0. We define the lattice dilated by L as LΓ := Span LZ {a 1 , . . . , a d } ∼ = LZ d . The finite torus of side L is defined as Γ L := Γ/(LΓ), that is: n i a i n i ∈ Z, 0 ≤ n i < L with periodic boundary conditions. The Euclidean distance on the torus Γ L is given by We shall denote by M ∈ N, M > 0 the total number of internal degrees of freedom of a particle. This might take into account the spin degrees of freedom, or sublattice labels. Setting S M := {1, . . . , M }, we define: We equip Λ L with the following distance, tracing only the space coordinates. For any x = (x, σ), y = (y, σ ) ∈ Λ L , we define: We shall describe fermionic particles on Λ L , in a grand-canonical setting. To this end, we introduce the fermionic Fock space, as follows. Let the one-particle Hilbert space be h L := 2 (Λ L ). The corresponding N -particle Hilbert space is its N -fold anti-symmetric tensor product H L,N := h ∧N L ; notice that the antisymmetric tensor product is trivial whenever N > M L d . The fermionic Fock space is defined as usual: where h L,0 := C.
For finite L, the fermionic Fock space is a finite-dimensional vector space. Thus, any linear operator on F L into itself is automatically bounded, and can be represented as a matrix. For any x ∈ Λ L , let a − x and a + x be the standard fermionic annihilation and creation operators, satisfying the canonical anti-commutation relations: {a For any subset X ⊆ Λ L , we denote by A X the algebra of polynomials over C generated by the fermionic operators restricted to X , {a − x , a + x : x ∈ X }. An example of operator in A Λ L is the number operator, defined as: The operator N counts how many particles are present in a given sector of the Fock space: given ψ ∈ F , it acts as N ψ = (0ψ (0) , 1ψ (1) , . . . , nψ (n) , . . .).
We shall denote by A N X the subset of A X consisting of operators commuting with N , also called gauge-invariant operators. Equivalently, these operators consist of polynomials in the creations and annihilation operators where the number of creation operators equals the number of annihilation operators.
Finally, let us define the notion of finite-range operators. Given X ⊆ Λ L , the diameter of X is defined as:

Definition 2.2 (Finite-range operators).
We say that O ∈ A Λ L is a finite-range operator if the following holds true.
(i) The operator O can be written as: (ii) There exists R > 0 independent of L such that O X = 0 whenever diam(X ) > R. Furthermore, there exists a constant C > 0 independent of X such that:

Dynamics
Hamiltonian and Gibbs state. For X ⊆ Λ L , an interaction potential Φ X is a selfadjoint operator in A N X . Being Φ X even in the number of creation and annihilation operators, it follows that if X ∩ Y = then [Φ X , Φ Y ] = 0. Given a sequence (Φ X ) X ⊆Λ L , we define the Hamiltonian H ∈ A N Λ L as: The Heisenberg time-evolution of an observable O ∈ A Λ L is, for t ∈ R: Later, we will also consider the Heisenberg evolution for complex times t , whose definition poses no problem due to the finite-dimensionality of the Hilbert space.
In the following, we shall suppose that H is a finite-range operator, in the sense of Definition 2.2. An example of Hamiltonian satisfying the above assumptions is: with λ ∈ R, |λ| small in a sense to be made precise, and V finite-range. We shall say that the non-interacting Hamiltonian H 0 is gapped if the spectrum of H has a spectral gap, uniformly in L. Given β > 0, µ ∈ R, the grand-canonical equilibrium state 〈·〉 β,µ,L associated with the Hamiltonian H , also called equilibrium Gibbs state, is defined as: where the trace is over the fermionic Fock space F L . Obviously, the Gibbs state is invariant under time evolution: It will also be convenient to define the imaginary-time, or Euclidean, evolution of O as: Notice that the imaginary-time evolution is no longer unitary, and the norm of γ t (O ) might grow in time. Finally, the following property holds true, also called KMS identity: For finite L, which is our case, this identity simply follows from the definition of Gibbs state, and from the cyclicity of the trace. Notice that the dynamics γ t in (2.7) trivially extends to all complex times t ; thus, the identity (2.8) actually holds replacing t 1 , t 2 by any two complex numbers z 1 , z 2 . Eq. (2.8) will play a fundamental role in our analysis.

Time ordering.
Let a x,t := γ t (a x ), where a x can be either a x or a * x . Let t 1 , . . . , t n in [0, β). We define the time-ordering of a 1 x 1 ,t 1 , . . . , a n x n ,t n as: Ta 1 x 1 ,t 1 · · · a n x n ,t n = (−1) π 1(t π(1) ≥ · · · ≥ t π(n) )a π (1) x π(1) ,t π(1) · · · a π(n) x π(n) ,t π(n) , (2.9) where π is the permutation needed in order to bring the times in a decreasing order, from the left, with sign (−1) π , and 1(condition) is equal to 1 if the condition is true or 0 otherwise. In case two or more times are equal, the ambiguity is solved by putting the fermionic operators into normal order. Other solutions of the ambiguity are of course possible; it is worth anticipating that in our applications this arbitrariness will play no role, since it involves a zero measure set of times. The above definition extends to operators in A Λ L by linearity. In particular, for O 1 , . . . , O n even in the number of creation and annihilation operators, we have: (1) ) · · · γ t π(n) (O π(n) ). (2.10) The lack of the overall sign is due to the fact that the observables involve an even number of creation and annihilation operators.

Euclidean correlation functions.
Let t i ∈ [0, β), for i = 1, . . . , n. Given operators O 1 , . . . , O n in A Λ L , we define the time-ordered Euclidean correlation function as: From the definition of fermionic time-ordering, and from the KMS identity, it is not difficult to check that: in the special case in which the operators involve an even number of creation and annihilation operators, which will be particularly relevant for our analysis, the overall sign is +1. The property (2.12) allows to extend in a periodic (sign +1) or antiperiodic (sign −1) way the correlation functions to all times t i ∈ R. From now on, when discussing time-ordered correlations we shall always assume that this extension has been taken, unless otherwise specified. Next, we define the connected time-ordered Euclidean correlation functions, or time-ordered Euclidean cumulants, as: where I is a non-empty ordered subset of {1, 2, . . . , n}, More generally, the following relation between correlation functions and connected correlation function holds true: where P is the set of all partitions of {1, 2, . . . , n} into ordered subsets, and J is an element of the partition P , J = { j 1 , . . . , j |J | }.
Driving the system out of equilibrium. We are interested in driving the system out of equilibrium, by adding a slowly varying time-dependent perturbation to the Hamiltonian. We define, for t ≤ 0: (2.14) where: H ∈ A N Λ L is a finite-range Hamiltonian, of the form (2.2); η > 0, ε ∈ R; g (·) is a smooth function vanishing at −∞, whose further properties will be specified later on; and P ∈ A N Λ L is a finite-range operator, As an example, we might consider: with µ(x) bounded uniformly in L. More generally, we will not require Ψ X to be quadratic in the fermionic operators.
The Hamiltonian H (ηt ) generates the following the Schrödinger-von Neumann non-autonomous evolution: We shall denote by U (t ; s) the unitary group generated by H (t ): (2.17) Using this unitary group, the solution of Eq. (2.16) can be written as Let O X ∈ A N X be a local operator. We will be interested in studying its expectation value on the time-dependent state: In particular, we will be interested in understanding the dependence of (2.19) in the external perturbation, and in establishing the validity of linear response, uniformly in the size of the system.

Main result
In what follows, we will consider Hamiltonians H (ηt ) = H + εg (ηt )P with H , P finite-range operators in A N Λ L , in the sense of Definition 2.2. We shall denote by 〈·〉 t the instantaneous Gibbs state of H (ηt ), Our main result holds under the following assumptions on the Hamiltonian H and on the switch function g (t ).
where: The next assumption specifies the class of switch functions g (t ) that we are able to consider. Assumption 3.3 (Properties of the switch function). We assume that g (t ) has the form, for all t ≤ 0: where the function h is such that: Otherwise, the function h(ξ) can be replaced by finite linear combination of Dirac delta distributions, supported on R + .

Remark 3.4. Thus, g is the Laplace transform of the function h. As discussed in Appendix A, the properties (3.5) are implied by suitable decay properties of the function g (z) for complex times. Our setting allows to include the function g (t ) = e t , a widely used switch function in applications, by choosing h(ξ)
Next, let us introduce a suitable approximation of the switch function, which will play an important role in our analysis. Definition 3.5 (Approximation of the switch function). Let η > 0 and suppose that g (t ) satisfies Assumption 3.3. We define: whereg β,η (0) := 0 and for ω ≥ 2π β : Remark 3.6.
(i) The approximation of the switch function satisfies the following key identity: The following estimate holds: (iii) Using that, for ξ 1 > ξ 2 and t ≤ 0, we have: Also, by Eq. (3.5): where the constant C > 0 only depends on the function h.
We are now ready to state our main result.

Remark 3.8.
(i) It is important to stress that ε 0 only depends on the constant c appearing in (3.2). In particular, it is independent of η and L. Furthermore, as commented after Assumption 3.1, for many-body perturbations of noninteracting gapped lattice models the constant c is actually independent of β: hence, for this class of models the range of ε for which Eq. (3.17) holds is independent of β as well.
(ii) For gapped systems, our result allows to take the zero temperature limit β → ∞ after the thermodynamic limit L → ∞. By Eqs. (3.13), (3.14), the existence of these limits is implied by the existence of the same limits for the equilibrium Gibbs state 〈·〉 β,µ,L . To the best of our knowledge, all previous works on many-body adiabatic dynamics considered the case in which the temperature is sent to zero before the thermodynamic limit.
(iii) Observe that in the limits β, L → ∞ and then η → 0 + , Eq. (3.17) implies that Tr O X ρ(t ) converges to the average over the zero temperature equilibrium state of the Hamiltonian H (0) = H + εg (0)P , thus recovering the zero temperature, many-body adiabatic theorem.
A rather direct consequence of the proof of the main result is the validity of linear response.

Corollary 3.9 (Validity of linear response).
Under the same assumptions of Theorem 3.7, the following is true: where the error term R β,µ,L (ε, η, t ) is bounded as: In Eq. (3.18), the function g β,η (t −i s) can be replaced by g (ηt ), up to replacing the error term R β,µ,L by R β,µ,L , such that: Furthermore, the main term in Eq. (3.18) is equal to the first order term in the Duhamel expansion, up to small errors: Remark 3.10. More generally, the method of proof allows to approximate all terms appearing in the Duhamel expansion for the real-time quantum dynamics in terms of the functions I (n) β,µ,L (η, t ) defined in (3.14).
The proof of the theorem is based on a precise control of the Duhamel expansion for the dynamics obtained after replacing g (ηt ) with g β,η (t ). As discussed in Remark 3.6, the lower the temperature, the better is the approximation of g (ηt ) in terms of g β,η (t ): the error term of order 1/β in our main result is produced by this replacement. The advantage of introducing this approximation is that the Duhamel expansion for the new dynamics can be expressed in terms of Euclidean correlation functions, for which precise decay estimates can be proved using tools from statistical mechanics, such as the cluster expansion. This is possible thanks to a rigorous version of the Wick rotation, at all orders in perturbation theory, which we prove by adapting an argument of [14]. With respect to [14], the presence of the switch function allows to avoid clustering assumptions on the real-time correlation functions of the system.
The proof of the main result will be given in Section 4, and it is organized as follows. In Section 4.1 we recall how to derive the Duhamel expansion for the many-body evolution. In Section 4.2 we introduce the auxiliary dynamics, obtained after replacing g (ηt ) with g β,η (t ) in H (ηt ), and we prove the closeness of the two dynamics in the sense of expectation of local observables using Lieb-Robinson bounds. In Section 4.3 we rewrite the Duhamel expansion for the auxiliary dynamics in terms of time-ordered Euclidean correlations, via the Wick rotation, and we use Assumption 3.1 to establish convergence of the (Wickrotated) Duhamel series, uniformly in the size of the system. In Section 4.4 we recall the cumulant expansion in ε for the instantaneous Gibbs state of H (ηt ). Finally, in Section 4.5 we put everything together, and we prove Theorem 3.7.

Proof of Theorem 3.7 4.1 Duhamel expansion
We start by recalling how to derive the well-known Duhamel series for the expectation of local observables. Given a time-dependent Hamiltonian H (t ) = H + εg (ηt )P , let us consider the associated unitary evolution: For ε = 0 one trivially has U (t ; s) = e −i H (t −s) . We are interested in deriving a perturbative expansion around the evolution generated by H . To this end, we define the unitary evolution in the interaction picture as: Clearly, U I (s; s) = 1, and: Next, we write, for T > 0 and for 0 ≥ t ≥ −T : where we used the cyclicity of the trace and the invariance of ρ β,µ,L under the dynamics generated by H . Finally, by Eq. (4.3): The procedure can be iterated. One gets: where R (m+1) β,µ,L (−T ; t ) is the Taylor remainder of the expansion, given by: On a finite lattice and for η > 0, the series is absolutely convergent. In fact, by using the boundedness of the fermionic operators, and the unitarity of the time evolution, we have the following crude estimate: for a universal constant C > 0. Thus, taking m large enough, uniformly in T , the error term can be made as small as wished. Hence, we have, in the T → ∞ limit: In order to extract useful information from this representation, we need estimates for the various terms that are uniform in the size of the system. In particular, we would like to prove that the series converges uniformly in ε, as L → ∞ and for η small. The main difficulty to achieve this is the control of the time-integral, uniformly in η. For fixed η, this problem might be approached using Lieb-Robinson bounds, see [41] for a review. This bound reads, for two operators O X and O Y supported on X , Y : for suitable positive constants C , c, v. Combined with the boundedness of the fermionic operators, this estimate (and its extension to multi-commutators, [15]) can be used to prove that the series in (4.9) is convergent uniformly in L, however only for |ε| ≤ ε(η) with ε(η) → 0 + as η → 0 + . In the next section, we shall study the time-evolution using a different approach, which gives estimates that are uniform in η.

The auxiliary dynamics
Let H (ηt ) := H + εg β,η (t )P , (4.11) with g β,η (t ) introduced in Definition 3.5. Here we will prove that the evolutions generated by H (ηt ) and by H (ηt ) are close, at small enough temperature, in the sense of the expectation of local observables. To compare the two evolutions, we will use Lieb-Robinson bounds for non-autonomous dynamics [12,15].
Proof. We start by writing: where g β,η (t ) is defined in Eq. (3.6), and the error term satisfies the bound, by Eq. (3.11): (4.14) Next, we write: where U (t ; s), U (t ; s) are the unitary groups generated by H (ηt ), H (ηt ), respectively. We estimate the argument of the limit as: where we used that ρ β,µ,L ≥ 0, Tr ρ β,µ,L = 1 and the unitarity of time-evolution. Next, we rewrite the argument of the norm as: where in the third line we used that U (t ; s) * = U (s; t ), and in the last line we used Eq. (4.13). Therefore, Next, we claim that: This inequality stems from the Lieb-Robinson bound for non-autonomous dynamics, see Theorem 4.6 of [12] for quantum spin systems, or Theorem 5.1 of [15] for the case of lattice fermions: for any two local operators O X , O Y . The proof of (4.19) is standard, and we give it here for completeness. Representing the perturbation P in terms of its local potentials, Eq. (2.15), we have: with D large enough to be chosen below. By the boundedness of the fermionic operators, and by the unitarity of the dynamics, the first term in the right-hand side is estimated as: where we used the fact that the sum is restricted to sets Y of bounded diameter. For the second term, we use the Lieb-Robinson bound (4.20), to get: Choosing D large enough, we have: This concludes the check of (4.19). Using the bound (4.19) in (4.18), we get: where in the last step we used the bound (4.14) and we exchanged the order of integration. Using that: we finally obtain: (4.27) Eq. (4.12) follows from assumption (3.5). This concludes the proof.

Wick rotation
Here we shall represent each coefficient in the Duhamel expansion (4.9) for the auxiliary dynamics generated by (4.11) in terms of Euclidean correlation functions, via a complex deformation argument. The advantage is that useful spacetime decay estimates for Euclidean correlations can be proved using statistical mechanics tools, such as the cluster expansion. This complex deformation is known in physics as Wick rotation, and here it will be established rigorously for the auxiliary dynamics. The next lemma is the main result of the section. Its proof is based on the adaptation of ideas of Section 5.4 of [14] to our adiabatic setting.

Lemma 4.2 (Wick rotation). Let
Let a(s) be a periodic function with period β, such that: Then, the following identity holds true, for all t ≤ 0:  (i) The function s → g β,η (−i s) satisfies the properties (4.28), recall Definition 3.5.
(ii) Notice that the function defined in (4.28) extends to a function a(z) on the lower-half complex plane, that is analytic for Im z < 0 and continuous for Im z ≤ 0.
The proof will be broken in a few intermediate steps.
In what follows, it will be convenient to use the following notations. We define inductively: Moreover, we set: Also, we shall introduce the n-dimensional symplex of side β as: ∆ n β := (s 1 , . . . , s n ) ∈ R n : β > s 1 > · · · > s n > 0 .   Proof. To start, let us prove the j = 0 case, which reads: (4.34) Let T > 0. By the KMS identity, Eq. (2.8), and using that B commutes with the number operator: (4.35) By assumption (4.28), we use the trivial but crucial fact a(i r ) = a(i r +β) to write:  Now, consider the function, for z ∈ C: For finite L and finite β, this function is analytic on Re z < 0, and it is continuous on Re z ≤ 0. In fact, the function z → 〈τ z (B )C 〉 β,µ,L is entire for finite L, β, while a(i z) is analytic for Re z < 0 and continuous for Re z ≤ 0, recall Definition (4.28) and Remark 4.3.
For ε > 0 small enough, let Γ be the complex path for (Re z, Im z): where every arrow corresponds to an oriented straight line in the complex plane. By Cauchy theorem, We start by writing: We claim that the last term vanishes as T → ∞. In fact: where we used the unitarity of the real-time dynamics. This estimate shows that the T → ∞ limit of the second term in the right-hand side of (4.41) vanishes, for β and L finite. Hence: which proves Eq. (4.34). Let us now discuss the j > 0 case, Eq. (4.33). By the KMS identity and a(i r ) = a(i r + β), we get: (4.44) We further rewrite this expression as, recalling the definition of B j (·), Eq. (4.31): (4.45) Let us now introduce the change of variables, for 1 ≤ k < j : Notice that β > s 1 > s 2 > · · · > s j > 0, we also have β > s 1 > s 2 > · · · > s j > 0, that is (s 1 , . . . , s n ) ∈ ∆ j β . In terms of these variables, for 2 ≤ k ≤ j : We then rewrite the term L T 1 in (4.45) as: (4.48) Let us now introduce the function: The function f (β,s 1 ,...,s j ) (z) is analytic for Re z < 0 and continuous on Re z ≤ 0. We have:  ,s 1 ,...,s j ) (r ). (4.51) As for the j = 0 case, we will use a complex deformation argument to rewrite L T 1 − L T 2 in a convenient way. To this end, let us now define the complex path for (Re z, Im z), for ε > 0 small enough: (4.52) By continuity of f (β,s 1 ,...,s j ) (z): where the second identity follows from Cauchy theorem and the last from the continuity of the integrand. The last term in the right-hand side of (4.53) vanishes as T → ∞. This is implied by the following estimate, recall Eq. (4.49): Consider now the first term in the right-hand side of (4.53). The integrand has the form, for a function g entire in all its arguments, recall (4.49): Let us introduce the change of variables, for 2 ≤ k ≤ j : We notice that β > s 1 > . . . > s k > . . . > s j +1 > 0. Thus, the second term in the r.h.s. of (4.53) can be written as the integral over the symplex ∆ β j +1 :  All in all, from (4.45), (4.53), (4.54), (4.57), relabelling the s variables as s variables: which concludes the proof of the proposition.
Next, we use Proposition 4.4 to rewrite the coefficients appearing in the Duhamel expansion in terms of imaginary-time correlations.  where the last equality follows from Proposition 4.4 for j = 0, applied to the r n integration. Next, using again (4.30), we write:  To conclude, recall that by Eq. (4.31):

Proposition 4.5 (Multiple complex deformation). Under the same assumptions of Lemma 4.2 the following identity holds:
where in the last step we used that B commutes with the number operator, which implies τ −i s (B ) = γ s (B ). Plugging (4.65) into the right-hand side of (4.64), and using the invariance of the Gibbs state under time-evolution, the final claim (4.59) follows.
Next, we rewrite the imaginary-time expressions appearing after the Wick rotation as connected correlation functions. We recall the following relation between the expectation value of a product of operators, and the truncated expectations: 〈O i 1 · · ·O i n 〉 = where the sum is over partitions P of {1, . . . , n + 1}, and J are the elements of the partition. In particular, J n+1 is the element of the partition that contains n + 1.
The right-hand side of (4.71) can be rewritten as: which we rewrite as, using again (4.66): The next proposition allows to rewrite the right-hand side of (4.59) in terms of truncated correlation functions, in Euclidean time.

Proposition 4.7 (Reduction to connected Euclidean correlations). Under the same assumptions of Lemma 4.2 the following identity holds:
a(i t + s j ) 〈γ s 1 (B ); γ s 2 (B ); · · · ; γ s n (B ); A〉 β,µ,L .    We claim that, for all k > 0: (4.79), we proceed as follows. First, we write: · 〈γ s π(1) (B )γ s π(2) (B ) · · · γ s π(k) (B )〉 (4.81) We claim that G is β-periodic in all its arguments: In particular, the function G extends to a periodic function on R k , with period β in all variables. With a slight abuse of notation, let us denote by G such periodic extension. Notice that the function s → a(i t + s) is also periodic with period β (recall definition (4.28)), which means that the whole integrand in the righthand side of (4.80) can be extended to a β-periodic function in all its arguments. Furthermore, we claim that, for all σ ∈ R: a(i t + s i ) G(s 1 , s 2 , . . . , s k ), (4.84) where S 1 β = R/βZ. We rewrite this expression as: where in the last step we used that e , as a function of s j , j = 2, . . . , k, is a function on (S 1 β ) k−1 . Then, the claim (4.79) follows from (recall that we can assume ω j ≥ 2π β , sinceã(0) = 0): This concludes the proof of Proposition 4.7. Proposition 4.7,Eq. (4.74) can also be rephrased as:

Remark 4.8. By the same arguments used in the proof of
a(i t + s j ) 〈Tγ s 1 (B ); γ s 2 (B ); · · · γ s n (B ); A〉 β,µ,L (4.87) where T denotes the time-ordering, as defined in Eq. (2.9). This is initially defined for operators whose imaginary-time arguments are in [0, β

Cumulant expansion for the instantaneous Gibbs state
In this section we shall review the well-known cumulant expansion for the Gibbs state of the Hamiltonian H (ηt ), that is:  For finite L and finite β, the series is norm convergent, thanks to the boundedness of the fermionic operators. Thus, the expectation value of a local operator on the Gibbs state of H (ηt ) can be written as: [0,β) n d t 〈Tγ t 1 (P ) · · · γ t n (P )〉 β,µ,L (4.96) which is analytic in ε for |ε| small enough. We would like to show that analyticity in ε extends to a ball whose radius is bounded uniformly in L. To this end, Eq.
Then, it is not difficult to see that the right-hand side can be written as a sum over time-ordered cumulants, defined as in Eq. (2.13). We have: (−εg (ηt )) n n! d s 〈Tγ s 1 (P ); · · · ; γ s n (P ); O X 〉. Under the assumption (3.2), the series converges for |ε| small enough, uniformly in L. By using Lemma 4.2, we will show that, for η small enough, the Duhamel series of the auxiliary dynamics is term-by-term close to the cumulant expansion of the instantaneous Gibbs state, Eq. (4.98).

Conclusion: proof of Theorem 3.7
We are now ready to prove our main result, Theorem 3.7.
To conclude the section, we discuss the proof of Corollary 3.9. (4.121) Next, we estimate the error introduced by replacing g β,η (s) with g (ηs). We have: where in the last step we used the Lieb-Robinson bound, as in (4.19). Finally, proceeding as in (4.25)-(4.27), we get: This concludes the proof of (3.21) and of the corollary.

A On the switch functions
Here we will discuss functions g (t ) that satisfy Assumption 3.3. This assumption holds true for the standard switch function g (t ) = e at with a > 0, where h(ξ) = δ(ξ− a), and more generally for the finite linear combinations of such functions. More generally, it is a natural question to understand under which conditions a function g (t ) can be represented as in (3.4), for a function h that satisfies the desired properties. Here we will give sufficient conditions for this to hold.
Let δ > 0 and let g (z) be analytic for Re z < δ. Suppose that, for all 0 ≤ k ≤ d + 2 and for all x < δ: Furthermore, suppose that for all 0 ≤ k ≤ d + 2 and for all y ∈ R: Examples of functions satisfying these assumptions are: with n ≥ d + 4 and a > δ. Let us check that functions satisfying (A.1), (A.2) verify Assumption 3.3. Let γ be the straight complex path on the imaginary axis γ : We have: which is well defined thanks to (A.5). We have: In the first identity we applied Fubini's theorem, and in the last identity we applied Cauchy integral formula, using that g (1 + ξ k )h 1 ≤ C k , (A. 8) which shows that the second assumption in Eq. (3.5) holds true. Let us now consider the first assumption in (3.5). We observe that since the function g (z) has no poles for Re z < δ: Eq. (A.9) follows from Cauchy integral formula, combined with the integrability properties (A.1), (A.2). Furthermore, for all 1 ≤ k ≤ d + 2: where we used that z k g (z) is integrable on the path γ, by (A.1). By the same assumption, ∂ k ξ h(ξ) is bounded uniformly in ξ for 1 ≤ k ≤ d + 2.
The function z k g (z) has the same analyticity properties as g , in particular it has no poles for Re z < δ. Thus, by Cauchy integral formula, thanks to the integrability assumptions (A.1), (A.2) we have: = 0.

B Properties of Euclidean correlations
In this appendix we shall give the proofs of the properties of Euclidean correlation functions that have been used in the proof of Theorem 3.7. These properties are well-known, and we collect the proofs here for completeness. Recall the notation for the n-dimensional symplex of side β: where in the last step we set s i = β, and used the fact that all the permutations contributing to the last sum are such that π(1) = i . The right-hand side of Eq.

C Decay of correlations for interacting models
In this appendix we shall discuss the validity of Assumption 3.1 for fermionic lattice models. The arguments presented in this section are well-known, and we reproduce them for completeness. For definiteness, we shall consider the following class of many-body Hamiltonians: where H (x; y) and v(x; y) are finite-range. The discussion that follows actually applies essentially unchanged to a larger class of Hamiltonians, obtained replacing the quartic interaction in (C.1) by finite-range interactions of arbitrary even degree in the fermionic creation and annihilation operators. Let 〈·〉 β,µ,L be the Gibbs state of the system: For |λ| small, the Gibbs state can be expanded in power series around the noninteracting one, following Section 4.4. The main advantage in doing this is that the non-interacting Gibbs state is quasi-free, which means that all Euclidean correlations can be computed starting from the Euclidean two-point function using the Wick rule. Let t , t ∈ [0, β), t = t . Let us denote by 〈·〉 0 β,µ,L the noninteracting Gibbs state (λ = 0). We define the non-interacting two-point function as: At equal times, the two-point function is defined by normal ordering: Eq. (C.3) is extended to an antiperiodic function over all t , t ∈ R with antiperiod β. The next proposition gives the explicit expression of the two-point function.
Proposition C.1 (Non-interacting two-point function). Let 0 ≤ t , t < β. Then: Proof. Eq. (C.5) can be proved by direct computation of the trace involved in the definition of the non-interacting Gibbs state, representing the Fock space in the basis of Slater determinants associated with the eigenstates of the Hamiltonian H . The computation can be found in standard textbooks in condensed matter physics, see e.g. [20,42].
Next, we collect useful decay estimates for the two-point function.
Proposition C.2 (Bounds for the two-point function). There exist C β , c β > 0 such that the following bound holds true, for all L > 0: Moreover, suppose that µ ∉ σ(H ) and dist(µ, σ(H )) ≥ δ with δ > 0 uniformly in L. Then, the constants C β , c β can be chosen uniformly in β: Proof. Let us start by proving (C.6). Suppose that t > t , the other case can be studied in the same way. The starting point is the following formula for the twopoint function: where the first identity follows from Proposition C.1, and in the second equality the complex path C is a rectangle that encircles the spectrum of H , and that crosses the imaginary axis at Im z = ±π/(2β). Thus, the path does not enclose any of the poles of the Fermi-Dirac function, which are given by z − µ = i 2π β (n + 1 2 ) with n ∈ Z, and it stays away from the spectrum of H . By construction, the only singularities encircled by the path C correspond to the eigenvalues of H . Since H is bounded uniformly in L, the length of the complex path is also bounded uniformly in L; in particular, we can choose C such that for all z ∈ C we have dist(z, σ(H )) ≥ π/(2β). The estimate (C.6) easily follows from the Combes-Thomas estimate for the Green's function, see e.g. [3].
Let us now suppose that dist(µ, σ(H )) > δ, and let us prove the estimate (C.7). We can use a complex representation (C.8), where now the path C splits into the disjoint union of two non-intersecting paths, C − and C + , that encircle the spectrum of H on the left or on the right of µ, respectively. The paths can be chosen so that the distance from any point z ∈ C to the poles is bounded below by a constant proportional to the spectral gap, uniformly in β. We write: (C.9) For z ∈ C + : (C.10) Instead, for z ∈ C − : ) . (C.11) All together, for z ∈ C : The exponential decay in space follows from a Combes-Thomas estimate for the Green's function in the complex integral, using that now the distance from z to the spectrum of H is bounded uniformly in β. This concludes the proof.

(C.22)
Since the number of sets X i with bounded diameter and containing a given point z i is bounded, we can estimate the X i sum as: where z n+1 = x ∈ X . A crude bound for the last sum is C n d 1 · · · d n , uniformly in the volume of the system (remember that the point z n+1 is fixed). Plugging this estimate in (C.23), we get: Thus, we obtain the following bound for the cumulant: (C.25) The integral in the right-hand side is bounded by a power of β, if no spectral gap is present. Instead, if µ lies in a spectral gap we have, recalling Eq. (C.21): with the understanding that t n+1 = 0. Recall that |t | β = i |t i | β . In the permutation π(1), π(2), . . . , π(n + 1), suppose that i = π(k) and n + 1 = π( j ) and that j > k (the other case is analogous). Then, we write: Plugging this into (C.26), we find: (C.28) Using this estimate in Eq. (C.25) we finally find: where the factorial comes from the sum over permutations. This concludes the check of Assumption 3.1 for non-interacting fermions.

C.2 Check of Assumption 3.1 for interacting fermions
Let us now discuss the case of interacting Fermi systems, with Hamiltonian H = H 0 + λV given by Eq. (C.1) with λ = 0. Let us denote by 〈·〉 λ β,µ,L the Gibbs state of the Hamiltonian H . Following Section 4.4, this Gibbs state can be written as a series in cumulants of the many-body interaction over the non-interacting Gibbs state: More generally, the interacting Euclidean correlation functions can be written in terms of the non-interacting ones, as: where the integral is over [0, β) m . All cumulants can be computed in terms of connected Feynman diagrams, using Wick's rule, and the bounds of Proposition C.2 allow to prove estimates for the Feynman diagrams that are uniform in L, and also in β if µ is in a spectral gap of H . With respect to the case discussed before, the main problem now is that the observables at the argument of the time-ordering might be quartic in the fermionic operators: this ultimately implies that the number of Feynman diagrams contributing to the order n grows as (n!) 2 , which beats the 1/n! factorial in Eq. (C.31). For fermionic models, this combinatorial problem is only apparent, as it can be solved keeping track of the minus signs arising from the anticommutation of the fermionic operators. The mathematical tool that allows to prove a bound for the cumulants that grows only as n!, and that is uniform in the size of the system, is the Brydges-Battle-Federbush-Kennedy (BBFK) formula [10,16,17,18], for the connected expectations, or cumulants, of a fermionic theory. See [21] for a review of recent applications to transport problems in condensed matter systems. Let us review its application to the problem at hand.
Let A(P i ) be a short-hand notation for a monomial in the creation and annihilation operators, A(P i ) = γ t i (a * x i ,1 ) · · · γ t i (a * x i ,k )γ t i (a y i ,k ) · · · γ t i (a y i ,1 ). (C.32) P i has to be understood as a set of points, labelled by a sign ε i = ±, which denotes creation operators (ε = +) or annihilation operators (ε = −), and by spacetime coordinates (x i , t i ) if ε i = + or (y i , t i ) if ε i = −. Without loss of generality, we can suppose that x i ,k = x i , and y i ,k = y i , for k = , since otherwise A(P i ) = 0 by Pauli principle. Monomials of this type appear when writing explicitly the operators V , O (i ) X i at the argument of the cumulant in Eq. (C.31) in terms of the fermionic operators. The BBFK formula provides a very useful identity for 〈TA(P 1 ); A(P 2 ); · · · ; A(P n )〉 0 β,µ,L (C. 33) in terms of the two-point function, Eq. (C.5). One has, if t i = t j for i = j : 〈TA(P 1 ); A(P 2 ); · · · ; A(P n )〉 0 Let us explain the various objects and symbols entering this formula. We view the monomial P i as being represented by a cluster of points, labelled by z, t variables, with a line attached: the line is incoming if the associated fermionic operator is γ t (a z ), or outgoing if the associated fermionic operator is γ t (a * z ). Each line is labelled by a label f , and we shall think of P i as being the collection of such labels. The T -sum in Eq. (C.34) is a sum over anchored trees between the cluster of points P 1 , . . . , P n , where the edges of the trees are associated with the contractions = ( f , f ) of incoming and outgoing lines. An anchored tree between the clusters associated with P 1 , . . . , P n becomes a tree between n points if one collapses the clusters into points. With each edge of the tree we associate a propagator g , g ≡ g 2 (x( ), t ( ); t ( ), y( )), (C. 35) where x( ), t ( ) are the space-time labels for the contracted outgoing line f , while y( ), t ( ) are the space-time labels for the contracted incoming line f . Notice that, in general, the lines forming the trees are only a subset of all the possible contractions that can be made between the lines associated with the clusters P 1 , . . . , P n (which are the contractions forming the connected Feynman diagrams). Informally, given a tree T the sum over the remaining contractions that exhaust all lines is taken into account by the integral in (C.34). There, s denotes variables (s i j ) i , j =1,...,n , and d µ T (s) is a T -dependent probability measure supported on a set of s i j ∈ [0, 1] such that s i j can be written as a scalar product (u i , u j ) for a family of vectors (u i ) with u i ∈ R n of unit norm. Finally, s i ( f ),i ( f ) g ( f , f ) is a matrix, labelled by the lines that are not part of T : f , f ∈ (∪ i P i ) \ P T , where P T takes into account all the f , f labels that are involved in the tree T . The notation i ( f ) indicates the label of the cluster P i to which the label f belongs to. Finally, α T is a suitable function of the tree T , and it takes the values ±1. Its value will not be important in the following.
Eq. (C.34) allows us to obtain the estimate: |〈TA(P 1 ); A(P 2 ); · · · ; A(P n )〉 0 β,µ,L | the product over the propagators associated with the branches of the tree introduces a decay factor as function of the space-time distance of the various P i 's.
To make good use of Eq. (C.36), we need a bound for the determinant. One possibility could be to express the determinant in terms of the matrix entries via Leibniz formula, but this would ultimately produce the same combinatorial growth observed in the naive Feynman graph expansion, which is useless for the purpose of proving convergence of the series in (C.31). A better estimate is obtained if g ( f , f ) are the entries of a Gram matrix, that is if g ( f , f ) = (a f , b f ) for a f , b f vectors of finite norm in a Hilbert space. In fact, in this case we could apply the Gram-Hadamard inequality, to obtain: where the product runs over all the labels in (∪ i P i ) \ P T . This bound grows as a power of the dimension of the matrix, instead of factorially, and can be ultimately used to prove convergence of the series in (C.31). The problem in our case, and in all applications to lattice fermionic models, is that the propagator g ( f , f ) cannot be expressed in a Gram form. This problem could be solved via an ultraviolet multi-scale decomposition of the Matsubara frequencies associated with imaginary times, as reviewed for instance in [24]. This analysis can be viewed as a warm-up for the multiscale analysis needed in order to tackle the infrared problem of interacting, gapless systems. Alternatively, a relatively simple way out to this problem to observe that g ( f , f ) can be expressed as the linear combination of the entries of Gram matrices [43,19]. We Notice that, by assumption on the model, the sums involve sets with bounded diameter, diam(Y i ) ≤ R. In order to control the time integrals and the sums over lattice subsets, we use the estimate (C.41); equal times give zero contribution to the integral, hence we can assume that all times are different. Let us consider: (C. 44) the T -sum is over the anchored trees connecting the cluster of points associated with Y 1 , . (1 + |s | β ).
The constant C n+m takes into account the sum over Y i z i , using that their diameter is bounded. For a given tree, the sum over the space-time coordinates of the points is performed via a standard pruning argument. One starts from the leaves of the trees, which are defined as the points attached to the rest of the tree by just one branch. If the leaf is labelled by x one does nothing, otherwise we integrate over the corresponding space-time variable and get a factor: which does not depend on t by time-translation invariance of the estimate for the two-point function, Proposition C.2. After having integrated out the leaves one deletes them, thus obtaining a new (smaller) tree. We then integrate over the new leaves, and repeat the process until all integrations are exhausted. By doing so, we obtain the bound: where Γ n+m+1 is the number of trees with n +m +1 vertices. By Cayley's formula, it is well-known that: Γ n+m+1 = (n + m + 1) n+m−1 .