Finite temperature and quench dynamics in the Transverse Field Ising Model from form factor expansions

We consider the problems of calculating the dynamical order parameter two-point function at finite temperatures and the one-point function after a quantum quench in the transverse field Ising chain. Both of these can be expressed in terms of form factor sums in the basis of physical excitations of the model. We develop a general framework for carrying out these sums based on a decomposition of form factors into partial fractions, which leads to a factorization of the multiple sums and permits them to be evaluated asymptotically. This naturally leads to systematic low density expansions. At late times these expansions can be summed to all orders by means of a determinant representation. Our method has a natural generalization to semi-local operators in interacting integrable models.


Introduction
As a consequence of the existence of extensive numbers of conservation laws with local densities the dynamical properties of quantum integrable models at finite energy densities are both rich and unusual. The two main settings of interest are finite temperature equilibrium response and time evolution after quantum quenches. In the first setting the aim is to determine two-point functions of the form where Z(β) = Tr(e −βH ) and A(x, t) is a Heisenberg picture operator, while in the quench setting one is interested in equal time expectation values Ψ|O(x, t)|Ψ , where |Ψ is an initial state that is a linear superposition of an exponentially (in system size) large number of energy eigenstates.

Finite temperature dynamics
Early work on determining (1) focussed on the spin-1/2 XX chain, which can be mapped to a non-interacting model of free fermions [1]. The long time and distance asymptotics of two-point functions was obtained from the solution of a Riemann-Hilbert problem in [2] arising from a Fredholm determinant representation [3,4]. A semiclassical approach to the low temperature regime in interacting integrable models was pioneered in [5] and has proved very useful [6][7][8] due to its relative simplicity. It is however limited in that it applies only to very low temperatures and cannot be easily extended. Perhaps the most direct approach to evaluating (1) or (2) is by introducing spectral representations, e.g.
where |n are eigenstates of energy E n and momentum P n . Early investigations of (3) focussed on integrable quantum field theories in the infinite volume [9][10][11][12][13][14], where the spectral representations need to be regularized. This problem was solved in [15][16][17] and a systematic low temperature expansion of dynamical two point functions in Fourier space was obtained [17][18][19][20]. For some correlators this expansion exhibits divergences close to the zero temperature mass shell and needs to be summed to all orders. Carrying out this summation exactly is an open problem. A similar approach was formulated for the case of the Ising field theory in [34] and used to obtain the late-time asymptotic behaviour of the order parameter two-point function [35].
In order to go beyond the low temperature regime in interacting integrable models it is useful to work in the micro-canonical ensemble and employ typicality ideas. This provides a more efficient spectral representation of the form where E β is a typical energy eigenstate at the energy density corresponding to inverse temperature β [21]. The representation (4) can be analyzed numerically for finite systems [21]. Moreover, in particular limiting cases it appears to be very efficient in that only a small number of states need to be summed over [22]. Very recently an axiomatic approach to formulating form factors between states at finite energy densities in the infinite volume limit was proposed [23] and used to formulate a spectral representation. Using this representation to obtain explicit results for dynamical two-point functions remains an open problem. An alternative approach to finite temperature dynamics is based on the Quantum Transfer Matrix approach [27,28]. The latter is highly efficient for determining static properties [29][30][31][32] and can be extended to dynamical correlation functions [33]. Very recently this method has been successfully applied to the XX model [24][25][26] and state-of-the-art results have been obtained. The generalization to determine dynamical two-point functions in interacting integrable models is an open problem.
The late time asymptotics of certain finite temperature two-point functions can also be accessed by applying generalized hydrodynamics [36,37] to the linear response regime, see [38,39].

