Quantum Quenches in Integrable Field Theories

We study the non equilibrium time evolution of an integrable field theory in 1+1 dimensions after a sudden variation of a global parameter of the Hamiltonian. For a class of quenches defined in the text, we compute the long times limit of the one point function of a local operator as a series of form factors. Even if some subtleties force us to handle this result with care, there is a strong evidence that for long times the expectation value of any local operator can be described by a generalized Gibbs ensemble with a different effective temperature for each eigenmode.


Introduction
In recent years, the unitary time evolution of an extended quantum system has attracted a lot of attention. The interest of the community on this topic was spurred by novel experiments with cold atoms that have shown the coherent evolution of a quantum system in a laboratory [1][2][3][4]. One of the simplest and most studied setting is the so called (sudden) quantum quench. In this situation, one initially considers an extended quantum system prepared in a pure state: this can be regarded as the zero temperature ground state of an Hamiltonian H 0 . At t=0 the Hamiltonian is suddenly changed from H 0 to H and then the system is let to evolve unitarily according to the new Hamiltonian H, without any coupling to the environment. Sure enough, there is a transient regime for such an abrupt change of the Hamiltonian but the main interest is in what happens in the long time limit. A relevant issue is whether the expectation values of local operators at t = ∞ can be derived by standard thermodynamical ensembles (e. g. Gibbs ensemble): if this is the case, it is said that the system thermalizes. The usual understanding is that thermalization should generally occur in extended systems but this notion has been recently challenged by experiments done on one-dimensional quasi-integrable systems, such as trapped 87 Rb atoms, which do not noticeably equilibrate ever after thousands of collisions [2]. Although the system is in a magnetic trap, this may be a manifestation of the quantum integrability of the Lieb-Liniger Hamiltonian for the atoms placed in the flat region of the trap. It was firstly conjectured by Rigol et al. [5] that, in dealing with integrable systems, one should employ not the usual density matrix (function of the energy only) but, instead, the density matrix of a generalized Gibbs ensemble which involves all the integrals of motionÎ l of the systemŝ ρ gen ∼ e − l α lÎl . ( This conjecture has been the subject of several studies, mostly done by using specific models [6][7][8][9][10], but also studied for initial thermal distributions [11]. An important step forward was taken in the work of Barthel and Schollwöck [12] in which they generalized a previous result by Calabrese and Cardy [13], proving rigourosly that for Gaussian initial states and quadratic (fermionic or bosonic) systems the conjecture does hold: a finite subsection of an infinite system indeed relaxes to a steady state described by a generalized Gibbs ensemble where the extra integrals of motions are simply the occupation numbers of each eigenmode. Moreover, they stated (postponing the proof to a future paper) that a similar result holds also for Bethe-ansatz solvable models, even if in this case the extra integrals of motions appearing in the density matrix do not have such a simple physical interpretation.
In this paper, we tackle the issue of thermalization in integrable systems from a different point of view: we focus our attention on continuous systems, i.e. integrable field theories in 1+1 dimensions. A reason behind this choice is that the notion of quantum integrability can be put on a firm ground only using, as criterium, the elasticity of the scattering amplitudes of the excitations of the system [14] and, indeed, this is a distinguished property of integrable field theories [15]. Moreover, the rich analytic structure of integrable field theories could help us to obtain results that are harder to obtain with other methods. At equilibrium, quantum field theory is a very useful tool in order to describe the scaling limit of real world condensed matter systems. If such a limit exists also out of equilibrium is still an open question: in fact, after a quantum quench a large number of high energy state is populated [16,17], while a field theory usually describes only the low energy excitations of the corresponding lattice model. However, studying quantum quenches in field theory is still an interesting topic, since we expect that the qualitative picture emerging from this analysis should at least give precious hints on the behavior of real world systems. The study of quantum quenches in integrable systems has also been addressed in [18,19].
The structure of this article is as follows. In section 2 we briefly remind the main properties of integrable field theories and we properly define the kind of quenches we are interested in. The main result of this paper is stated in section 3, where we analyze the long time limit of the one point function of a local operator. Using a form factor expansion, we show that its asymptotic value could be obtained from a physically transparent generalized Gibbs ensemble, where the conserved charges are simply the occupation number of each eigenmode. Since the argument in section 3 may appear a little abstract, in section 4 we compute the exact time dependent one point function of the energy density operator of the Ising model, showing that its asymptotic value agrees with the general formula. Finally, we present our conclusions in section 5. In Appendix A we gather all the technical details of our derivation.

Basic facts about integrable quantum field theories
In this paper we focus our attention on integrable quantum field theories in 1+1 dimensions (for a review see e.g. [15], [20], [21]). These theories are characterized by an infinite set of local conserved currents that greatly constrain the dynamics. Indeed, only elastic scattering events can occur in these theories: no particle production is allowed and the sets of the final momenta coincide with the initial one. Moreover, due to Yang-Baxter equation, the amplitude of the n → n particle scattering event can be written in terms of the 2 → 2 scattering amplitudes. Using these remarkable properties of the scattering processes and some additional assumptions on the analytic structure of the S matrix, one is able to determine all the scattering amplitudes and, from their pole structure, the mass spectrum of the theory. For sake of clarity in the following we focus our attention on the simplest integrable field theories, i.e. those with only one excitation of mass m (e.g. Ising model or Sinh-Gordon Model, where the latter is directly relevant for one-dimensional bose gases [22,23]), even if our results could be easily extended to any integrable field theory with richer spectrum of excitations.
A convenient basis can be constructed by the action of a creation operator Z † (θ) on a vacuum state |0 , where θ is the rapidity of the particle (with energy and momentum given by E = m cosh(θ), p = m sinh(θ) respectively). Since the theory is not free, the operators Z(θ) and Z † (θ) do not simply commute or anticommute: they satisfies instead the Zamolodchikov-Fadeev algebra where S(θ) is the two particles S matrix. Consistency of this algebra is ensured by the relationship S(θ) = S(−θ) and by the unitarity condition S(θ)S(−θ) = 1. Except for the free bosonic theory, where S(0) = 1, the S-matrix of all other interacting theories satisfies S(0) = −1. Multi-particle states are given by acting with an ordered string of the creation operators Z † (θ i ) on a vacuum state |0 , i.e. where These orderings select a set of linearly independent vectors that form a basis in the Hilbert space. A convenient resolution of the identity is given by . . dθ n 2 π |θ 1 , . . . , θ n θ n , . . . , θ 1 | , where the integration is extended to all values of the rapidities. The structure of an integrable theory is so powerful that permits to compute also the matrix elements of any local operator O(x) on the basis (3). Indeed, by solving a set of monodromy and recursive equations (based on the S-matrix), one can determine the form factors [21] F O (θ 1 , . . . , θ n ) = 0|O(0)|θ 1 , . . . , θ n .
Consider, for instance, the Ising field theory at T > T c . This model can be described by a free fermionic theory with S-matrix S = −1 and only one kind of particle of mass m ∝ (T − T c ). There are three physically interesting operators in such a theory: the energy operator ǫ (that is proportional to the trace of the energy momentum tensor), the spin operator σ (the order parameter) and the disorder operator µ related to σ by the Kramers-Wannier duality. With a proper choice of the overall normalization, their form factors are given by the formulas [24] F ǫ (θ 1 , . . . , F σ (θ 1 , . . . , θ n ) = σ i (n−1)/2 n l<m tanh( θ l −θm 2 ) if n is odd 0 otherwise , where σ = 2 1 3 e − 3 4 A 3 m 1 4 and A = 1.282427 . . . is the Glasher constant †. It is clear that, from formula (5), we can easily restore the x and t dependance of the matrix element 0|O(x, t)|θ 1 , . . . , θ n , since the energy and the momentum operators are diagonal in the basis (3). However, it is not so obvious how to obtain from (5) the generic matrix element θ n , . . . , θ 1 |O(0)|θ ′ 1 , . . . , θ ′ n : this task can be accomplished by exploiting the crossing symmetry. If A and B are two sets of rapidity, we have where the sum is over all the possible ways of splitting the sets A/B in two subsets A 1 / B 1 and A 2 /B 2 while S AA 1 and S BB 1 are the products of S(θ) we need to rearrange the rapidities in the proper order, namely The symbol A + 1 in (9) denotes that each rapidity θ 1 . . . θ r in A 1 is shifted by an infinitesimal amount ǫ i so that A + 1 |O|B 1 is simply related to the form factors (5) Here comes however the tricky point. If the ǫ i are finite, the form factors are (for real rapidities) regular functions, while if we take the ǫ i → 0 limit (as we need to do at the end of the calculations) the form factors usually diverge, because we are in the kinematical situation in which some of the rapidities of the bra and the ket states coincide: the simplest case of this circumstance is provided by the 2-particle matrix element θ|µ|θ ′ where µ is the disorder operator of the Ising model (8), which indeed diverges when θ = θ ′ . This discussion shows that a prescription is needed for handling these kinematical divergencies. The one proposed in [25,26] consists of taking only the regular part of (12) and discarding all the terms proportional to an inverse power of ǫ i = Finite Parts lim It should be stressed, however, that this prescription alone is not enough to properly take care of all the divergencies and, usually, it must be supplemented with extra corrective factors coming from the Bethe-ansatz technique [25]. Without entering in too many details, to show this aspect, consider for instance as a density matrix of the system ‡ wheren(θ) = Z † (θ)Z(θ) and λ(θ) is an appropriate function of θ. Notice that, if λ(θ) = 1 T m cosh(θ) the density matrix (14) describes the familiar canonical ensemble, † Here we are using the conformal normalization of the operators, such as in the limit while, in the more general case, it can be associated to the generalized Gibbs ensemble (1): indeed, in an integrable field theory the (infinite) set of local continuity equations determines the conserved chargesQ s labeled by an integer index s [15], and each chargê Q s can be written aŝ where q s is a constant. Obviously, it would be nice to write down the average of a local operator O w.r.t. the density matrix (14) in term of these form factors. But, if we do it applying the prescription (13) alone we end up with the result which is wrong since it does not agree with the thermodynamic Bethe ansatz. It was firstly conjectured by LeClair and Mussardo [25] that the correct expression is instead where theλ are dressed according to the integral equatioñ ϕ being the derivative of the phase shift Originally, (17) was only an educated guess but its correctness was confirmed by subsequent checks [27][28][29]. Even if presently there is still no complete proof of the exactness of (17), a step forward its rigorous derivation has been taken in the recent papers by Pozsgay and Takacs [30][31][32], who showed that (17) might be obtained in a rigorous way by putting the system in a finite volume, thus discretizing the rapidities to regularize the kinematical singularities, and taking the infinite volume limit only at the end of the calculations.

Integrable boundary condition and quantum quenches
Our aim is to study the time evolution of the expectation value of a local operator O(x) on a pure state |B that is not an eigenstate of the integrable Hamiltonian H, i. e.
As shown by Calabrese and Cardy [13,33], performing a Wick rotation, this dynamical problem can be mapped into a statistical problem defined in a strip geometry, where the initial state |B plays the role of a boundary condition. In an integrable field theory, the most natural boundary states are the ones that do not spoil the integral of motions. These integrable boundary states were originally studied by Ghoshal and Zamolodchikov in [34] and their properties can be summarized as follows: these states have a peculiar form, given by a coherent superposition of pairs of equal and opposite rapidity (Cooper pairs) § where the indices a and b run over all the particles. The amplitude K ab (θ) satisfies a set of equations that depends on the S matrix (such as boundary Yang Baxter, boundary unitarity and crossing equations): the different solutions of these equations provide the set of integrable boundary conditions of the theory in question. As an explicit example, let's consider once again the Ising model, in which there are three integrable boundary conditions associated to the amplitudes where k = 1 − h 2 2m and h is a real number (the boundary magnetic field). Clearly, the magnetic boundary condition interpolates between the free and the fixed ones. So, we see that only particular choices of the amplitude K give us an integrable boundary state.
Our interest in the following is to study a particular class of quenches, i.e. those where the initial state is expressed by a coherent superposition of particle pairs, as the one given in eqn. (21). Among these quenches there are, of course, the integrable boundary states but the whole class of these quenches is actually larger than the genuine integrable ones: indeed, for theories with only one particle, it suffices that the amplitude K(θ) satisfies the "boundary cross-unitarity equation" This equation simply comes from the request of invariance of the expression (21) under a change of variable θ → −θ.
Notice that, at least for a free fermionic theory as the Ising model, a boundary state of the form (21) has a very simple meaning within the theory of quench processes. Suppose in fact that H 0 is a free fermionic theory, with its ground state |B identified by the condition Z 0 (θ)|B = 0, where Z 0 (θ) is the annihilation operator related to H 0 . If, at t = 0, the Hamiltonian is suddenly changed from H 0 to a new almost free Hamiltonian H (for instance, making a quench of the mass from m 0 to m), the new sets of creation/annihilation operators Z † , Z will be related to the initial ones by a Bogoliubov transformation of the form § Here, for simplicity, we ignore the presence of additional zero-rapidity terms in the expression of |B . For theories with many particles, K ab (θ) must also solve an appropriate boundary Yang-Baxter equation to ensure commutativity of the various terms in the exponential expression (21).
Hence, the algebraic equation that identifies the boundary state in terms of the new creation/annihilation operators is given by whose solution is Summarizing, in the following we consider those quantum quenches where the initial state is given by a superposition of Cooper pairs according to 21), where the amplitude K is a (regular) function that satisfy (23). However the boundary state as written in (21) is typically not normalizable. Indeed, if the amplitude K does not goes to zero for large rapidities (as in eq. (22)), it means that we are exciting modes with arbitrary high energy, hence the state has an unphysical divergent energy density. So, as done for the critical systems [13,33], it is therefore convenient to introduce an extrapolation time τ 0 that make the norm of the state (21) finite: this parameter plays the role of an ultraviolet cut-off, in any case present in physical systems. So, we are interested in where

Long time behavior
In this section, we are concerned with the computation of and in the following we will show that O can be expressed as an average over a density matrix like (14). At first sight, it seems very unlikely that we can do so: our boundary state (28) has a very peculiar structure, since it is the superposition of pairs of opposite rapidity, while in an average over an ensemble (17) there is no sign of such a structure. How comes that the system retains no memory of this pair structure in the long time limit? First of all, a trivial remark: since |B is translational invariant, O(x, t) B do not depend on x and therefore, from now on, we will set x = 0. In principle, it is quite clear what we have to do in order to compute (30): first we expand the exponential in (28) and, taking into account that the Hamiltonian is diagonal in the particle basis, we thus arrive to the double sum where we used the short-hand However it is difficult to compute the long time limit directly from (31). The reason is that the matrix element in (31) are not regular functions (rather they have delta-like contributions) and, in such a case, we cannot apply a stationary phase argument. In order to isolate the singular parts, we have to employ eqn. (9). Consider, for instance, the term with n = l = 1 in the numerator of (31), i.e.
Applying the crossing relation (9), we can recast this term as where we make use of the symmetry properties (29) of G. Actually, the inner product θ 1 , −θ 1 | − θ ′ 1 , θ ′ 1 above is divergent in the infinite volume limit and it should be regularized by putting the system in a box of length L, obtaining This, however, is not an important contribution since, as shown in Appendix A, all these inner products cancel out with the denominator of (31). What is really crucial is that, apart from these infinite volume divergencies that we can easily regularize, the integrands in (34) are all well-behaved functions: hence we can now easily take the infinite time limit t → +∞, so that (34) simply becomes because the first term in (34) vanishes for the fast oscillation of its integrand.
In the light of this example, the strategy to compute the expectation values of local operators can be stated as follows.
(i) We first expand the exponential in the numerator of (28), ending up with the double sum (31).
(ii) Then, we use (9) in order to isolate the delta-like terms and, after having done that, we take the infinite time limit, where all terms that explicitly depend on time go to zero, due to the fast oscillation of the integrand. This is a simple consequence of the stationary phase argument, that can also be seen in the following way. If the infinite time limit exists (and the stationary phase argument assures us that it does exist), then it must coincide with the temporal average We know that the time-dependent part of the numerator of (31) consists of a sum of terms as where F is the regular function obtained by applying (9). So, the only contributions to the infinite time limit comes from the region E n (θ) = E l (θ ′ ), whose Lebesgue measure is zero, so the integral goes to zero since F has no delta-like term.
With these steps in mind, it is a nice combinatorial exercise to show that O can be finally expressed as In order to keep our exposition clear, we will present the details of the combinatorics behind (39) in the Appendix A. However, (39) is still a meaningless expression, since we have to regularize it in a proper way. One way to do it is by analogy with LeClair and Mussardo formula, discussed in section 2.1. When we perform an average over a density matrix (14), we end up with the following expression LeClair and Mussardo suggested that the proper way to regularize the ǫ i → 0 limit of this expression is to take the the connected part of the form factors (13) and to dress λ(θ) according to the integral equation (18). This regularization scheme holds for every function λ(θ). The situation is the same in equation (39), with |G(θ)| 2 that plays the role of e −λ(θ) . So, the natural way of regularize (39) lead us to where |G (θ)| 2 is dressed in the same way as the term e −λ(θ) entering the thermodynamic Bethe ansatz The above dressing formula is based on the LeClair and Mussardo conjecture (it has actually the same mathematical structure) and on the possibility to exchange the ǫ i → 0 and t → +∞ limits. Assuming that such a regularization scheme is indeed correct, it turns out that that the long time limit of (27) could be described by a generalized Gibbs ensemble (14), where the constant of motions are simply given by the occupation numbern(θ) and the function λ(θ) is fixed by the conditions B |n(θ)|B B |B = T r (ρ λn (θ)) , thus proving Rigol et al.'s conjecture for integrable field theory. If we look at our starting point (27), this result is quite unexpected: when we expanded the exponential in (27) we had a double sum, and it was only thanks to the infinite time limit that we could rewrite it as a single summation. Moreover, while the boundary state (28) is formed by pairs of particles with opposite rapidity, this feature is completely lost in the final expression (41).

The simplest example
It is instructive to see how the general ideas of the previous section apply in the simplest case provided by the one point function of the ǫ operator of the Ising model (6). Indeed, for this operator we can calculate exactly its one point function for any time with elementary techniques. From the form factors (6), it follows that the operator ǫ is a quadratic form in the creation -annhilation operators Since the theory is free, we can easily calculate the expectation value of binomials of the creation-annihilation operators on the boundary states (introducing, for instance, a generating functional) and we have that Hence ǫ(x, t) B = + 2 π m dθ 2π The results (47) has all the features we expect to hold in the general case. First of all, the time dependent part goes to zero as a consequence of the fast oscillation of the integrand. The long time asymptotic value obviously agrees with our general result: in this case, the structure of the operator is so simple that the entire sum reduces to a single term. More important, since we are able to calculate exactly the time dependence, we can also estimate the approach to the t → +∞ limit value using a stationary phase approximation. It turns out that (47) approaches its asymptotic value as an inverse power law, in contrast with the exponential decay of the massless case [33] [13]. Heuristically, the different behavior of the two cases (massive and massless) can be understood in the following way. We know that the relaxation to the stationary value is due to the speed of the quasiparticle excitations: in a massless theory, all the quasiparticles move at the fastest possible velocity, so we can expect that the decay to the limit value is the fastest possible (exponential decay), while, in a massive theory, the time decay has to be slower for the spread of the velocity distribution of the quasi-particles. It may be interesting to compare this continuum formula (47) with the one point function of the transverse magnetization of the quantum Ising model [35]: a quench of the transverse magnetization in the lattice corresponds to a mass quench in the continuum. In both cases there are oscillations modulated by a power law with exponent − 3 2 . However, it must be noted that the long time leading behavior of the transverse magnetization on the lattice is not only dictated by the small momenta, and therefore it is not a surprise that the continuum theory is unable to capture all the features of the lattice model.

Conclusions
In this paper we studied the quantum quenches in the context of integrable field theories in 1+1 dimension. We considered a specific class of quenches, those whose initial state is a coherent superposition of Cooper pairs (28), a structure that is in agreement with the integrability of the theory in the bulk. For such quenches, we analyzed the one point function of a local operator and argued that, in the long time limit, it can be neatly expressed in terms of the formulas (41 ) and (42): the final result employs a series over the connected form factors (13), which are finite expressions.
As discussed extensively in the text, our final result must be presently considered only as a solid conjecture since it is based on a regularization scheme of the kinematical singularities of the original form factors. This regularization is closely related to the one used in LeClair and Mussardo formula for the one point function at finite temperature [25] -formula that has successfully passed many checks although it is not yet rigorously proved. In this respect, we have performed an additional but independent test of our result for a particular one point function of the Ising model. If the full correctness of eqn. (41) will be confirmed by a future rigorous analysis, it follows that the asymptotic value of the one point function of a local operator could be computed as an average over a density matrix (14) with a temperature different for each eigenmode: in short, our result is a strong evidence that Rigol et al.'s conjecture does hold for integrable field theories.
In this paper we limited ourselves to integrable theory with only one massive particle but our results could be extended to theories with more than one particle without any physical difficulty, although with a combinatorics much more involved. Finally, our results seems to suggest that in a generic massive integrable theory the decay towards the asymptotic value is dictated by a power law (in contrast with the conformal exponential decay) whose leading exponent can be determined by a stationary phase approximation.

Acknowledgments
We would like to thank G. Brandino, M. Burrello, P. Calabrese, E. Canovi, G. Gori, T. Macrì, D. Rossini, A. Silva and S. Sotiriadis for many fruitful discussions. This work is supported by the grants INSTANS (from ESF) and 2007JHLPEZ (from MIUR).

Appendix A. A diagrammatic proof of eq. 39
In this appendix we would like to show that, for t → +∞, In this appendix, we will call the form factors like θ + n . . . θ + 1 |O(0)|θ 1 . . . θ m regular form factors, in order to distinguish them from the complete form factors θ n . . . θ 1 |O(0)|θ 1 . . . θ m , that have also delta-like contributions. For ǫ i finite these regular form factors are continuous functions but, despite their name, they can have a singular ǫ i → 0 limit, hence the need of the regularization procedure previously discussed. Let's firstly briefly anticipate the main steps of the proof: after we expand the numerator of the l.h.s. of (A.1) (as done in (31)), we obtain B |O(x, t)|B = +∞ n,l=0 Then, applying repeatedly the crossing relation (9) to the matrix elements in (A.2), we arrive to an expression in which we can take the infinite time limit. In this limit, (A.2) reduces to the r.h.s. of (A.1) times the denominator of the l.h.s thus completing the proof. Figure A1. In this figure we present the building blocks of our diagrammatic representation. On the left, we can see a bra and a ket with n = 3 Cooper pairs. In the middle, we see one of the diagrams that represent the inner product of the bra and the ket. The figure on the right shows a term proportional to the regular form factor θ + 2 , θ + 1 |O|θ ′ 1 , θ ′ 2 . The circle in the middle stands for the operator O, that has 2r = 4 legs, half connected to the bra and half to the ket.

Introductory remarks
Here we would like to highlight some basic facts that we will find useful for our proof. First of all, we notice that, when we recast (A.2) in terms of the regular form factors, the only terms that survive in the infinite time limit are the time-independent ones. As a consequence, we can discard all the terms in the double sum in (A.2) where n = l. Moreover, from (9) it is clear that the term 0|O|0 in the r.h.s. of (A.1) is correct, so in the following we will not consider this contribution.
In order to follow more easily our ideas, it may be useful to develop a diagrammatic representation. In fig. A1 it is shown a bra state with n pairs (a ket state can be introduced in a similar way). An inner product between a bra and a ket both made of n pairs can be represented as the sum of all the possible ways to link together the particles of the ket with the particles of the bra, each link meaning a contraction between the corresponding particles. Of course, one should be careful about the permutation of particles and the corresponding S matrix, but we will take care later of these details. We are interested in the matrix elements of the operator O between states made of Cooper pairs. In particular, our aim is to reduce the full matrix element θ n , −θ n , . . . θ 1 , −θ 1 |O| − θ ′ 1 , θ ′ 1 , . . . − θ ′ n , θ ′ n to the regular form factors θ + ir , . . . , θ + i 1 |O|θ ′ j 1 , . . . , θ ′ jr . These regular terms can be diagrammatically represented as in fig. A1, where the operator O has 2r legs: r of them are connected to the particles in the bra, while the other r are connected to the particles in the ket. So, our combinatorial problem reduces to a problem in which we have to connect the legs of the operator to the bra/ket and the remaining particles together, according to (9). In order to get the right combinatorial coefficients, we have to remember what follows.
• As we stated before, only the matrix elements with the same numbers of particles in the ket and in the bra give a non vanishing contribution in the long times limit. Therefor, the number of legs connected to the bra is always equal to the number Figure A2. Two different diagrams: in the left one all the particles are connected to the operator, while in the right one there is a contribution that factorizes.
of legs connected to the ket and, from now on, we will only specify the number of particles connected to the bra.
• In principle, when we consider the matrix element between two states with n pairs, i. e. θ n , −θ n , . . .
we could expect to end up with a sum of regular terms with r legs linked to the bra, with r ≤ 2n. However, it turns out that r ≤ n: if we connect r particles to the operator, we are left with 2n − r delta functions, and to suppress all the time dependencies, we have to eliminate n integration variables, hence r ≤ n.
• Finally, if we link a particle to the operator, its pair partner cannot be connected to O, otherwise their time dependence survives.

The disconnected terms
In this subsection we will show how some contributions to (A.1) (the disconnected terms) cancel out with the denominator Z (A.3). In order to understand this point, we analyze the matrix elements As it is shown in fig. A2, we have essentially two types of diagram. Let's focus our attention on the second one: it is clear that the contribution in the dotted box factorize from the integral with the form factors. Hence, if we have a set of particles that is completely disconnected from the operator O (this means that no one of these particles or their partners is connected to O or to a particle whose partner is connected to O), its contribution factorizes.
What we want to show now is that these disconnected pieces cancel out with the denominator Z. Let us consider the term in (A.2) with n pairs in the bra (let's denote it as O n ), and let's focus our attention on the contribution O [n,k] , such that k pairs in the bra (as well as k pairs in the ket) are disconnected and so only n − k pairs in the bra are connected to O. Of course, a similar term can also be obtained from O n−k , when no particle is disconnected. With our notation, this term can be written as O k! 2 is the right one for Z k , thus concluding the proof of A.5.

The last step
Here we conclude our proof, showing that O [m,0] has actually the right structure to give (A.1). Let's call O {m,r} the contribution where the operator has r legs connected to the bra, the bra has m pairs and no particle is disconnected from O. We have already pointed out that for t → +∞ where the summation i 1 ,...,ir ′ is over all the positive integers i j such that j i j = m.
In order to prove (A.8), we need to be a little careful with the ordering of the particles and the labeling of the rapidities. However, if we exchange two particles, the contribution is same, since (as we already know) the pairs do commute while the exchange of two particles forming a pair is equivalent to change of the integrable variable θ → −θ. This is a consequence of the symmetry of G (29). Before concluding our proof of (A.1), we need to understand how to label the rapidity. We start with 2m integration variables and the delta-functions reduce them to r. So, we use the convention to label the rapidities as in (A.8): we call θ 1 the rapidity of the particle of the bra closest to O and we take advantage of delta functions in such a way that the rapidity of the particle in the ket nearest to O is also θ 1 , and so on.
We can now finally show that we obtain exactly the structure (A.8). First of all, from our previous reasoning about the long time limit, it is clear that, performing all the the exchanges done. However, if we remind that at the end of our calculation we have θ i = θ ′ i for i = 1, . . . r, it is easy to see that this product of S matrices reduces to S(0) . We can repeat this procedure until the overall coefficient is (S(0)|G(θ 1 )| 2 ) i 1 −1 . Then we contract −θ 1 with −θ ′ 1 , obtaining another |G(θ 1 )| 2 and we go on doing the same manipulations on −θ 2 . It is clear that at the end we obtain exactly A.8.