Quench dynamics
Early work on quench dynamics again focussed on models that can be analyzed by means of free fermion techniques [40][41][42][43][44][45][46][47]. Notably, in [44,45] exact results for the late time behaviour of one and two point functions of the order parameter in the transverse field Ising model (TFIM) after quantum quenches were obtained. One way of going beyond free theories is to employ a spectral representation Ψ|O(x, t)|Ψ = n,m Ψ|n m|Ψ n|O(0, 0)|m e i(En−Em)t−i(Pn−Pm)x .
This was used to obtain the late time behaviour for small quenches in the TFIM [44,45,[48][49][50][51] and the sine-Gordon model [52]. The small quench regime is also accessible by semiclassical methods [53][54][55][56][57], which have the advantage of being significantly simpler to implement. A much more efficient spectral representation is provided by the Quench Action Approach [58]. For translationally invariant initial states this allows one to express expectation values of local operators after a quantum quench from an initial state |Ψ as where |Φ s is a representative state fixed by two requirements: first, it is a simultaneous eigenstate of the Hamiltonian and of the (quasi)local conservation laws I (n) of the theory under consideration, and, second, it correctly reproduces the expectation values Expression (7) affords a more efficient spectral representation involving only a single sum over energy eigenstates as The behaviour in the steady state reached in the limit t → ∞ is given by the expectation value in the representative state and this has been analyzed in a number of cases [62][63][64][65][66][67]. The time dependence is significantly more difficult to obtain. So far results are restricted to a particular one-point function for small quenches in the sine-Gordon model [52] and density correlations at late times after a quench in the repulsive Lieb-Linger model [68].

Local vs semi-local operators
Locality properties of the operator of interest have important implications in both finite temperature and quench contexts. For quantum quenches this was emphasized in [42,43] and clarified through explicit calculations in Refs [44,45,52,58]. A precise definition of the mutual locality index ω(A, B) of two operators exists in the context of relativistic integrable quantum field theory, see e.g. [69]; specifically, the product of operators A(x, τ )B(0, 0) as a function of (x, τ ) has the property where A C denotes the analytic continuation along a counter-clockwise contour C around zero. Let us for simplicity consider the case of a diagonal scattering theory with only a single "elementary" particle excitation created by the field Ψ(x). A convenient basis of energy eigenstates is given in terms of scattering states of elementary excitations where θ j are rapidity variables related to the energy and momentum of a single-particle excitation by ǫ(θ) = M cosh(θ) , p(θ) = M v sinh(θ). Spectral representations of correlation functions (in the infinite volume) involve form factors like As we will see below the case M = N is of particular interest. Local operators have vanishing mutual locality index with Ψ(x). As a consequence of kinematic poles [70] the form factors become singular when rapidities in the set {θ j } approach those in {θ ′ j }. In the case N = M the structure of singularities is [16] In contrast, for semi-local operators B with ω(B, Ψ) = 1 one has instead This shows that form factors of such semi-local operators are much more singular than the ones for local operators. Form factors in integrable lattice models have analogous structures of singularities.

One and two-point functions of semi-local operators
The nature of singularities for semi-local operators (13) has been exploited previously to obtain results for 1-point functions after small quantum quenches [52,58]. The aim of this work is to extend this approach to general quantum quenches as well as to dynamical twopoint functions at finite temperatures. We focus on the case of the order parameter in the TFIM because the form factors are particularly simple in this case. This allows us to exhibit in considerable detail which states in the respective spectral representations contribute to the late time asymptotics of one and two-point functions. These considerations can be generalized to interacting integrable models, as will be shown in a following publication.

Transverse Field Ising Model
The TFIM Hamiltonian on a ring with L sites reads where σ α j acts like the corresponding Pauli matrix at sites j and like the identity elsewhere. We assume J, h > 0 and consider periodic boundary conditions. We refer the reader to appendix A of [44] for details about the diagonalization of this Hamiltonian. We simply recall here that it can be expressed in terms of free fermions α k as and that the Hilbert space is divided into a Neveu Schwartz (NS) sector with states of the form and a Ramond (R) sector where the energy E ({p i }) and momentum P ({p i }) of such states are given by with an identical relation for the Ramond sector. We will be interested in two problems involving the summation of form factors of the order parameter over the full Hilbert space, which are given by [72][73][74][75] Here ξ = |1 − h 2 | 1/4 , m is even (odd) for h < 1 (h > 1), and The terms ξ L and η k do not depend on the momenta (except k) and for large L approach one with exponential accuracy ξ L ≈ 1 , e η k ≈ 1 ; thus, they will be set to 1 in the following.

Quenches in the quench action framework
We consider the following quantum quench setup [44,45]: at time t = 0 we prepare the system in the ground state of the TFIM (14) at a magnetic field h 0 < 1 At times t > 0 we evolve the system with Hamiltonian H(h) with h 0 = h < 1. As |Ψ is not an eigenstate of H(h) this results in interesting dynamics. The order parameter one-point function at time t > 0 is given by (6), where the representative state |Φ s is characterized by the root density [58] ρ (k) = 1 − cos ∆ k 4π , For later convenience we define the density of particles in the representative state In a large finite volume L we may choose [58] where the momenta q i are distributed according to the root density ρ(k). The time-evolved initial state is given by with K (p) = tan(∆ p /2). Eqn (6) thus provides the following representation for the order parameter one-point function that is a sum of form factors over states that are expressed in terms of pairs of momenta.

Dynamical correlation functions at finite temperature
In the TFIM, the density of momenta q of the representative state |E β in (4) is It is very useful in the following to define the corresponding density Inserting a resolution of the identity between the two spin operators in (4) leads to the following spectral representation (30) The terms in the sum depend on the regime of the TFIM: in the ordered phase, h < 1, M has the same (even/odd) parity as N , whereas in the disordered phase, h > 1, it has opposite parity. Moreover, in contrast to the quench case, (30) involves squares of form factors, and the intermediate states do not have a structure where momenta only appear in pairs {−p i , p i }.

Systematic approach to form factor expansions for semilocal operators
In this section we present a general framework for carrying out the form factor sums (27) and (30) analytically at late times Jt ≫ 1. It is based on decomposing the form factors (19) into partial fractions so that the sums over the p's decouple and can be evaluated exactly.
The key observation is then that an oscillatory sum with a pole of order d like n e int (n+1/2) d grows as t d−1 , so that the leading poles give the leading time behaviour, and the terms in the partial fraction decomposition can be organized according to the total number of poles. This naturally leads to an expansion in the number of particles per unit site -N/L -in the representative state, which has already been proven very efficient for simpler quantities such as the free energy in Bethe ansatz solvable interacting models [59,60].
In sections 3.1 and 3.2 we consider the application of this framework in the context of the finite temperature case (30), due to the more canonical sum over form factors that it involves. Since the form factors differ for h < 1, h > 1 and h = 1, we will treat these cases separately.
We will be interested in large time or space asymptotics of correlation functions, generically defined as requiring the phase it(E ({q}) − E ({p})) + il(P ({p}) − P ({q})) in (30) to be large. This is in particular the case of the large time and distance asymptotics at fixed on which we will focus. However, the static correlations case t = 0 and large l is also covered by our calculations; we refer the reader to section 3.5.3 for details on this case. For later convenience we introduce the following notations where v max is the maximal group velocity of the elementary fermion excitations in the TFIM. According to whether there exists an x 0 such that ε ′ (x 0 ) = 0 ('time-like region', tv max > ℓ for t, ℓ ≥ 0) or not ('space-like region', tv max < ℓ for t, ℓ ≥ 0), the next-to-leading terms in our expansions differ. We will in the following compute these terms in the space-like region. Their calculation in the time-like region is more involved and will be reported elsewhere.
In section 3.4 we briefly present the application of our framework to the dynamics after quantum quenches.

Case h < 1: identical number of particles
In this subsection we treat the case h < 1, for which the sum (30) includes intermediate states with the same number N of particles as the representative state. We will show in section 3.2 that the contributions of intermediate states with different particle numbers M = N is always subleading in time. We exploit this result right away and focus on states with M = N in this subsection.

Partial fraction decomposition of form factors
We recall that the partial fraction decomposition of a ratio of two polynomials with P 0 (X) a polynomial of degree deg(P ) − i a i and B i,ν independent of X, given by The squared form factor appearing in (30) can be written as The ε(p 1 ) factors are not polynomial in p 1 but are nevertheless bounded and without zeros. Seen as a function of p 1 , it can thus be written i A i , B i independent of p 1 and C a bounded function of p 1 . Repeating the operation for the other momenta, one can write where the second sum is over a complete set of functions  (36) is that each p appears at most once (however, the q's may appear several times).

Carrying the sum over the momenta p i
Let us briefly anticipate the method that we will use to carry out the the sum over p in (30). If there is a ν i = 0 then the sum over p i is an oscillatory Riemann sum of a bounded function, hence it decays to zero with time. Thus the leading behaviour is obtained for ν i > 0 for all i. Then (and only in this case) the coefficients A ({q} , {p}, {ν} , f ν ) are independent of the p's and the maps will no longer depend on {ν j }. These coefficients will be denoted by A ({q} , {ν} , f ) and are obtained as 1 In F V U the sets U and V can contain arbitrary momenta and are not meant to be each in the sectors R and N S (this freedom will be indeed useful in section 3.1.6 below). For this reason we impose in the denominator that the momenta are different u = v, which is automatically satisfied if they are in different sectors, but not otherwise.
This follows from the partial fraction decomposition, with a factor 2 because of the 1/2 inside the sin.
From here on we will only consider terms in the partial fraction decomposition such that ν i > 0 for all i, so that (37) applies. We denote the corresponding contribution to the spectral representation (30) where the third sum is over a complete set of functions f : {1, ..., N } → {1, ..., N }.
In this form one can perform the sums over p j using the following relations proven in Appendix A with q ∈ NS, to obtain Eqns (40) are valid only when ε ′ (q) = 0. If there is a point where ε ′ (q) = 0, ie if we are in the time-like region, corrections in time to (40) have to be taken into account, and they are expected to modify significantly the subleading corrections in the correlation function. We leave this matter of discussion for future work.

Constraints on the functions f
The set of functions f over which we need to sum is actually quite constrained. First, since F We also need f ..,N , but as there is only one derivative with respect to p i and none with respect to p j this factor will make the coefficient More generally, if f takes k times the same value at points with ν i = 1, then all the k(k−1)/2 terms sin 2 p i −p j 2 contribute to a zero of order k(k−1); since the number of derivatives is equal to k, we must have k = 2.
These arguments show that the sum over f can be replaced by a sum over three disjoint subsets I 0 , I 1 , I 2 ⊂ {1, . . . , N }, where I k is the set of points with ν = 1 attained k times by f . The remaining points {1, ..., N } − (I 0 ∪ I 1 ∪ I 2 ) all have ν = 2. There is a combinatorial factor N ! 2 |I 2 | corresponding to the number of such functions with this precise ouput. It follows that A ({q} , {ν} , f ) depends only on the sets I 0 , I 1 , I 2 and we have The expression for the coefficients A(I 0 , I 1 , I 2 ) can be simplified as follows. We observe that, whenever we have ν i = 2 in (37), the various factors depending on p i and q f (i) precisely compensate one another. Hence we can work with a reduced form factor involving only momenta in I 0 , I 1 , I 2 for any function f such that The set I 2 does appear twice by construction, and I 0 does not appear at all. The decomposition of {1, ..., n + 2m} into {1, ..., n}, {n + 1, ..., n + m}, {n + m + 1, ..., n + 2m} is arbitrary, and it needs only to involve one set with n elements and two sets with m elements.
The sum of all the terms in (42) with |I 1 | = n, |I 0 | = |I 2 | = m will be denoted by S n,2m . We can factorize S 0,0 and, using the explicit expression for χ 1 (q) and χ 2 (q), write If n and m stay finite in the limit L → ∞ we have i∈I 0,1, Since the momenta selected by the sets I 0,1,2 are drawn from the momenta {q j |j = 1, . . . , N } of the representative state with density ρ, the term S n,2m is of order (N/L) n+2m times S 0,0 . Hence this expansion naturally leads to an expansion in N/L.
where, for a function f (x), we used This term is the correlation function at order 1 in ρ uniformly in large t.
This term leads to and Although they are individually both divergent in L in the scaling limit L → ∞, their sum is not divergent and is, see appendix A, with the following value in the space-like regime where sgn ( Contrarily to the previous case, this order ρ 2 of the correlation function at fixed large t cannot be uniform in t, since it diverges for t → ∞. 3.1.6 Recursive structure of A(I 0 , I 1 , I 2 ) In the general case the amplitudes A(I 0 , I 1 , I 2 ) are obtained from (43), but the sums over momenta associated with the index sets I 0 , I 2 in (44) cannot be carried out as simply as in the cases treated above. In fact, as we noted earlier, the derivatives corresponding to I 2 , I 0 in (43) have to be applied on the double zero sin 2 p i −p j 2 in the numerator to give a non-vanishing result, so that one actually has The form factor F can itself be decomposed into partial fractions. One obtains We observe that the sum over {q 2 } in (44) will play a similar role to the sum over p's in (30), with however the important difference that they are drawn from the original q's of the representative state and are not arbitrary momenta as is the case for the p's in (30).
3.1.7 Partial fraction in A(I 0 , I 1 , I 2 ) leading in density Relation (54) reveals a recursive structure in the calculation of A(I 0 , I 1 , I 2 ). However, we will not develop this recursion further here, but will rather focus on the leading partial fractions of (54) obtained with a set {ν} such that ν i = 2 for i ∈ I 2 and ν i = 1 for i ∈ I 1 , and functions f : I 1 ∪ I 2 → I 1 ∪ I 0 that map I 2 to I 0 in a one-to-one fashion and fulfil f (i) = i for i ∈ I 1 . Figure 5: Sketch of the two functions f after one step of recursion. In red are indicated the p's, in blue the q's that play the role of the p's after the first step of the recursion, that are those with index in I 1 or I 2 . The 'new' q's on the last row are those with index in I 1 or I 0 .
The coefficient then reads Let us select one i ∈ I 1 and introduce a reduced set I ′ 1 = I 1 − {i}. Performing the derivative with respect to p i gives (56) We observe that the first factor in the second line has to be differentiated precisely one more time for the result not to vanish (the fact that f is one-to-one on I 1 to I 1 is essential for this), so that with I ′′ 1 = I 1 − {i, j}. Applying the same reasoning to the remaining momenta yields 3.1.8 Result: correlation function at O(ρ 2 β ) uniformly in t at large t In analogy with our notations for the first level of the recursive structure we denote by S n,2m|0,0 all contributions to (44) that arise by specifying ν i = 2 for i ∈ I 2 in (54). We observe that in (53) the form factor vanishes if there are coinciding momenta in I 0 and in I 2 , or momenta that occur both in I 1 and I 0 or I 2 . Hence one only has to impose that momenta in I 1 are distinct among themselves, and that momenta in I 0 are distinct from those of I 2 . This gives The two sets of sums factorize. The first set can be calculated using the results of Appendix A (60) As for the terms in the sum over {q 1 }, they vanish for n odd, while for even n = 2p they are calculated in Appendix A where c is defined in (52). We note that these two calculations involve the same sums as in S 0,2 and S 2,0 above. It follows that the infinite volume limit of the sum of all S n,2m|0,0 reads n,m≥0 where S 0,0 has been computed in (45). We have thus obtained the order O(ρ 2 β ) contribution to the correlation function uniformly in t with In the time-like region, the exponent would be the same, but the constant would differ. This result should be compared to the semiclassical approach of Sachdev and Young [5,71], which gives As expected our result reduces to the semiclassical one in the limit βJ ≫ 1.

Case h > 1: different numbers of particles
We now turn to the case h > 1. Here the sum (30) involves only intermediate states with numbers of particles that are different from that of the representative state. Hence we must study form factor sums with M = N . We will compute the prefactors at order O ρ 1 β in the space-like region, and because of saddle points effects only at order O ρ 0 β in the time-like region.

General structure
As in the case M = N , the form factor | NS q 1 , ..., q N |σ x l |p 1 , ..., p M R | 2 can be decomposed into partial fractions The important difference to the case M = N is that we cannot always neglect the contributions with ν j = 0. Indeed, let us denote by k the number of ν j = 0. The corresponding contributions give rise to k oscillatory bounded integrals that will each decay with time. On the other hand each of the M − k sums over the other momenta p i will generate an oscillating factor e −itε(q f ν (j) ) according to (40), while N factors e itε(q k ) are already present. The resulting oscillatory sums may have singularities, but according to (40) summing these singularities does not consume any oscillatory factor, it only lowers the number of singularities. Hence in the end we will be left with |N − M + k| oscillatory bounded integrals to perform. In total, there are thus k + |N − M + k| of such integrals. Hence if M ≤ N the case k = 0 is still dominant at late times, but if M > N then all the cases 0 ≤ k ≤ M − N are a priori of the same order. In both cases these leading terms involve |M − N | oscillatory bounded integrals, so that we can conclude that the terms M = N are exponentially smaller than the case M = N in the space-like regime, and typically around 2 t −|M −N |/2 smaller in the time-like regime.
It follows that for h > 1 the dominant terms in (30) are obtained for M = N ± 1, which we now consider in turn.

Case
Using (40), the corresponding contribution to (30) is This term is of order ρ, so in the time-like region it vanishes at the order of our computation.
In the space-like region we have There, ε(x) is monotonous and e itε(x) is periodic (because the distance ℓ is an integer) so the first integral decays with time faster than any power-law and cannot be simplified further.  (40), the corresponding contribution to (30) is In the infinite volume limit we obtain the following result in the space-like regime where the prefactor is accurate to first order in the density ρ β . In the time-like region the prefactor is only accurate to O(ρ 0 β ), because saddle point effects arise at order O(ρ β ). A saddle point approximation gives in this regime where SP = {s, ε ′ (s) = 0} denotes the set of saddle points and ± the sign of tε ′ (s).
. Using (40), the corresponding contribution to (30) becomes Taking the infinite volume limit we obtain In the space-like region χ 1 (x) = i sgn (ε ′ (x)) and one has with the prefactor at order ρ 1 In time-like region, this term vanishes at order ρ 0 .

Result: correlation functions for h > 1 at leading order in ρ β
Putting everything together, in the space-like regime and at leading order in density we obtain the following result In the time-like regime we have instead where the sum is over the saddle points s of ε(x). This result should be compared to the semiclassical approach of Sachdev and Young [5,71], which gives As expected our result reduces to the semiclassical one in the limit βJ ≫ 1.

Case h = 1
In the case h = 1 the structure of the form factor (19) is modified compared to the case h = 1: since ε(0) = 0, there are additional poles. The nature of these poles is moreover different from those appearing in the partial fraction decomposition of the previous sections, since they involve pairs of momenta: for ǫ pp ′ to vanish we must have p = p ′ = 0. We note that at zero temperature the model becomes critical at h = 1 and correlation functions should exhibit power-law decays. These issues are beyond the scope of this work and will be examined elsewhere.

Quantum quench case
We now turn to the time evolution of the order parameter one-point function after a quench of the transverse field within the ordered phase. Our aim is to evaluate the spectral representation (27) obtained in the framework of the quench action approach.

Generalities
The sum over form factors appearing in the context of quantum quench dynamics (27) differs from the finite temperature case (30) notably because now in both states of the form factors the momenta come in pairs p i , −p i . One can then write Focusing again on M = N in (27), we have with

Differences to the finite temperature case
We apply a partial fraction decomposition to the second line of (79), seen as a ratio of polynomials in the cos p j . The procedure is the same as in section 3.1. More precisely, we define which we decompose into partial fractions with cos u taken to be the relevant variables. This gives There are however some noteworthy differences from section 3.1 brought by the presence of factors involving the function f (p) and the fact that there is still a p i dependence in the sum (79) outside the partial fraction decomposition. In place of (40) we now have as shown in appendix A. The analog of eqn (44) is given by with Here we used that f (p) (80) fulfils The possibility of differentiating f (p) also modifies the derivation of (58), which now takes the form Because of (86), however, only K = {} remains after the sum over q k .

Result
: O(ρ 2 Q ) uniformly in t at large t The final result of the above calculation for the time evolution of the order parameter onepoint function after a quench of the transverse field within the ordered phase is where This result holds at the second order in the density O(ρ 2 Q ). The reader may have noted that, since ε ′ (0) = ε ′ (π) = 0, the quantum quench dynamics is a 'time-like region' case, so the saddle point effects might modify this prefactor at order ρ 1 . It turns out, however, that ρ(0) = ρ ′ (0) = ρ(π) = ρ ′ (π) = 0, so the saddle point corrections are higher order in time and do not affect the prefactor.

Comments on the form factor summation
Having carried out form factor summations to obtain the leading late-time asymptotics in the low-density regime in both the finite temperature and the quantum quench contexts it is useful to take stock and stress some features that we expect to be of a general nature.

Which states govern the late time dynamics?
The leading late time behaviour follows from eqns (40) that enter the sum over p's of the partial fraction decomposition (36). These formulas are obtained by isolating the singularity and dropping the integral of an oscillatory bounded function, cf. Appendix A. The singular part involves momenta p j in (30) that are at distance O(L −1 ) of any of the q i . The aim of this section is to quantify more precisely the number of p j that have to be summed in order to recover the late time dynamics. In other words, we would like to know the smallest function η(L) such that with p n = q + 2π L (n + 1/2). We first observe that if Lη (L) → N 0 remains finite when we take L → ∞, then As our spectral representation involves N ∝ L sums of this kind we obtain an infinite product over factors that are strictly smaller than 1 and obtain a vanishing answer. Hence retaining only a finite number of p j 's in our sum is clearly insufficient.
Conversely, if η(L) = η stays finite we do have (90) at leading order in time because the terms we drop compared to having limits ±∞ contribute to an oscillatory integral of a bounded function and ∞ η e iθx obtained from the expansion of the exponential integral E 2 (x) = ∞ 1 e −xu u 2 du. If we want this term to become negligible compared to (40) at late times in the scaling limit, then we need Lη n+2 (L) → ∞ for all n in order to ensure that both the Euler MacLaurin correction terms in (93) and (94) are negligible. Hence any power-law η(L) = L −ν with ν > 0 will not suffice. Stated differently, the number Lη(L) has to be larger than any L ν for 0 < ν < 1, but any macroscopic fraction ǫL with ǫ > 0 is sufficient. We will denote by mesoscopic this number of states (in contrast with microscopic O(1) and macroscopic O(L)).
Finally, it is clear from e.g. (45) that in order to recover an exponential decay in time, one should multiply a O(L) number of terms like (40), and to recover the right exponent one should take into account 'almost all' these terms at each momentum q.
We conclude that if the sum over the p j in (30) is viewed in terms of particle-hole excitations over the q i , the leading late time behaviour emerges from a mesoscopic number (i.e. larger than any L ν for 0 < ν < 1, but smaller than ǫL for any ǫ > 0) of particle-hole excitations around each q i . It represents an exponential number of states, but still sub-entropic, in the sense that it includes only states whose macroscopic state is the representative state itself. We expect this to be a general feature of form factor expansions of semi-local operators that holds also in interacting theories.

Non-vanishing low-density limit of form factors
The leading order term both in time and density (45), obtained by keeping only the double poles in the partial fraction decomposition (36), has an interesting physical interpretation in terms of low-density limit of the form factor (19). This limit consists in assuming that ρ(x) is small everywhere, i.e. that L(q i+1 − q i ) → ∞ when L → ∞ in the scaling limit 3 . One can observe then that the square of the form factor (34) vanishes in the scaling limit L → ∞ unless each p is at a distance of order O(L −1 ) from a q, and because of the low-density assumption this q becomes unique in the limit L → ∞. Thus one can write with n i an integer. In the low-density limit we have and so in the scaling limit | q 1 , ..., q N |σ x l |p 1 , ..., which constitutes the low-density approximation of the form factor in the case it is not vanishing in the scaling limit. Now, we trade the 1/N ! in (30) for an ordering p 1 < ... < p N and choose an ordering q 1 < ... < q N . Such a constraint translates into a constraint on the n j in (95), namely n j+1 > n j − L 2π (q j+1 − q j ) ≈ n j − 1 2πρ(q j ) . At leading order in ρ, this constraint becomes − 1 4πρ(q j ) < n j < 1 4πρ(q j ) , so that (30) becomes Hence in the low density limit ρ(q) → 0 we have This expression factorizes and can be computed along (45) with formulas (155).

Quantum quench dynamics beyond low densities
The framework presented in section 3 is general and permits us to compute the late time behaviour and subleading corrections of form factor sums, order by order in ρ Q = N/L. In particular, it yields the expression of the observables of interest as Ce −t/τ , where both C and τ are the exact expressions at a given order (here second order) in ρ Q . The computation of a generic order in the density is however rather involved. In this section and in the following one we focus on the exponent τ of the exponential decay, for which another more efficient but less general approach can be used to calculate it at all orders in ρ Q , i.e. writing the correlation function as e −t/τ where τ includes all orders in ρ Q . This first section treats the quantum quench problem introduced in subsection 2.1.

Determinant representation
As shown in section 3.2, the leading contribution in (27) is obtained for N = M , on which we will focus. The starting observation is that the last term of (78) can be written as a Cauchy determinant, indeed we have (100) Let us defineM = C T C and M j ik = C ij C kj , so thatM ik = j M j ik . The determinant ofM can be expanded as follows The whereas the sign of the permutation changes. Hence the sum over τ vanishes unless all the j's are distinct, i.e., if they are a permutation of 1, ..., N . In conclusion we find ik . This relation permits one to eliminate the 1/N ! factor in (79), which then reads (103) where f (p) is defined in (80) and M is explicitly given by
The first two properties (i) and (ii) follow immediately from the definition (20). The third one (iii) means that if some p's are set to some q's in a one-to-one fashion, then one recovers the same function g with the remaining p's and q's. As for property (iv), it means that (ii) holds at order (p i − q σ(i) ) 2 .
We will now argue that by virtue of these properties setting g to 1 does not affect the leading behaviour at late times t. First, property (i) ensures that g does not modify the general structure (36) by allowing e.g. for higher order poles, or poles at other momenta. Property (ii) ensures that g does not modify the values A(I 0 , I 1 , I 2 ) when I 1 = {}, because in these cases there is always a double zero to differentiate and g is then evaluated at a permutation of the momenta. When I 1 = {}, the function g does change A(I 0 , I 1 , I 2 ), but, because of property (iv), it will not modify the pairing structure of (58), and will only modify the factors in (58) by an extra additive term. Finally, property (iii) allows one to repeat these steps recursively in (53). We now observe that the resulting partial fraction decomposition will always boil down to evaluating sums of the form (170) and (178). These are divergent as L → ∞, cf.
(172) and (181), and hence can only appear together, which cancels the divergence. Sums of type (170) emerge from I 2 , which is not affected by g, and the ones of type (178) come from I 0 , which is instead affected by g. Since the t dependence is carried by (170), the leading time behaviour will never depend on g.
Based on these arguments we now make the approximation of setting g = 1 where the matrix A is given by e it(2ε(p)−ε(q i )−ε(q j )) sin p sin q j f (p)/f (q j ) (cos q i − cos p)(cos q j − cos p) , i, j = 1, . . . , N . (108)

Asymptotic forms of the matrix elements
In the next step we work out the large-L asymptotics of the matrix elements A ij .

Diagonal matrix elements
The diagonal matrix elements were already computed in (83), cf. Appendix A. They read where we have defined

Off-diagonal matrix elements
To compute the off-diagonal elements, we write 1 (cos q i − cos p) (cos q j − cos p) and use the first equation in (83), see also appendix A. We obtain

Approximate determinant representation
Combining (109) and (112) provides the following approximate determinant representation where

Evaluating the determinant
We now write det M = exp tr log M and use the expansion (115)

First order
The first order gives

Second order
where These sums are computed in appendix A. We obtain (119) Once exponentiated in the determinant, we recognize the terms obtained earlier with the partial fraction decomposition approach (88), (89), but without the contributions involving the ε in the prefactor C. This difference is a direct consequence of our approximation g = 1.

Leading late-time contribution at all orders
In order to compute the O(t) term at higher orders in ρ, we first notice that it can only arise from the second order poles in the matrix entries, hence by pairs of momenta separated by o(L 0 ). In this regime the matrix elements become with ∆ ij = i − j. We obtain in this regime (121) The sums over ∆ j can be carried out using that for ∆ ′ = 0 which is obtained from and relations (155) by carefully treating the cases ∆ = 0, −∆ ′ . Using that θ = O(L −1 ) to neglect the |θ| term in (122) we arrive at

Influence of the boundaries
In the discussion above the momenta p are constrained to be positive. In order to use eqns (40) we therefore had to neglect possible boundary effects for q's close to zero. We now verify that this does not influence the result. We have with E 2 (x) = ∞ 1 e −xt t 2 dt the exponential integral function. Lη 1,2 are the number of vacancies between 0, π and q i . They are η 1 = q i 2π and η 2 = π−q i 2π . Hence the correction to tr Ξ is which goes to zero for large t because the density vanishes quadratically at 0 and π.

Result: late-time asymptotics of the order parameter after a quench
Substituting (116), (119) and (123) into (115) we arrive at the following result for the late-time asymptotics of the order parameter one-point function after a quench within the ferromagnetic phase with C given in (89) at order O(ρ 2 Q ). The decay rate reproduces the exact result obtained in [44,45]. However, the prefactor C differs from the one conjectured in [44,45]. We address this difference in subsection 4.5 below.

Numerical Checks
We now present some numerical checks of eqn (126) for the time evolution of the order parameter. In the limit of large separations ℓ ≫ 2Jt, using the Lieb-Robinson bound and the clustering properties of the initial state |Ψ , the two-point function factorizes into the product of two one-point functions that are identical by translational invariance with γ a constant of order 1. We can then obtain the one-point function Ψ|σ x 1 (t) |Ψ as the square root of the two-point function σ x ℓ+1 (t) σ x 1 (t) in the limit ℓ ≫ 2Jt, which can be efficiently computed numerically, as it can be expressed as the determinant of a block Toeplitz matrix even in the thermodynamic limit [45].
Numerical checks of the exponential decay in (126) have already been reported in [45]. Since our prediction (89) differs from the one conjectured in [45], we will focus on the prefactor C x FF (α) of the asymptotic behaviour of the two-point function in the limit ℓ, t → ∞, α = t/ℓ fixed In [45] it was assumed that the constant C x FF (α) is independent of α. Calculating the asymptotics of the correlator for α → ∞ then leads to [46] From (89) it however follows C = C x FF (∞), suggesting in turn that C x FF (α) is in fact αdependent. This is indeed supported by our numerical results, even though the difference |C − C x FF (∞)| is tiny. In Figs 10 we show that our results (126) and (89) are in agreement with numerical calculations of the order parameter one-point function. 40 60

Dynamical correlation functions at arbitrary finite temperatures
In this section we follow the same reasoning as in section 4 but for dynamical correlation functions at finite temperature. We again treat the two cases h < 1 and h > 1 separately.

Ordered phase h < 1
In the ordered phase, the sum in (30) involves states with the same number of particles as in the quench case (27). However, it is not possible to express each term as a Cauchy determinant as in section 4. This can however be overcome for the leading late time behaviour, where one can work with an approximate version of the form factor (19). We focus again on intermediate states with M = N in (30), which were argued in section 5.2 to give the leading late time behaviour. As discussed in section 4.1.1, the dominant contribution to the correlation function arises from the sums (170), and the term proportional to t on the right-hand side of (170) has its origin in the isolated singularity 1 x 2 in 1 sin 2 x . Hence, replacing terms of the form 1 sin 2 x by 1 x 2 in the form factor (19) will not affect the leading behaviour at late times. It follows that where we have defined The form factorF has the same Cauchy determinant structure we encountered in the quench case (100). Hence we can follow through the same steps and obtaiñ The leading time behaviour of the coefficients of this matrix can be computing as in section 4.2 with formulas (155). One obtains with Since one is interested only in the late time dynamics, we can focus on terms i, j such that q i − q j = o(L 0 ) as in (120). Then one obtains the same formula as in (120)

Two-point dynamical correlation functions in the ordered phase
Following through the same steps as in the quench case we arrive at where, using the results of section 3, dxdy in the space like-region at order ρ 2 β ξ in the time like-region at order ρ 0 β .
As far as we know the result (137) has not previously been obtained in the literature.

Numerical checks
In order to check the accuracy of (137) at finite times we have carried out numerical simulations following Ref. [76], where the two finite temperature dynamical two-point function for a finite open chain is computed exactly as a Pfaffian of a known matrix. As long as the two points are sufficiently far from the boundaries, then they take almost the same values as in an infinite chain. In Fig. 11 we compare (137) to numerical results in the space-like region by considering the logarithm of the correlator For simplicity we take the extreme case α = 0, which corresponds to setting t = 0. We recall indeed that static correlations are also covered by our calculations, as explained in section 3.5.3. In the left panel we plot (137) as a function of distance ℓ and mark by red crosses numerical results obtained for a chain of L = 200 sites. We see that the asymptotic result (137) is in excellent agreement for α = 0. In the right panel we test the accuracy of the O(ρ 2 β ) value of the prefactor (138) by considering the quantity where C is given by (138). We see that Γ(T ) approaches 1 for small values of T , which means that the prefactor (138) is indeed correct to order O ρ 2 β . The linear increase in ρ β shows that the correction to our result for C is O ρ 3 β . We now turn to the time-like region. In Fig. 12 we compare numerical results for L(ℓ, t) as a function of t for different values of ℓ to (137). We recall that in the time-like region eqn (137) only gives the leading time behaviour, i.e. the exponent of the exponential decay. Hence only the slope of our analytic result for L(ℓ, t) has to match the numerics. This is indeed seen to be the case in Fig. 12. In the space-like regime, i.e. at sufficiently short times, (137) is again seen to be in very good agreement with the numerical results. for v max t > ℓ.
In order to have a more quantitative check on the exponential factor in (137) it is useful to consider the difference where χ xx (ℓ, t) is given by (137). If our exponential factor is exact, the difference should decay at late times more slowly than exponential, i.e. supposedly as a power law. In Figure  13 we plot (141) as a function of log Jt and indeed observe a linear behaviour. This confirms that the exponential factor in (137) is exact and moreover establishes that the subleading asymptotics is power law in time. In this case the partial fraction decomposition (66) of | NS q 1 , ..., q N |σ x l |p 1 , ..., p M R | 2 necessarily involves one ν j = 0. To obtain the p j dependence of the corresponding A's in (66), it suffices to decompose the form factor | NS q 1 , ..., q N |σ x l |p 1 , ..., p M R | 2 by starting to write the partial fraction decomposition with respect to p j and retaining only the ν j = 0 part (which corresponds to P 0 (X) in (33)), and decomposing with respect to the other momenta in the usual fashion. This part is . This latter factor satisfies the hypotheses of (106) which allows us to replace it by 1 as long as we are interested only in the leading time behaviour. It follows that under this approximation the ν j = 0 terms in the partial fraction decomposition (66) of | NS q 1 , ..., q N |σ x ℓ |p 1 , ..., p M R | 2 contribute to (30) as Taking into account all the possible j = 1, ..., N + 1 for which one can have ν j = 0 in (66), we obtain where {p} = {p 1 , ..., p N } has now the same number of momenta as in the representative state. The second factor in (143) is precisely of the same form as the one we considered in the ordered phase.

Two-point dynamical correlation functions in the disordered phase
Putting everything together we obtain the following result for the leading late time behaviour of the two-point function where (cf. section 3.2) l, t ≥ 0, and v max is defined in (32).

Numerical checks
We have checked the accuracy of (144) by comparing it to numerical calculations following Ref. [76]. In the following we show results for In the space-like region v max t < ℓ we furthermore check our result for the prefactor C(ℓ, t) (145) by computing If and only if the prefactor in (145) is correct at order O(ρ β ), C(ℓ, t) → 1 when T → 0.
In Fig. 14 we compare our analytic expression (144) in the space-like region to numerical results for R(ℓ, t). Our analytic expression is seen to be in good agreement with the numerical results and the remaining discrepancy is due to O(ρ 2 β ) corrections to the prefactor. In Fig. 15 we present results for R(ℓ, t) in the time-like region v max t > ℓ at low temperatures. We see that the analytical result is in excellent agreement with the numerics.  In order to check the accuracy of our result for the exponential decay of the two-point function for intermediate and high temperatures we compare (144) to numerical results ln |R(ℓ, t)| in Fig. 16. As our result for the prefactor C(ℓ, t) only holds at low temperatures we expect the numerical results to differ from our analytical prediction by an essentially constant offset. This is indeed what we observe.

Summary and Discussion
In this work we have considered two problems in the transverse field Ising model: (i) The time dependence of the order parameter after a quantum quench; (ii) The dynamical order parameter two point function in equilibrium at finite temperatures. Using the quench action approach for (i) and a micro-canonical formulation combined with typicality ideas for (ii) these problems can both be formulated in terms of spectral representations using Hamiltonian eigenstates, i.e. sums over form factors. These highly intricate sums over a macroscopic quantity of momenta of a form factor with many singularities represent however a considerable technical challenge. We showed that in the case of semi-local operators such as σ x j in the TFIM, this difficulty can be addressed by decomposing the form factor into partial fractions, which permits the sums over momenta to be decoupled and the late time behaviour to be determined. These partial fractions can be organized in terms of the degree and the position of theirs poles, which naturally leads to an expansion in the density of particles of the representative state. The leading behaviour at late times can be then computed at all orders in the density through a determinant representation of these poles, which leads invariably to an exponential decay.
Our analysis provides a precise characterization of the excitations over the representative state that contributes to the late time behaviour of correlation functions of semi-local operators. We find that simultaneous particle-hole excitations of all particles in the representative state contribute to the correlation function, and we altogether have to sum over O (ǫL) N excited states, where ǫ is a fixed number that can be taken as small as desired and N/L is the density. In particular, this implies that the appealing picture of an expansion in terms of a finite number of particle-hole excitations over the full momentum space fails for semi-local operators at finite temperatures. The form factor sum is dominated by mesoscopic excitations (in the sense given in section 3.5.1) around each single particle, which is an exponential number of states, but subentropic in the sense that includes only states whose macroscopic state is the representative state itself.
We have compared our analytic results to numerical computations in both the finite temperature and quench contexts. In the absence of saddle points, thus where the sums (40) are valid at all momenta, we find remarkable agreement at all times. In cases where saddle points occur our numerical results indicate the presence of a multiplicative power-law behaviour as a subleading correction to the exponential decay. We believe that this will emerge from saddle point effects of these classes of excitations, due to the fact that (40) is not valid anymore close to the saddle point. Higher-order corrections in time will also arise from contributions with ν i = 0 in (36); in this case there is no singularity and the full momentum space should be summed over to take them into account. We believe that partial fraction decompositions for form factors of semi-local operators are not only suited for extracting the late time asymptotics, but should also be useful for determining intermediate and short time behaviours.
In interacting models, the singularity structure of the form factors of semi-local operators is not fundamentally changed and a partial fraction decomposition will provide a useful organizing principle as well. A notable difference is that the Bethe equations link the different particles so that the momenta sums can never fully decouple. This will be the subject of a subsequent paper.
For local operators however, the story is radically different. The singularity structure of the form factors (12) is completely dissimilar and one cannot use partial fraction decomposition related to an expansion in powers of particle density to extract the late time behaviour. Excitations over the full momentum space and not only near the singularities will play equally important roles. This will be discussed elsewhere.

A Riemann sums of singular functions
Sums of form factors lead to Riemann sums of functions with a quadratic singularity of the Since the integral of f (x) is divergent, the limit L → ∞ cannot be directly taken as for regular functions and needs special treatment. The purpose of this appendix is to explain techniques to compute them.

A.1 Oscillatory sums
Oscillatory Riemann sums of functions with a quadratic singularity cannot be estimated with the usual results on stationary phase approximation to obtain their large time behaviour. The principle is then to add and remove an elementary function with the same singularity, but whose Riemann sum can be computed directly, so that the remaining function has no singularity and thus has an oscillatory Riemann sum that vanishes at large times.

A.1.1 Generic example
For concreteness, let us illustrate this procedure with the Riemann sum for f (x) a function such that f (x) = 1 The second term on the right-hand side can be turned into an integral, while we rewrite the first as follows (151) The second and third Riemann sums can now be turned into integrals without any problems, which gives The three integrals are oscillatory integrals of bounded functions and hence vanish for large θ. We conclude that (153)

A.1.2 Elementary oscillatory sums with singularities
The above analysis shows that the leading asymptotics of Riemann sums with simple or double poles involves the following sums n∈Z e iw(n+α) (n + α) = π sin πα e iπα sgn (w) , These are readily obtained by computing the Fourier series coefficients of the right-hand sides multiplied by e −iwα and seen as a 2π-periodic function of w. In particular we have n∈Z e iw(n+1/2) (n + 1/2) m = π 2 1 − |w| (40) Following the steps outlined above we obtain for q ∈ NS

A.1.3 Equations
At leading order in time, one can Taylor expand the ε in the remaining sum to fall back on an elementary oscillatory sum. The three integrals involve oscillatory bounded functions and hence decay to zero with time.
Similarly we obtain for q ∈ NS where the integrals are again integrals of oscillatory bounded functions and vanish at late times.
Finally we use (155) to obtain the asymptotic values of the sums in (156) and (157)  (158)
(162) Taylor expanding ε(p) around p = q and using (155) we finally arrive at A.1.5 Two-dimensional sums In this subsection we apply the method described above to the calculation of the twodimensional sums arising in section 3.1.5. We first consider (cos q i − cos q j ) 2 e 2it(ε(q i )−ε(q j )) .
where for large L the q's are distributed according to the density ρ(q). Starting with the sum over q i and using that sin xf (x) (cos x−cos q) 2 − f (q) sin q(x−q) 2 − f ′ (q) sin q(x−q) is a bounded function of x, we conclude that Next we Taylor expand ε(q j ) around q j = q i and use to obtain where ϑ j = 2tε ′ (q j )/(Lρ(q j )). The sums over i can then be carried out using (155) Finally we turns sums into integrals and use (86) to obtain Next we consider Ω 2 (t) = 1 L 2 i =j e it(ε(q j )−ε(q i )) sin 2 q i −q j 2 sgn (ε ′ (q j )ε ′ (q i )) , The sum over j can be carried out by proceeding in the same way as for Ω 1 (t), which gives Finally we turn sums into integrals to arrive at our final result

A.2 Non-oscillatory sums
We saw in the previous section that the non-decaying part of an oscillatory Riemann sum with singularity only depends on the singularity. However, if it is not oscillatory the Riemann sum will depend on the full function and its calculation is a bit more subtle. The principle is to separate the sum over two regions, one that is over points at distance less than AL from the singularity, and the other one at distance larger than AL, for a fixed arbitrary A > 0. By adding and removing the corresponding singularity in the first sum, one can take the limit L → ∞ and obtain the result as a seemingly A-dependent but in fact constant expression. To obtain a simpler expression for this constant one then typically looks at the limit A → 0 of the A-dependent expression.

A.2.1 Generic example
Let us again illustrate the procedure with an example: for f (x) a function such that f (x) = 1 x 2 + O(1) for x → 0, f being regular otherwise. Let us fix A > 0 and decompose with f (x) = f (x) − 1/x 2 a regular function. Using B k=1 k −2 = ∞ k=1 k −2 − B −1 + O(B −2 ), one obtains at fixed A and L → ∞ Then, writing f (x) = (x 2 f (x))/x 2 , integrating by part the first integral We note that the second integral involves a function that is bounded for |x| ≤ A and so vanishes for A → 0. Since the relation holds for any A > 0, in particular for A → 0 we get

A.2.2 Two-dimensional sums
We now apply these ideas to the calculation of in the space-like region, where sgn (ε ′ (q j )ε ′ (q i )) = 1 everywhere. We first decompose the sum into divergent and finite parts We evaluate the second term by splitting it into contributions from |q i −q j | > B and |q i −q j | < B. The latter vanishes in the limit B → 0, while the former becomes where we used in the first line BL k=1 k −2 = π 2 /6 − (BL) −1 + O(L −2 ), and where in the last line we have integrated by parts. Finally we take B → 0 to obtain Next we consider the sum We again separate it into divergent and finite parts l =j sin q l sin q j (cos q j − cos q l ) 2 − L 2 ρ 2 (q j ) (l − j) 2 = 8π 3 3 j ρ 2 (q j ) + 8 L 2 l =j sin q l sin q j (cos q j − cos q l ) 2 − L 2 ρ 2 (q j ) (l − j) 2 .
We evaluate the second term by splitting it into contributions from |q l −q j | > B and |q l −q j | < B. The latter vanishes in the limit B → 0, while the former becomes where in the last line we have integrated by parts. Finally we take B → 0 to obtain