Overlap singularity and time evolution in integrable quantum field theory

We study homogeneous quenches in integrable quantum field theory where the initial state contains zero-momentum particles. We demonstrate that the two-particle pair amplitude necessarily has a singularity at the two-particle threshold. Albeit the explicit discussion is carried out for special (integrable) initial states, we argue that the singularity is inevitably present and is a generic feature of homogeneous quenches involving the creation of zero momentum particles. We also identify the singularity in quenches in the Ising model across the quantum critical point, and compute it perturbatively in phase quenches in the quantum sine-Gordon model which are potentially relevant to experiments. We then construct the explicit time dependence of one-point functions using a linked cluster expansion regulated by a finite volume parameter. We find that the secular contribution normally linear in time is modified by a $t\ln t$ term. We additionally encounter a novel type of secular contribution which is shown to be related to parametric resonance. It is an interesting open question to resum the new contributions and to establish their consequences directly observable in experiments or numerical simulations.


Introduction
Understanding the out-of-equilibrium dynamics of isolated quantum many-body systems is one of the most challenging problems in contemporary physics. Due to the direct insight into quantum statistical physics provided by the experimental realisability of closed quantum systems, significant progress has been made both on the experimental and theoretical side in the study of non-equilibrium behaviour. Using cold atomic gases it has become possible to engineer and manipulate isolated quantum systems [1,2,3,4,5,6,7,8,9], and recent studies have led to a series of interesting discoveries such as the experimental observation of the lack of thermalisation in quantum integrable systems [1,2,3,10].
A paradigmatic framework for studying non-equilibrium dynamics is provided by quantum quenches [11] which correspond to a sudden change of some parameters of a system prepared in an equilibrium state, typically in its ground state. For a long time, the focus of the theoretical investigations was the description of the late time asymptotic steady state. To explain the stationary state of integrable quantum systems, the concept of the generalised Gibbs ensemble (GGE) was proposed [12] and later experimentally confirmed [4]. However, specifying the complete set of conserved charges for the GGE in interacting integrable models proved to be a non-trivial problem [13,14,15,16,17,18,19] .
Beyond the steady state it is also of interest to describe the actual time evolution and identify universal features of the non-equilibrium dynamics. Theoretical description of the out-of-equilibrium time evolution is much more difficult and less understood than the characterisation of the steady state, and analytical results have mainly been obtained in systems that can be mapped to free particles [20,21,22,23,24,25,26,27,28,29,30,31,32,33,34] and in conformal field theory [11] and in a few cases in interacting integrable systems [35,36,37,38,39,40,41].
Quantum field theory (QFT) provides an effective and universal description of quantum systems near their critical point. Therefore small quenches in the vicinity of the critical point are expected to be described by a non-equilibrium QFT, capturing universal physics even out of equilibrium. Quenches in quantum field theories, on the other hand, are interesting also in their own right being relevant to high energy physics, and for certain experiments as well [42,43].
In suitably small quantum field theory quenches, the semi-classical approach [44,45,46,47] and approaches based on form factor expansions [48,49,50,51,52] lead to analytical predictions for the time evolution of certain observables. Whereas the perturbative approach in [51,52] can be applied to any quench in which the pre-quench Hamiltonian is integrable, analytical results have only been obtained by perturbation theory up to first order in the quench amplitude. The method developed in [47,48,49,50] can be applied whenever the post-quench Hamiltonian is integrable and the post-quench particle density is suitably small. This latter approach also requires to consider specific, so-called integrable quenches, for which the initial state can be cast in a squeezed coherent form in the post-quench basis as written in terms of the post-quench Faddeev-Zamolodchikov creation operators Z † a (ϑ) for particle species a, the post-quench vacuum |0 and the pair-amplitude K ab , which we will often refer to simply as the K function. It must be stressed that, for models with one particle species as well as for the repulsive regime of the sine-Gordon model, the smallness of the quench essentially guarantees the structure of (1.1), whereas for a generic quench protocol in a generic integrable QFT it is not known and probably not true [52,53] that the initial state obeys this particular structure in the post-quench basis. For interacting integrable field theories, the determination of the K function for a specific quench protocol assumed to be integrable is a difficult and unsolved problem, with some progress in the numerical determination [53,54] and analytic approximations based on form factors [55,56].
Nevertheless, integrable quenches are known to exist for lattice spin systems [57], and are an ideal starting point for analytic considerations in field theory due to their close resemblance to the integrable boundary states introduced by Ghoshal and Zamolodchikov [58]. In fact, following the case of boundary states it is possible to extend (1.1) with zero-momentum particles yielding where g a is called the one-particle coupling which corresponds to the quench breaking (particle number) parity. Such an initial state is also motivated by experiments [43] which show the presence of oscillations with a frequency corresponding to the energy of the zero-momentum particle (i.e. the particle mass).
For an initial state (1.2), time evolution of the vertex operator e iβφ/2 in the attractive regime of the sine-Gordon model was studied in [50]. In the homogeneous (translationally invariant) quenches considered here, the presence of zero-momentum solitons or antisolitons is excluded if the initial state is annihilated by the topological charge which is a typical situation including e.g quenches by local operators starting from the ground state of the model with a different coupling. For stationary breathers B n , however, the one-particle couplings g Bn can appear in (1.2). The conclusion of [50] was that one-particle oscillations of time dependence e −imt show an exponential decay with several relaxation rate contributions. In particular, with only one breather species in the model the contribution of the first breather to the relaxation rate of one-particle oscillations is For boundary field theories, however, it is known that the existence of a one-particle coupling implies a first order pole in the corresponding K function at the origin, the correct relation of the residue to the one-particle coupling was first found in [59] and later proven in [60]. In this paper we show that this relation extends to the case of integrable quenches corresponding to an initial state of the form (1.2) as well. As an immediate consequence, the integral in (1.3) becomes divergent since S B 1 B 1 (0) = −1 . Even though much of the derivation in [50] remains valid, the singular expressions clearly need to be corrected. In this paper we first demonstrate the presence of the singularity (1.4) whenever a zero-momentum particle is present in the initial state. Starting from this observation, we follow a linked cluster expansion similar to the one performed in [48,49,50] but now supported by a finite volume regularisation to treat the singularities and compute the time dependence of one-point functions O(t) up to terms containing five particles focusing on the one-particle oscillations e −imt . As a result of these calculations, we show that the linearly time dependent secular contribution in [50] is modified by a mt ln mt term and a novel type of secular term is encountered which has analogies with the phenomenon of parametric resonance.
The detailed structure of our paper is the following. In Section 2, first a general argument is presented for the existence of the singularity (1.4), and then we consider two particular quench protocols that illustrate our arguments. In the Ising model, quenches from the ferromagnetic to the paramagnetic phase are discussed and the singular structure of the K function is explicitly demonstrated in the continuum limit. Then we consider a quench protocol in the sine-Gordon model which consists of shifting the field by a constant. Computing the quench overlaps using an expansion in the quench magnitude, we demonstrate explicitly that the pair amplitude for the first breather possesses the anticipated singularity up to leading order. In Section 3 we perform a linked cluster expansion for the time-dependent one-point function, using a finite volume regularisation which was first introduced in [61]. This allows us to refine the argument for the existence of the singularity (1.4), and also the explicit construction of the contributions for O(t) up to four particle terms, from which the terms corresponding to one-particle oscillations are extracted. We then consider the fiveparticle terms, but only present the leading order emergent secular contributions. In Section 4 we collect and present our formulas describing the time evolution, address the question of resummation in the linked cluster expansion and generalise our results to cases involving more than one particle species. Finally we discuss a class of secular terms linked to a mechanism analogous to parametric resonance, and conclude in Section 5. Due to the large amount of tedious calculations involved, most of their details are relegated to appendices. In Appendix A, the finite volume formalism is reviewed, while Appendices B and C are devoted to the quenches in the Ising and sine-Gordon models. The technical details of the linked cluster expansion and the calculation of the time evolution can be found in Appendices D, E, G, H and I. To confirm the validity of our calculations they were numerically cross-checked at several points; details of these checks are presented in Appendix J.
2 Integrable quenches with one-particle coupling and the singularity of K In this section we show that in integrable quenches with non-zero one-particle coupling g the K function necessarily possesses a first order pole of the form (1.4) at the origin. For simplicity, we focus on models with one particle species but the argument carries over to systems with several particle species. First we present a general argument based on an analogy with the one-point functions of bulk operators in the presence of boundaries discussed in [62]. The core of the argument is based on the cancellation of divergent terms in the expectation value which is evaluated using finite volume regulators. We then proceed to a concrete example of the quench in the Ising field theory crossing the phase boundary, and discuss an interesting quench in the sine-Gordon model, for which the one-particle coupling and the singular part of the K function can be calculated to lowest order in the quench parameter and the presence of the singularity in K and the relation of its residue to g can be verified explicitly.

Boundary one-point functions
In the following we briefly review the boundary problem discussed in [62]. Let us consider an integrable field theory with a single massive particle constrained on a finite line x ∈ where the coordinate x plays the role of the Euclidean time variable, H is the infinite volume Hamiltonian in the crossed channel and |B is the boundary state corresponding to the boundary condition B. When the boundary state contains zero-momentum particles associated with a nonzero coupling of a single particle state to the boundary in the original channel, |B can be expanded in the asymptotic multi-particle basis [58] as where for simplicity we assumed that the model has only one species, in terms of the boundary reflection factor R B (ϑ), and Z † (ϑ) are the Faddeev-Zamolodchikov creation operators satisfying the commutation relations g B is the one-particle coupling to the boundary and the amplitudes K B (ϑ) satisfy the boundary crossing-unitary equation [58] K which serves as a consistency relation of (2.3) and also implies K B (ϑ) * = K B (−ϑ). The single particle coupling implies a pole for the reflection factor which then yields a pole at the origin for the boundary K function. Whereas it was argued thatḡ B = g B in [58], it was later demonstrated in [59, 63] that the proper normalisation isḡ B = g B /2, with a general proof given in [60]. The relation between the residue of K and the boundary one-particle coupling is crucial for the consistency of a number of theoretical constructs, such as the boundary form factor bootstrap considered in [64]. For the one-point function, the approach of [62] is to put the theory in a finite volume L in the crossed channel (with periodic boundary conditions) and consider the limit where H L is the finite volume Hamiltonian with periodic boundary conditions and |B L represents the boundary state in finite volume. The finite volume is introduced here as a regulator for disconnected contributions arising from the matrix elements of the operator O which appear once |B L is expanded in terms of the corresponding finite volume multi-particle eigenstates of H L . The disconnected terms lead to positive powers of the dimensionless volume variable mL which only cancel if the singularity of K B is exactly as in (2.4). Therefore the singularity of K B and the relation of its residue to the one-particle coupling are consistency criteria for the existence of well-defined one-point functions in the infinite volume theory. The time evolution of the expectation value of a local operator O after a quantum quench with initial state and post-quench Hamiltonian H is given by the expression which is just a real time analogue of the boundary expectation value. The only difference is that the K function appearing in |Ω is not related to any reflection factor; in fact, Ghoshal-Zamolodchikov boundary states are not normalisable, while for quench initial states the factor N is chosen to ensure Ω|Ω = 1. However, it is clear from the calculations performed in [62] that the condition for the cancellation of singularities in unaffected by these details, therefore a well-defined expectation value O(t) after an integrable quench with one-particle coupling g/2, the K function must have a first order pole at the origin with a residue equal to −ig 2 /2.
In the next subsections we discuss two examples where this relation can be verified explicitly, while in Section 3 we consider the real time evolution and show that this condition must hold for consistency.

Quench in the Ising field theory from the ferromagnetic to the paramagnetic phase
In this subsection we study the Ising field theory and quenches across the two phases of the model, as in the Ising model only quenches from the ferromagnetic to the paramagnetic (PM) phase can account for a zero-momentum particle in the initial state. In our approach, we perform the continuum limit of various quantities obtained in the lattice model and in particular, we show that the pair-amplitude possesses a pole with the residue expected from the above considerations.
Quenches in the lattice model were discussed in [25,26,27] and for particular quenches within the ferromagnetic phase, calculations also in the continuum model [48] were performed and numerically checked [65]. Although for the F M → P M quench such calculations in the QFT were not carried out, the scaling limit of the analogous quantities in the lattice model make perfect sense, which we regard as the characteristics of the QFT quench problem. Throughout the subsection, we only discuss the most important features and formulas, and for a more detailed treatment of the topic, Appendix B is recommended to consult.
We first recall, that the transverse field quantum Ising model (TQIM) is defined by the Hamiltonian where σ α i denotes the Pauli matrices at site i, J > 0, h is the transverse field and the boundary conditions are assumed to be periodic. By applying the Jordan-Wigner transformation, the Hamiltonian (2.6) can be mapped to spinless Majorana fermions with dispersion relation [66, 67] and with an energy gap ∆ = 2J|1 − h|. The model possesses a quantum critical point at h = 1 separating the paramagnetic or disordered phase for h > 1 and the ferromagnetic, ordered phase for h < 1. In the disordered phase, the expectation value of σ x i , i.e. that of the magnetisation operator vanishes, while in the ferromagnetic phase its value is non-zero. The Hilbert space of the model consists of two sectors with respect to fermion number parity. In the Neveu-Schwarz and Ramond sectors states with even and odd number of fermions are present, respectively, resulting in the quantisation condition for the wave numbers where n is a positive integer. Performing a quench in the transverse field h, the pre-and post-quench excitations can be related via a Bogoliubov transformation if the initial state is the pre-quench vacuum. As a consequence, the squeezed-coherent form of the initial state in the post quench basis (1.1) is guaranteed. Focusing on quenching from the ground state of the FM phase to the PM phase, ( h 0 → h with h 0 < 1 and h > 1), one can write [25,26,27] 9) where N N S and N R are normalisation constants (2.10) a k and a † k are fermionic operators and K can be found in Appendix B in (B.8). In the scaling limit of the TQIM, J is sent to infinity together with h → 1 such that the gap associated with the fermion mass remains finite M = 2J|1 − h| . (2.11) In addition, the lattice spacing is sent to zero as a = 1 2J . In particular, for the F M → P M quench, we send δh → 0 in the following way: 12) which ensures that the dispersion relation in the post-and pre-quench model is ε h (pa) → M 2 + p 2 and ε h 0 (pa) → M 2 0 + p 2 respectively, i.e. the mass in the PM and FM phase is M and M 0 . Upon the substitution k = pa, the continuum limit of the square of the lattice K function (B.8) in (2.9) is (2.13) From (2.13) it is easily seen that this function has a 1/p 2 singularity at the origin with the coefficient This singularity corresponds to the presence of a zero-momentum particle in the Ramond contribution to the initial state (2.9). As p = M sinh ϑ, the coefficient of the singularity of |K(p)| 2 equals g 4 4 M 2 , therefore K(ϑ) can be written in the form (1.4) with (2.14) Now we show that the one-particle coupling expressed with M and M 0 in (2.9) equals g. To calculate the latter, we put the theory into finite volume, where (cf. also Section 3) the finite volume expansion of the integrable initial state reads where the I denote quantum numbers labelling the finite volume states and N 2 can be found in (A.6). Then from (2.9) must hold. It is convenient to calculate the logarithm of their ratio: when L → ∞, so (2.14) indeed holds.

Phase quenches in the sine-Gordon model
Consider the sine-Gordon model defined by the action in a finite volume L with quasi-periodic boundary conditions The quench protocol consists of abruptly shifting the sine-Gordon field Φ → Φ + δ/β at t = 0, i.e. changing the phase of the cosine potential regarding the pre-quench vacuum as initial state for the post-quench evolution. The peculiarity of this protocol is that, as shown in Appendix C, it is possible to relate the pre-and post-quench ground states by a unitary transformation is the zero mode of the conjugate momentum field Π =Φ. This allows one to derive a form factor expansion for the overlaps with an arbitrary state |χ by expanding the exponential into power series, whereP is the momentum operator and the α i index l − 1 complete sets of asymptotic multiparticle eigenstates. As the initial state is the ground state of a translational invariant Hamiltonian, the overlaps are non-zero only for states |χ of total spatial momentum zero. Due to the integrals over x i , this also restricts the intermediate states |α i L to have zero total momentum so we can write where the tildes over the sums mean that only zero momentum states are included. We now compute the overlap for the first breather to first order in δ/β. Matching the finite volume expression of the initial state (2.19) with the general case (2.15), and using (A.6) and the finite volume form factors (A.14), we obtain where F B 1 is the infinite volume one-breather form factor of Φ (C.11).
Since the form factors of Φ with even number of breathers B 1 vanishes, the lowest non-trivial order for the pair amplitude K is (δ/β) 2 : where ϑ is related to the quantum number I via the Bethe quantisation condition , and θ 1 , ..., θ n are the particle rapidities in the state |α 1 L determined by finite volume quantisation relations analogous to (2.27). Using the expression for finite volume form factors (A.14) it can be written as (2.28) In Appendix C it is shown that in the limits L → ∞ and ϑ → 0 only the n = 1 term survives where the single particle in the intermediate state is also a B 1 , from which Using the form factor kinematical singularity equation (A.12) one can extract that for small ϑ where we used (2.25) which establishes the relation (2.4) between the one-particle coupling of the first breather and the singularity of its pair amplitude for this particular quench in the sine-Gordon model. We remark that for this particular quench protocol we do not know whether it leads to an integrable initial state of the generalised squeezed form (1.2). However, we recall that the proof of the analogous relation in boundary quantum field theory presented in [60] does not depend on integrability either (in fact, it works in general D + 1 space-time dimensions).
Note also that the above argument straightforwardly generalises to a much larger class of quench protocols, the "exponential quenches" when the initial state is related to the post-quench ground state via where Ψ(x) is a local field which breaks particle number parity.

Linked cluster expansion in finite volume
To describe the time evolution of expectation values of local operator, we follow the approach developed in [48,49,50] and apply a linked cluster expansion, combined with the finite volume regularisation scheme used in [62,68]. The latter is based on the finite volume form factor formalism developed in [61]; the ingredients necessary for our calculations are described in Appendix A. For a quench starting from an initial state written in terms of post-quench multi-particle states as in (1.2), a natural approach to compute the one-point function of a local operator O is to decompose into contributions from states with different number of particles, which results in an expansion in terms of form factors of the local operator. However, in infinite volume form factors possess pole singularities due to (A.12), and for quenches with one-particle coupling the K functions also possess singularities. As a result, the contributions are ill-defined and need to be regularised which can be done by putting the theory in finite volume, where due to the quantisation of the particle momenta neither the kinematic singularities of the form factors nor the singularities of the pair-amplitude contribute. The finite volume L can be considered as a physical regulator and the expectation value is then written as which is first evaluated for finite L where one can verify the cancellation of singular terms explicitly and then take the limit L → ∞.
To perform the calculation, one needs an expression for the initial state in finite volume which was derived in [62]: where I, J are Bethe quantum numbers and we used the notations introduced in Appendix A. To simplify notations, one can write the following shorthand for where |Ω (0) = |0 L , |Ω (1) = g 2 N 1 (L)|{0} L , etc. denote the contributions with a fixed number of particles.
To ensure the convergence of the expansion for high energies one can introduce a regulator parameter R and consider where R > 0. When recasting the quantum number sums in terms of contour integrals, R ensures that the integrals themselves are convergent and therefore allows appropriate manipulations of the contours. At the end of the calculation the parameter R is sent to zero.
Following the procedure introduced in [69], one can separate contributions indexed by particle number as follows and for a proper normalisation of the state one must also divide by the "partition function" In particular for Z, the first few terms are and Let us turn to the issue of the expansion parameter. Whereas in our calculations R is eventually sent to zero at the end, for O(t, R) is expected to be well-defined for any finite R. Therefore let us first treat R as a large positive quantity, and introduce the parameters z = e −m(R/2+it) and z = e −m(R/2−it) . Then the order of C kl isz k z l , and that of Z n is (zz) n . The inverse of the partition function, thus can be expanded in powers of zz as where the first few terms read Putting these ingredients together, in a finite volume L we obtain where analogously to to [62] and [69],D nm is introduced as with the first few terms having the form Since theD kl are of order z kzl , they must separately be well-defined as L → ∞: and so the infinite volume limit can be written as For convenience and later use, we also introduce the quantities G n = n l=0D n−l,l , (3.9) whose infinite volume limit is denoted by G n .
The expressions individual C kl contain finite volume form factors, which in general are given by where it is understood that the rapidities {ϑ 1 , . . . , ϑ k } and {ϑ ′ 1 , . . . , ϑ ′ l } are solutions to the corresponding Bethe-Yang equations with quantum numbers {I n }, {J n }. Formula (3.10) is valid whenever there are no coinciding rapidities; otherwise a more complicated formula taking into account disconnected contributions is necessary. In this paper we are interested in contributions to one-particle oscillations, for which coinciding rapidities cannot occur, the numbers of particles in the two multiparticle states differ by an odd number which excludes the two possible cases with disconnected terms (cf. [61]).
Note that the equality (3.10) is valid up to a suitably chosen phase factor which can be changed by redefining the phases of the finite volume eigenstates |{I 1 , . . . , I n } L . This includes also the fact that the ordering of the particles is not determined by first principles and any exchange leads to an Smatrix factor according to (A.10). It is clear that all such ambiguities cancel in the expectation value (3.3); however, for a practical calculation one must fix the phases of the multi-particle contributions to the matrix elements consistently. Here we make use of the consistent prescription introduced in [62]: any time the amplitudes K(ϑ i ) and K * (ϑ i ) appear with some ϑ i , the explicit order of the rapidities substituted into the relevant form factor is given by (−ϑ i , ϑ i ) and (ϑ i + iπ, −ϑ i + iπ), respectively. Exchanging any two pairs of rapidities does not make any difference, therefore the phase of the form factors is completely fixed by the above rule. Note that the presence of zero-momentum particles does not produce any additional ambiguities.

The singularity of K revisited
Now we are ready to complete our arguments why the pair-amplitude must be singular in integrable quenches with one-particle coupling. Recall that in the boundary problem [62] the ordering of the terms was performed according to powers of e −m(R−x) and e −mx , which resulted in the same expressions for the boundary version ofD kl that we obtained in (3.7) with the expansion parameters z = e −m(R/2+it) andz = e −m(R/2−it) . In the boundary problem, the existence of the infinite volume limit (3.8) and eventually O(x) B requires the presence of the singularity g 2 ϑ must hold for the quench problem as well. The easiest way to see that is to consider the one-particle contribution −Z 1 which behaves as mL. To make C 12 − Z 1 C 01 finite in the infinite volume limit, C 12 must have a similar volume-dependence, which is ensured by the singularity of K involved in C 12 .
At the end of the calculations the regulator R is sent to zero. The contribution G n (3.9) is of order |z| n = (e −mR/2 ) n and just as in [62] it turns out that the coefficient of the largest power of the mL term is always of order g n . Focussing on the behaviour of the singular term, the small parameter of the linked-cluster calculation can be identified with g. As the singularity of the K is of order g 2 one can formally treat K as a term of order g 2 . This way of counting the orders results in the same classification of contributionsD kl (3.9) that we obtained by considering z andz as the expansion parameters.
This counting of orders is clearest for a perturbative quench corresponding to changing the Hamiltonian as where Ψ is a purely odd operator (i.e. whose form factors with an even number of particles in the pre-quench system is zero). Using perturbation theory to compute the overlaps following [51,52] one obtains that the one-particle coupling is of order λ, while the pair amplitude K is always of order λ 2 .
When the perturbing operator has even matrix elements as well, the pair-amplitude can also possess a λ order term. Whereas this term is not singular due to the regular behaviour of F (iπ + ϑ, iπ − ϑ) at ϑ = 0, at order λ 2 a singular term similar to the one found in Section 2 is always present. Therefore the presence of the zero-rapidity pole singularity of K is generic.

Contributions up to 4 th order: analytic continuation of the boundary result
In this section we present all the terms up to fourth order using the Euclidean quantities computed in [62] and continuing them to real time. The Euclidean one-point function was computed up to contributions D kl with k + l ≤ 4 which are collected in Appendix E. For the analytic continuation, we apply the R → 0 and x → −it substitutions together with K B → K and g B → g which give G 0 : 0|O|0 , (3.12a) Note that these integrals remain well-defined even when there is a pole in the amplitude K(ϑ) at ϑ = 0 because the form factors possess a zero ϑ i = 0 as a consequence of the exchange axiom (A.10) and the general property S(0) = −1. Concerning the large rapidity behaviour, normalisability of the initial state requires that K tends to zero fast enough for large ϑ, again ensuring the existence of the integrals.
We now turn to analysing the time dependence originating from (3.12). Due to the oscillatory integrands it is convenient to apply the stationary phase approximation (SPA) briefly discussed in Appendix D.4. Both SPA and direct analysis leads to the following type of terms expected from (3.12): where α is either integer or half integer and n is an integer. For terms D kl the lower bound of the oscillation frequency is always nm = (k − l)m.
As the main objective of this paper is to study the time dependence of one-particle oscillations, we concentrate here on terms with n = 1 but we will also briefly comment on the time dependence of the non-oscillatory part of O . The non-oscillatory parts include the static contributions G 0 : 0|O|0 , whereas for the only time dependent integral, the SPA (D.16) can be applied, yielding where C = lim which is finite. For terms with one-particle oscillation, one finds (3.18) For the last term, SPA cannot be applied directly, therefore we shift the contour off the real axis where (as shown in [62]) the contribution from the term cosh ϑ sinh 2 ϑ vanishes and reintroduce the regulator R . We rewrite the resulting expression using (D.14) after which the SPA (D.16) can be applied, and finally perform the R → 0 limit. The result is For a more detailed derivation, the interested reader is referred to in Appendix F. While the first term has the standard ∼ 1/ √ t time dependence, the second one behaves as ∼ √ t for large time. We return to this peculiar finding in Sec. 4.3.

Leading order time dependence from G 5
The contributions involving five particles are D 05 , D 14 , D 23 and their complex conjugates D 50 ,D 41 , D 32 . Based on (3.7), the expressions to evaluate are (3.20) The calculation of these expressions is very tedious; even in Appendices G, H and I where the details of the calculations are presented, we focused only on leading secular contributions to the one-particle oscillations, i.e. terms of the form e −imt t α with the highest power α. As one-particle oscillations originate exclusively from D kl with |k − l| = 1, we can focus on D 23 and its conjugate D 32 . The leading order secular terms have two origins: a residue contribution from encircling the poles of the form factors F 5 (iπ + ϑ 1 , iπ − ϑ 1 , −ϑ 2 , ϑ 2 , 0) when ϑ 1 ≈ ϑ 2 and a contribution from these poles when a contour integral is performed with an integration contour just above the real axis.
Concerning the former term, the explicit expression reads D Res whose derivation can be found in Appendix G. Unlike time dependent terms discussed so far, deriving (3.21) involves no SPA, hence it is valid also for small times. Note that the coefficient of the oscillatory factor e −imt is purely imaginary and linear in time.
For the other term denoted by D CInt 23 (t), the explicit formula reads where and with Their derivation can be found in Appendices H and I. It is shown in Appendix I that in the large time limit mt ≫ 1 the integral kernels behave as and where F S and F C are the Fresnel sine and cosine functions, respectively. This leads to the final result (3.32)

Discussion of the results
Before discussing the time evolution of the one-point function, we collect the leading order timedependent contributions for the non-oscillatory and one-particle-oscillatory part of O from each term D kl we considered in the previous section. In the long-time limit 1 ≪ mt, these are and where C is given in (3.17) and γ 1,2 in (3.31).
In the preceding subsections we found that the long time asymptotics of the leading order contributions to the one-particle oscillations contain, besides the original oscillation e −imt from G 1 , two new types of terms: one with time dependence √ te −imt from G 3 and terms of time dependence te −imt and t ln te −imt from G 5 .
Since these are secular terms growing for large t, it is necessary to sum up higher order contributions coming from G 2n+1 . Computing terms G 7 and higher is extremely tedious and has not yet been performed, therefore in this work we can only present a limited discussion of their resummation based on insights gained from earlier works [48,49,50].
In addition to these works, one could also try to compare with the F M → P M quench in the Ising model considered in Subsection (2.2); however, we show that this is unfortunately not possible. In Appendix B it is shown that in the continuum limit the time evolution of the magnetisation is where the relaxation time is given by for M < M 0 . At first sight, since the K function is known for this particular quench, the expansion of (4.3) can be matched with the form factor expansion evaluated in our work. However, note that (4.3) and (4.4) are non-analytic functions of the pre-quench mass M 0 around the origin, and therefore are also non-analytic in g around g = 0 due to the relation (2.14). This is not surprising since a quench across a quantum critical point is expected to be a large quench with possibly non-analytic behaviour, and therefore the form factor expansion is not expected to be valid at all. This situation is in marked contrast with the phase quenches in the sine-Gordon model considered in Subsection (2.3) where the shift δ/β can play the role of a small parameter.

Connection with previous results and discussion of possible resummation of G 4n+1
First we consider terms G 4n+1 . Contributions to one-particle oscillations originate from the pieces D 2n,2n+1 . In [50] these terms served as the sole source of secular terms which were shown to sum up to an exponential function to order t 2 , i.e. the result including the two leading order corrections had the form which is the expansion of The real part of the above integral is the relaxation time, while the imaginary part is a frequency shift. It is easy to see that sending g → 0 in our expressions (3.30), which is equivalent to removing the singularity of K, reproduces the result for τ above as the integrand is non-singular and The origin of (4.6) is essentially the kinematic singularity of the form factors. In our calculation we have an additional ingredient, namely the singularity of K which gives some new contributions according to (3.30). Assuming an exponentiation similar to that observed in [48,49,50], the leading order time dependence from G 4n+1 is which can be resummed into Therefore, besides the frequency shift, the relaxation for late time is naively expected to be superexponential with a dependence of the form e −t ln t . This also means, that under the assumption of exponentialisation, τ −1 in (4.6) is to be replaced with m g 4 4 log(mt) π − γ 1 − iγ 2 in the exponential function. However, this cannot be concluded safely without computing at least D 45 (t) and checking whether one obtains the correct combinatorial coefficients for the terms involving higher powers of − g 4 4 log(mt) π mt which is the condition for exponentialisation as established in [48,49,50]. Based on analogy with [50], one can argue that terms containing γ 1,2 exponentialise and their resummation leads to a relaxation rate and a frequency shift; however, the fate of the logarithmic term cannot be decided without examining higher order contributions, whose straightforward evaluation is extremely complicated. Assuming that the relaxation is of the usual exponentially decaying form (i.e. the logarithmic part does not exponentialise) also means that it is not clear at the moment what part of the real terms (involving γ 1 ) exponentialises and determines the relaxation rate.

Multiple species
It is easy to generalise our results to the case of multiple species following the reasoning in [50] which studied relaxation in the attractive regime of the sine-Gordon model for operators semi-local to soliton excitations (and consequently local with respect to the soliton-antisoliton bound states, i.e. the breathers B n ). In the following we write down the result for semi-local operators where and where j, k, l, o index breather excitations, and γ (n) 1,2 are obtained by substituting S in (3.31) with the B n − B n scattering amplitude S nn and g with the B n one-particle coupling g n . The expressions for γ (llk) 1 and γ (llk) 2 are obtained from (3.31), by replacing g with g j and also in (3.23) and (3.24); in addition, in the residue term D Res 23 (3.21) it is necessary to replace (4.14) Since solitons can have no one-particle coupling in the quench due to topological charge conservation, the amplitude K ss is regular at the origin. Note that (4.11) are identical to the results in [50]. It is easy to understand this from the fact the expressions given in [50] for τ −1 jk with j = k and for τ −1 jkl with j = l, k = l, are regular even if the K j are singular. The only terms to be revisited are τ −1 jj and τ −1 lkl where the naive application of (4.11) results in divergent contributions and must be modified using our calculations performed above.
Results for theories with multiple species with fully diagonal scattering can be obtained by omitting the solitonic contributions from the above results and replacing the breathers with the actual particle spectrum.

Parametric resonance
Finally, we turn to contribution (3.19) The origin of this term is the integral which is derived from a finite volume contribution of the form where ϑ is determined by the quantisation rule Let us recall the phenomenon of parametric resonance, for which the simplest example is the Mathieu equation describing a system with only one degree of freedom, This equation has a region of instability in which the solution of the equation oscillates with an exponentially growing amplitude. This region is when the ratio ω p /ω 0 is sufficiently close to 2, where the width of the region of proportional to q. It is clear that the term (4.17) couples the one-particle mode with the two-particle modes and satisfies the condition for resonance on the threshold; however, the interplay between the integration over rapidity and the singularity of the form factor F O 3 (iπ, −ϑ, ϑ) results in a growth of √ t instead of being exponential. This cannot be the final story: since a quench pumps a finite energy density in the system, the oscillations cannot grow without bound and therefore higher terms G 4n+1 must modify this behaviour to keep the amplitude bounded. Since at present we do not have control over these higher order terms, we cannot predict the eventual fate of this class of contributions.
A heuristic analogy can be drawn by noting that in the Mathieu equation the driving oscillator is external, while in the quench the two-particle modes are also dynamical and the total energy stored in the system is conserved. A closer analogy for this dynamics is provided by the differential equation describing coupled modes The mode x is parity odd and is the analogue of the one-particle mode, while the mode y is parity even and is the analogue of the two-particle mode. Depending on the choice of the parameters, this system has solutions in which the energy stored in the modes shows an oscillatory behaviour in time.
With the choice ω 2 = 2ω 1 the condition of parametric resonance is satisfied, but in contrast to the Mathieu equation, the system has a total energy which is conserved since the driving and driven modes now form a closed system. An example of a resonance solution is shown in Fig. 1; notice the long plateau in the amplitude of the driving mode y showing the non-linearity of this system. This is in contrast to linearly coupled oscillators where the energy transfer would itself have a harmonic dependence on time.
It is an interesting question whether in the full quench situation such a non-trivial behaviour can be observed in the time evolution for some choice of the parameters. Since the √ t term is of order g 2 and the next secular contribution is of order g 4 , it is quite possible that for small enough quenches there is some intermediate time window in which the parametric resonance dominates. To study this requires the detailed analysis and resummation of higher order terms which is left for further investigation.

Conclusions
In this work we studied integrable quenches with zero-momentum one-particle states in the initial state and the subsequent time evolution of one-point functions. We found that (similarly to the case of integrable boundary conditions) in the presence of a non-vanishing one-particle coupling the two-particle amplitude K must have a singularity at the origin ϑ = 0, where g is the one-particle coupling determining the overlap of the post-quench zero-momentum one-particle state with the initial state. We gave a general argument based on the relation of the quench time evolution to a boundary field theory problem discussed in [62]. In addition, we presented two explicit examples, the quenches from the paramagnetic to the ferromagnetic phase in the Ising QFT and phase quenches in the sine-Gordon model, in which the relation between the one-particle coupling and the residue of K at zero rapidity could be established using explicit calculations and perturbation theory. We also pointed out that the sine-Gordon phase quench is a particular member of a more general class of "exponential quenches", where the same argument applies. A general proof establishing the relation (5.1) also follows from our explicit computation of the time evolution using finite volume as a regulator, which shows that the infinite volume limit only exists when (5.1) holds. We stress that for the results in [50] the presence of a pole in K is problematic as the expressions for the relaxation times are not well-defined due to divergence of the integrals involved in the formulas.
We then proceeded with the explicit computation of the time evolution of one-point functions focusing on the one-particle oscillations with time dependence e −imt . We used a modification of the linked cluster expansion introduced in [48,49,50], where instead of rapidity space point splitting, we applied a finite volume regulator first proposed in [61]. The advantage of this regulator is that it uses a physical parameter, i.e. the system size L, and since the thermodynamic limit must be well-defined, the computed one-point function must have a finite limit as L → ∞. The cancellation of terms containing positive powers of the volume (resulting from kinematical singularities) provides an important consistency check for the computation.
We found two important secular contributions to the one-particle oscillations. The first one takes the following form at leading order The terms γ 1 and γ 2 are directly analogous to the contributions found in [50], with two essential differences. First, the integrals expressing the coefficients γ 1,2 are now entirely well-defined. If there was no pole in K(ϑ) (which we argued to be impossible for g = 0), these terms would reduce to the expressions given in [50]. The second difference is the presence of the logarithmic timedependence for g = 0. In analogy to [50] it is expected that either higher order terms G 4n+1 get resummed to an exponential function e −t/τ with τ −1 → m g 4 4 log(mt) π − γ 1 − iγ 2 , or alternatively, only terms containing γ 1,2 are resummed into a relaxation rate and a frequency shift of the oneparticle oscillations. We stress that the fate of the logarithmic term is unclear at present, which also implies that the part of γ 1 entering the resummation is not defined until this issue is dealt with, which requires further investigation.
Our results can be easily extended to local operators in theories with more than one species with diagonal scattering, and also to operators semi-local with respect to solitons in the sine-Gordon model.
The other class of secular contributions to one-particle oscillations is a novel one having the form and its leading order is g 3 . The origin of this secular term is a physical effect analogous to parametric resonance, and is caused by the effective coupling between one-and two-particle modes. This coupling is established by the corresponding form factor of the operator and the K function resulting in a singularity at the threshold of the two-particle continuum. At this point the ratio between the frequencies is exactly two satisfying the condition of parametric resonance. Note that the singularity of the K function is an essential ingredient as it is the origin of enhancement in the effective coupling between the modes.
Unfortunately, it is rather difficult to calculate the next secular contribution from G 7 and even having an expression for the next term, it is not guaranteed that they can be resummed in an effective manner to extract the long time behaviour. Therefore it is presently unclear how this phenomenon influences the fate of the one-particle oscillations. However, it is clear that despite the initial growth indicated by the presence of √ t , the amplitude must eventually saturate as only finite energy density is injected in the system during the quench.
In most of our discussions we tacitly assumed that the quenches are integrable, so that the initial state has the generalised squeezed state form (3.21). However, we wish to note that the arguments establishing the singularity of the pair amplitude are valid even for non-integrable quenches whenever a one-particle coupling is present. A generic initial state can always be written as [56] Integrable initial states of the form (5.1) satisfy K n = 0 for n > 2, but the main conclusions of this paper are unaffected by the presence of the higher K n , which would only result in additional higher order corrections to the frequency shift and relaxation terms. Finally, we wish to comment on the relevance of the phase quenches in the sine-Gordon model introduced in Subsection 2.3. Being both experimentally realisable and analytically tractable, the sine-Gordon theory has attracted a lot of attention. The sine-Gordon model emerges as an effective description in cold atom experiments involving tunnel coupled quasi-one-dimensional condensates, where the sine-Gordon field corresponds to the relative phase of the condensate [42,43]. Therefore, phase quenches can be an ideal protocol to compare experimental and theoretical results. This is especially true for moderate sized quenches, where the small post-quench density of excitations makes the form factor series valid and it is possible to extract analytic results about the post-quench time evolution.

A Finite volume formalism
A.1 Multi-particle states in finite volume Excited states of a massive integrable quantum field theory with one particle species in a large, but finite volume can be described as scattering states consisting of n particles with rapidities ϑ n given by the solution of the Bethe-Yang equations Using the fact that the effective statistics is fermionic, i.e. S(0) = −1, the two-particle phase shift function is defined as S(ϑ) = −e iδ(ϑ) , δ(−ϑ) = −δ(ϑ) , and the following prescription is obtained for the quantum numbers of the particles: I k ∈ Z for odd n , I k ∈ Z + 1 2 for even n .
The state corresponding to quantum numbers {I 1 , . . . , I n } is denoted as and it is independent (up to a possible phase ambiguity) of the ordering of the I k . They are normalised so that their scalar products are . . δ In,I ′ n , with the quantum numbers ordered by convention as I 1 < · · · < I n and I ′ 1 < · · · < I ′ m . The total energy and momentum can be expressed as up to exponential corrections governed by some mass scale µ. For a systematic treatment of exponential corrections to excitation energies see [70,71] and also [72,73]. The rapidity space density of states in the n-particle sector of the theory is given by the Jacobian

A.2 The initial state in finite volume
The integrable initial state (1.1) can be written in finite volume as [62] where the rapidities ϑ are the solutions of the appropriate Bethe-Yang equations with a constraint of zero overall momentum. For the two-particle states the constrained Bethe-Yang equation is whereas for the four-particle case the quantisation condition for the rapidities {−ϑ 1 , ϑ 1 , −ϑ 2 , ϑ 2 } is given by the system of equations The normalisation factors N 1 (L), N 2 (θ, L), N 3 (ϑ, L) and N 4 (ϑ 1 , ϑ 2 , L) in (A.3) are not determined by first principles and were calculated in [62] up to finite size effects exponential decaying with the volume. For the one and two-particle states one finds Note that the total two-particle density ρ 2 satisfies and so For the three-particle state the normalisation of the three-particle states is given by and in the four-particle case the normalisation reads

A.3 Form factors and their properties
In integrable models, the elementary form factors are meromorphic functions of the rapidities and satisfy a series of form factor bootstrap equations (for a review we refer to [74,75]). In a theory with only one particle species and in the absence of bound states the form factor equations are: I. Lorentz invariance: where s is the Lorentz spin of the operator O. In this work we only consider scalar operators corresponding to s = 0. II. Exchange: III. Cyclic permutation: where σ is the mutual locality index of the operator O with respect to the interpolating field φ that creates the particle excitation and is defined via the condition O(x, t)φ(y, t) = e 2πiσ φ(y, t)O(x, t) .  For the case of coinciding rapidities this relation must be modified to account for disconnected contributions [61], but we do not need the corresponding expressions here.
B The transverse field quantum Ising chain and its scaling limit

B.1 Definition of the model
The transverse field quantum Ising model (TQIM) is defined by the Hamiltonian where σ α i denotes the Pauli matrices at site i, J > 0, h is the transverse field and the boundary conditions are assumed to be periodic. By applying the Jordan-Wigner transformation, the Hamiltonian (B.1) can be mapped to spinless Majorana fermions with dispersion relation [66, 67] having an energy gap ∆ = 2J|1 − h|. The model possesses a quantum critical point at h = 1 separating the disordered or paramagnetic phase (PM) for h > 1 and the ordered, ferromagnetic phase (FM) for h < 1. In the disordered phase, the expectation value of the magnetisation operator vanishes, while in the ferromagnetic phase its value is given by The Hilbert space of the model consists of two sectors with respect to fermion number parity. In the Neveu-Schwarz and Ramond sectors states with even and odd number of fermions are present, respectively, resulting in the quantisation condition for the wave numbers where n is a positive integer. In particular, the Fock space of the model in the paramagnetic phase can be written as 5) and the ground state is the Neveu-Schwarz vacuum, |0 P M N S . In the ferromagnetic phase, the zero momentum excitation has negative energy, therefore the Ramond vacuum is redefined as |0 R → a † 0 |0 R after which a particle-hole transformation is implemented on the zero-mode a † 0 → a 0 . The Fock space of the model in the FM regime is hence where the ground state becomes degenerate in the thermodynamic limit and the finite volume states corresponding to macroscopic magnetisation are 1

B.2 F M → P M quench and σ x i (t)
Performing a quench in the transverse field h, the pre-and post-quench excitations can be related via a Bogoliubov transformation if the initial state is the pre-quench vacuum. As a consequence, the squeezed-coherent form of the initial state in the post quench basis (1.1) is guaranteed. Focusing on quenching from the ground state of the FM phase to the PM phase, ( h 0 → h with h 0 < 1 and h > 1), one can write [25,26,27] |0, h 0 where N N S and N R are normalisation constants and The late time evolution of the magnetisation operator is given by the expression [26]: where k 0 is a solution of the equation cos ∆ k = 0, α(h, h 0 ) is an unknown constant and (B.11)

B.3 Continuum limit of the model and the F M → P M quench
In the scaling limit of the TQIM, J is sent to infinity together with h → 1 such that the gap associated with the fermion mass remains finite In addition, the lattice spacing is sent to zero as a = c 2J , where c is the speed of light that we set to 1. It is easy to see that the dispersion relation (B.2) under scaling limit transforms as with {ψ(x, t),ψ(y, t)} = 2πδ(x − y). The lattice magnetisation operator σ x i is related to the continuum field σ with the conformal normalisation via σ(na) =sJ where A=1.282427129... is Glaisher's constant.
In the following we perform the continuum limit for quantities relevant for the F M → P M quench. It is important to note, however, that unlike for quenches within the ferromagnetic phase, where explicit calculations in the field theory framework [48] and numerical studies [65] were carried out, for the F M → P M quench no such investigations have been done. On the other hand, performing the naive scaling limit, the resulted QFT quantities make perfect sense, hence expected to be the results for the field theory quench problem. In the scaling limit sending δh → 0, we write Note that ζ has a finite limit − ln 2 when M 0 = 0, but its derivative is infinite at the origin as M 0 → 0, hence ζ is not an analytic expression for small values of M 0 .

B.4 The singularity of the Ising K function
From (B.18) it is easily seen that this function has a 1/p 2 singularity at the origin with the coefficient The appearance of this singularity is consistent with our consideration presented in paper, as from the expansion of the initial state (B.7) it is seen that there is a zero-momentum particle in the Ramond sector. In the following, we show that the pole strength in the K or more precisely in |K| 2 is g 4 /4, if the one particle coupling is g. For this goal, we first determine g from the pole as

(B.24)
It is convenient to calculate the logarithm of the ratio, which we write as When M L → 0, the expression denoted by B in (B.25) is zero, which can be easily seen by considering the Euler-Maclaurin formula: for the two terms in B the difference between the lower endpoints for the integration and for the boundary terms consisting of the functions and their derivatives goes as 1/L and the expressions in the logarithms are smooth, well-behaved functions. The term A in (B.25) can explicitly be computed, therefore when L → ∞,

C Phase quenches in the sine-Gordon model and exponential quenches
In this appendix we give more details of the calculation presented in Section 2.3. Recall that the phase quenches in the sine-Gordon model consist of abruptly shifting the sine-Gordon field φ → φ + δ/β, i.e. changing the phase in the model and regarding the pre-quench vacuum as the initial state for the post-quench evolution. To obtain the relation between the pre-quench and post-quench ground states, consider the sine-Gordon model as the perturbation of the compactified free massless bosonic conformal field theory in finite volume with the Hamiltonian where the V a are the so-called vertex operators (normal ordered exponentials of the boson field). Note that the free bosonic part of (C.1) commutes with the zero mode of the conjugate momentum field Π 0 = ∂ t Φ, whereas due to the canonical commutation relations Therefore the initial state in the finite volume theory with PBC reads For a form factor expansion for the overlaps in finite volume, one can use the expression Expanding the exponential into a Taylor series, the general term is in whichP is the momentum operator and α i indexes a complete set of asymptotic eigenstates. Due to translation invariance, the overlaps are only non-zero when the total momentum p of the state |χ is zero, and so only intermediate states with zero total momentum contribute. We can then write where the tildes over the sums signal that only zero momentum states are included. Now for a state containing a zero-momentum first breather and to first order in δ, we can use (C.8), (A.6) and (A.14) to obtain which gives where F denotes the form factor of the field Φ in the sine-Gordon theory which has only has nonvanishing matrix elements with states composed first breathers when their total number is odd. In particular [76], (C. 12) and we remark that the form factors of the sine-Gordon field φ can be obtained from that of the vertex operators : e iαφ : by differentiating with respect to α. For the pair overlap the lowest nontrivial order is δ 2 and one obtains where the particles A k are either breathers or solitons and their is an implicit summation over all possible choice of their species. Now we are interested in the singular behaviour of K which can only originate from those of the form factors. The infinite volume limit itself is finite, since the numerator contains an explicit L 2 , and in the n-particle term the denomination contributes 1/L and also 1/L n from the density factor, while rewriting the discrete summation in terms of integrals results in a state density factor of behaviour L n−1 , so the leading term is independent of L.
First let us focus on terms where all particles A 1 . . . A n are first breathers. We now proceed to analyse the singularity of these terms by setting ϑ to a value ǫ and to consider the limit of small ǫ. The form factor has a kinematical singularity when some subset of m particles among the A 1 . . . A n has similarly small rapidities, which we take to be a multiple of ǫ. The dependence of the most singular term can be obtained from the kinematical residue equation (A.12) but also taking into account that the form factor having a first order zero when two rapidities coincide due to (A.10) and S(0) = −1. There are three cases: • m < n − 1: The behaviour of F (iπ + ϑ, iπ − ϑ, θ 1 , ...θ n ) is ǫ 1+(m−1)m/2−2m , where the 1 comes from the coincidence of iπ − ϑ and iπ + ϑ , the m(m − 1)/2 comes from the coincidence of the m rapidities in the set θ 1 , ...θ n , and the −2m are from the kinematical singularities "activated" in the limit. For F * (θ 1 , ...θ n ) one obtains a factor ǫ (m−1)m/2 , therefore the overall behaviour is ǫ 1+(m−1)m−2m which is only singular for m = 1 or 2.
• m = n − 1: one must be aware that the condition of zero total momentum, that sending m = n − 1 rapidities to zero results in all the n of them going to zero, leading to the behaviour ǫ 1+(n−1)n−2n which is singular for n = 2.
Note that whenever there is a singularity it is always of order one. For the case when all the A 1 . . . A n are first breathers, n must also be odd for the form factor not to vanish due to parity invariance. The n = 1 contribution is (C.14) Using the form factor equations 15) this contribution alone gives the expected pole contribution ϑ .
(C. 16) Now we demonstrate that contributions with n ≥ 3 are regular at the origin ϑ = 0. Observe that in finite volume the rapidity ϑ in finite volume is eventually quantised according to with some half-integer quantum number I, so it is always displaced by an amount of order 1/L from the origin. Fixing I results in the parameter ǫ being essentially 1/L, therefore the singularity 1/ǫ manifests itself as a divergence of (C.13) proportional to L when L goes to infinity with I fixed. A simple power counting then gives that the contribution is where L r is the state density factor resulting for the r particle rapidities among the θ 1 , ...θ n whose summation is left free once fixing the positions of those needed for the singularity and also taking into account the zero-momentum constraint. Clearly one obtains r ≤ n−2 resulting in a cancellation of any divergence for L → ∞.
To finish this discussion, let us consider the case when the set of particles A 1 . . . A n contains other species (higher breathers or solitons) as well. Let us suppose that the total number of first breathers among A 1 . . . A n is k with k < n. The counting of the degree of singularity only involves first breathers, so we obtain that the only possible cases are again k = 1 or 2. However, in the denominator of (C.13) now one has a density ρ n with a behaviour L n with n > k, which leads to a regular limit for L → ∞.

D Some useful relations
Here we collect some useful formulae regarding the K function and form factors, identities from the theory of distributions and relations for the stationary phase approximation (SPA) that are useful in the text.

D.1 The K function
The K function possesses a singularity at the origin, if the one particle coupling to the boundary denoted by g is non-zero. In this case, the singular term in K is Due to the relation K(−ϑ) = S(−2ϑ)K(ϑ) the constant term in the expansion of K is also expressible with g. Writing K as and expanding Also note that due to real analyticity K(−ϑ) = K * (ϑ), all even/odd coefficients in the expansion of K around ϑ = 0 are purely real/imaginary, respectively.

D.2 Form factor singularities and expansions
Consider form factors for a single species, such as the the multi-B 1 form factors of the sine-Gordon model. Using the form factor equations (A.10) and (A.12) one can derive the universal expression i.e. the first derivative of the S-matrix at the origin. Another form factor we need is F 5 (iπ, −ε, ε, −ϑ, ϑ), where ϑ is non-zero. Using (A.13), one can explicitly compute which is universal as well. For F 5 (iπ + ϑ + ε, iπ − ϑ − ε, −ϑ, ϑ, 0) the most singular term is O(ǫ −2 ), however, we also need the ε −1 terms. Denoting the coefficient of the sub-leading singularity by F ε 5 (ϑ), it can be shown based on (A.13) that

D.3 Some distribution identities
Suppose that f (x) is a well behaved function with vanishing at infinity. Then it is well-known that where P´denotes the principal value. This identity has the following counterpart for second order singularityˆ∞ One can also writeˆ∞ hold as well. A useful way of evaluating the principal value integral is where g(x) is an appropriate mask function that is even, grows at the infinity and satisfies g(0) = 1, a convenient choice being g(x) = cosh x.

D.4 Stationary phase approximation
Consider the following integral with oscillatory argument: in which g(x) has one global extremum at x 0 , and f (x) is regular and decays fast enough for large |x|. The asymptotic behaviour of this integral for large t can then be evaluated as and

F Evaluating D 12
In order to analyse the time dependence of the term we reintroduce the regulator R enabling us to shift the contour off the real axis where (as shown in [62]), the contribution from the term cosh ϑ sinh 2 ϑ vanishes, i.e., we end up with Using now (D.14), we find where we placed the integration contour back to the real axis as the integrands are now free of poles for real rapidities and also got rid off the regulator by setting its value to zero. Performing the SPA (D.16) results in While the first term has the standard ∼ 1/ √ t time dependence, the second one behaves as ∼ √ t for large time. This peculiar finding is discussed in Sec. 4.3.
G Evaluating G 5 , part I. Notations, D 05 , D 14 and residue terms from D 23 In this and the following appendix we describe how to compute the five-particle contribution G 5 . According to the discussion in Section 3 it can be obtained as (the real part of) the infinite volume limit ofD 05 +D 14 +D 23 , (G.1) The calculation ofD 05 andD 14 is easy; we show explicitly the cancellation of terms with positive powers of mL and extract the time dependence based on stationary phase approximation (SPA). It turns out that this time-dependence is only of sub-leading order. The calculation ofD 23 is, however, so long and tedious that we only extract terms accounting for the leading time-dependence and then carefully argue why other terms produce only sub-leading effects. In particular, in this Appendix we only focus on residue terms and leave contributions resulting from a contour integration for Appendix H. Finally, we provide numerical evidence for the existence of the infinite volume limit demonstrating the cancellation of terms behaving as mL and (mL) 2 in Appendix J. We use the following notation to abbreviate formulas involving form form factors: whenever a vertical line is seen in the argument of the form factors, rapidities in the argument refer to rapidity pairs, and a single rapidity is indicated by putting it into brackets {}. Rapidities to the right of the line are understood to be shifted by iπ. For example, Let us start with

In a similar spirit, we introduce
which is free of divergences, therefore its infinite volume and R → 0 limit is simply Concerning the long time asymptotics of this term, the stationary points are ϑ 1 = ϑ 2 = 0 where the product of the form factor and the K factors can be expanded in non-negative powers ϑ n 1 ϑ m 2 . Applying a Gaussian approximation´d ϑ 1 2π´d ϑ 2 2π ϑ n 1 ϑ m 2 e −itm(1+2ϑ 2 1 +2ϑ 2 2 ) yields time dependence of the form t −1−m−n which we neglect.
Turning now to the corresponding form factor F 5 (iπ, −ϑ 1 , ϑ 1 , −ϑ 2 , ϑ 2 ) has a ϑ 2 behaviour around the origin when ϑ 1 = ϑ 2 = ϑ, therefore it remains regular even when multiplied with the two K functions. However when only one of the rapidities is close to zero then it has a first order pole F 5 ∝ 1 ϑ , hence taking into account the singularity of the appropriate K function leads to a second order pole. Following the formalism introduced in [68], one can write the sum using contour integrals where the two-dimensional product contour C I ×C J encircles the solution of the Bethe-Yang equation determining ϑ 1 and ϑ 2 with quantum numbers I and J . To open the contours it is necessary to subtract the residue terms when ϑ 1 = 0 or ϑ 2 = 0. Hence where we used that on the contours with imaginary parts −iε e iQ 4,k → ∞ in the infinite volume limit. Now we split C 14 according to the lines in (G.8) as C 14 = C int 14 + C res1 14 + C res2 14 . Clearly, C res1 14 = C res2 14 =: C res 14 and the infinite volume limit of C int 14 is regular: since e iQ 4,k → 0 in the infinite volume limit on the contours with imaginary parts +iε. The residue contribution is e iQ 4,1 + 1 e iQ 4,2 + 1 h({0}|ϑ 1 , ϑ 2 ) R .

(G.9)
For ϑ 2 ≈ 0 based on (D.8) in the contour integral around the ϑ 2 = 0 point the 1/ϑ 2 2 term acts as a differentiation on the other regular ϑ 2 dependent factor, leading to Turning to −Z 1 C 03 we have Therefore the O(L) term in (G.11) is cancelled, resulting in 2ϕ(ϑ 1 )) .

(G.12)
Pulling the contours back to the real axis by sending ε to zero, one can then send R to zero as well with the final result 2ϕ(ϑ 1 )) .

(G.13)
Addressing the time dependence of this term, note that the structure of D 14 is reminiscent of D 12 discussed in Section 3.2 resulting in a √ t type behaviour again due to a mechanism analogous to parametric resonance. The difference from the case of D 12 is that the oscillations are of the form cos 3mt instead of cos mt. Since our primary focus is on one-particle oscillations, we do not discuss this term further here.

G.2 Evaluation of
(G.14) where I and J are the quantum numbers specifying ϑ 1 and ϑ 2 , i.e.
Note that I takes half-integer, while J takes integer values, according to the discussion in Appendix A.1. The expression is singular when ϑ 1 = ϑ 2 or when ϑ 1 = 0 and ϑ 2 is finite, and these singularities are of second order. One can first write the sum over J as a contour integral: where C J surrounds the positions ϑ 2 corresponding to J. Opening the contour leads to 16) where in the first two terms we perform the mL → ∞ only for the rapidities ϑ 2 , from which only the integration along the contour above the real axis survives. We thus have Consider the residue terms , , where (D.9) was made use of. Note, that both the second and first order singularities are to be taken into account. Adding C res1A

23
and C res1B

23
, one has , 20) or after some manipulations , The four terms of C res1 23 have singularities at ϑ 1 = 0. However, these singularities can only produce terms with positive powers in mL but no non-trivial time dependence since h(ϑ 1 |ϑ 1 , {0}) R contains no ϑ 1 dependent function multiplied by t. We therefore have a single secular term, namely C res/sec 23 and taking the limit R → 0 results in C res/sec 23 H Evaluating G 5 , part II. Contour integral terms from D 23 In this appendix we evaluate the contour integral from (G.16) which reads Let us start with a summary of our method first. Manipulating (H.1) leads to various terms of which only (H.10) and (H.31) give rise to interesting time dependence. The origin of these terms was made clear at the beginning of this Appendix, but their actual evaluation needs further non trivial integral manipulations discussed in Appendix I. The largest part of this Appendix is dedicated to showing that apart from (H.10) and (H.31) no other term yields any interesting time dependence. Starting from (H.1) we can subtract and add back the singularities of the five particle form factors. Using (D.9), this leads to where Ω(ϑ) = (1 − S(ϑ)) (1 − S(−ϑ)).
To make the rather technical evaluation more transparent, we introduce an additional simplification which does not affect the end result. Apart from (H.10) and (H.31), many other terms also possess singularities that in principle could lead to non-trivial time-dependence, but cancel each other in the end. These singularities emerge from the region where ϑ 1 is around zero. But at zero rapidity S(0) = −1, and all other quantities behave similar to the Ising model which has a constant S = −1 everywhere, resulting in K(ϑ) being an odd function of ϑ. A good example is the sub-leading singularity of the form factor F ε 5 (ϑ 1 ) defined in (D.9), which for small ϑ 1 behaves as whereas its Ising counterpart is exactly Therefore to keep the reasoning as simple as possible, from now on we perform our calculations for the Ising scattering matrix S = −1 and show the cancellation of certain singularities. It turns out that in the only nontrivial time dependent terms (H.28) and (H.31) the original S matrix can be easily restored.

23
and its descendants

23
, C intBII2 23 , C intBII3 23 and C intBII4 23 from which we have the first term which is regular, and the second term where we used (H.3) and (D.13). This can be further simplified to (H. 19) H.3 The term C intA As a function of ϑ 2 this term has no singularity at the origin so we can pull the contour back to the real axis at once: , thus a stationary phase evaluation using (D.16) yields a 1/ √ t behaviour multiplied by the function value in ϑ 2 = 0. This expression as a function of ϑ 1 has singularities of 4th and 2nd order, but the 4th order ones just cancel due to as F ε 5 (ϑ 1 ) ∝ 8F 1 ϑ 1 around the origin. Hence in ϑ 1 no 4th order singularity is present and when the sum over I is converted to an integral the remaining 2nd order singularity can only produce terms of the type mL but no higher power of L. These are expected to be cancelled by terms from Z 1 D 12 .

H.5 Terms with non trivial time dependence
As seen in (H.4), a large number of terms originating from (H.1) involving integration with respect to ϑ 2 can be evaluated directly using the SPA (D. 16). Once the SPA is performed, these terms have singularities in ϑ 1 , which are of either 4th or second order, and those of 4th order behaviour cancel each other, whereas the second order singularities cannot produce time dependence stronger than O(t 0 ). For some terms, however, the ϑ 2 integration cannot be easily performed. The first example is provided by (H.10) which we write again as (H. 28) and split into C intBI2a and C intBII2b 23 provide the second example, which is written again for better transparency as Ker a stac (ϑ 1 , t, R) =e imt Ω(ϑ 1 ) √ πmt e 2imt(cosh ϑ 1 −1) e −iπ/4 2 cosh ϑ 1 e −mR(cosh ϑ 1 +3/2) .

(I.5)
As for Ker the SPA, (D.16) cannot be directly applied, we proceed as follows: we first differentiate the integrand with respect to t and apply the SPA which becomes now possible: Now integrating with respected to t one ends up with the Fresnel sine and cosine functions, denoted here by F S and F C respectively: where f is an integration constant that is independent of time, but is a function of ϑ 1 and R. We do not determine the precise form of f (ϑ 1 , R), only quote its R = 0 limit: which is determined by noticing that the t → ∞ limit for Ker(ϑ 1 , t, R) is proportional to | sinh ϑ 1 |.
Note that √ cosh ϑ 1 −1) cosh ϑ 1 f (ϑ 1 , R) ≈ ϑ 2 1 around the origin, i.e. its second derivative is continuous at the origin. Hence integrating it with |K| 2 gives a finite and well-defined result. In (I.2), however, Ker(ϑ 1 , t, R) is integrated with a 1/ sinh 4 (ϑ 1 ) type of function due to K(ϑ 1 )F ε 5 (ϑ 1 )+F 1 Ω(ϑ 1 )K ′ (ϑ 1 ), where the non-analytic behaviour of f (ϑ 1 ) must be carefully handled. The origin of the non-analytic term can be summarised as follows: using the SPA, (D.16), a term proportional to 1/ √ t is obtained, but in the asymptotic expansion of the oscillatory integral, terms proportional to t −1/2+n are also present. For any finite t, these lead to analytic behaviour as expected from (I.4), but in the t → ∞ limit keeping only the leading terms, non-analyticity can emerge.
Substituting Ker a stac into (I.1), the discrete sum to evaluate reads that has a 1/ϑ 2 1 singularity at the origin which we treat with a contour integral representation: (I.10) As Ker a (ϑ 1 , t, R) is even in ϑ 1 , its derivative vanishes at the origin, and so only the derivative of e iQ 2 (ϑ 1 ) + 1 gives a pole contribution resulting in and hence expected to be cancelled by the appropriate counter-term from Z 1 D 12 . The contour integral in the L → ∞ limit reduces to the upper contour and can be rewritten using (D.14) as which is integrable since Ker a stac (ϑ 1 , t, R) is an even and regular function with respect to ϑ 1 . Then the derivative of |K(ϑ 1 )| 2 sinh 2 ϑ 1 cosh ϑ 1 gives a √ t contribution which we are not interested in at the moment. The contribution linear in t comes from |K(ϑ 1 )| 2 tanh(ϑ 1 ) (Ker a stac (ϑ 1 , t, 0)) (I. 13) I.2 Time dependence from C intBI2b 23 via Ker: imaginary part As Ker ≈ ϑ 2 1 around the origin, one can naively try to evaluate the integral However, as time t grows, due to the asymptotics of the Fresnel function (F S → 1 2 and F C → 1 2 ), the interval in which Ker b ≈ ϑ 2 1 holds shrinks as [− 1 √ t , 1 √ t ], outside of which Ker b ≈ 0 since ℑm Ker b includes the difference of F C and F S . On the other hand, at the endpoints of the intervals the integrand behaves as t/ϑ 1 → t √ t, so one expects that the integral is linear in t and its coefficient is determined by the small ϑ behaviour of the K function. One can check this assumption numerically and conclude that in the long time limit Similarly to the imaginary part, ℜe Ker ≈ ϑ 2 1 around the origin, hence one can once again try to evaluate the integral directly Again, as the time t grows, Ker ≈ ϑ 2 1 holds in the interval [− 1 √ t , 1 √ t ] , while at the endpoints of the intervals the integrand behaves as t/ϑ 1 → t √ t, so one expects a term linearly dependent on t whose coefficient is given by the small ϑ behaviour of the K function. However, outside this interval the overall behaviour of the integrand is now of the type 1/ϑ 1 accounting for a logarithmic dependence and a t ln t type behaviour with a coefficient related to the the small ϑ behaviour of the K again. Differentiating Ker stac with respect to t and neglecting (−imt − Rm/2) in (I. 16), the SPA (D. 16) can be directly applied resulting in a 1/t term whose coefficient can be identified with that of the logarithmic term. Explicit calculation shows that log(mt) π .

23
To calculate the time-dependence from C intBI−II

23
we cannot use Ker stac directly, i.e. the long-time approximation (I.7) of (I.4) in the sum The reason is that the singularity of |K(ϑ 1 )| 2 F ε 5 (ϑ 1 ) is of order four, whereas for any finite t, Ker stac (and also Ker) behave as ϑ 2 1 around the origin and the resulting 2nd order singularity is sensitive to the non-analyticity of f in the long-time approximation of Ker stac (ϑ 1 , t, R) in (I.7).
With Ker(ϑ 1 , t, R), which is an analytic function for any finite t without such a singular behaviour, one can formally express the time dependence using the contour manipulations and eventually (D.14), yielding where Ker(ϑ, t, 0) ′′ is the second derivative with respect to ϑ 1 and ... refers (mL) 0 terms with no time dependence. Ker(ϑ, t, 0) ′′ can be approximated using the SPA, (D.16) giving π 2 (−1 − i) √ mt, hence the pole term yields an allowed mL √ mt factor. The way to actually evaluate C intBI−II 23 is done by performing numerically the integration in (I.4) and its derivative needed for (I. 20) or to calculate only (I.4) and perform the sum in (I. 19).
After performing the numerical integral, we evaluated the sum for the Ising case and came to the conclusion that the leading order time dependence has the form g 2 F 1 e −imt (Υ 1 mt + iΥ 2 log mt) , (I. 21) for large mt. It turns out, however, that for the linear term one can even keep the real part from Ker stac in (I.7) which has a milder behaviour than f in (I.8) and the contour integral manipulations can formally be performed. The end of the analysis by generalising the Ising result is (I.22) Note, that for the Ising model is a real function, whereas for an arbitrary interacting IQFT it has an imaginary part as well.
As argued at the beginning of Appendix (H), the singular structure at the origin of is the same for the Ising an for any interacting theory, and since the source of time dependence is attributed to the singularity, we can keep the real, i.e. singular part in (I.22). As a consequence, it is possible to extract the linear time dependence from (I.22), and one is allowed to use again the Ising model, where For Υ 1 it is useful to differentiate (I. 22) with respect to t because for the resulting function the SPA (D.16) can be applied yielding a constant and a 1/t type term. Elementary calculation shows that the former term is π . (I.24) The details of the numerical study can be found in the next Appendix.

J Numerical checks for the calculations
Our quite lengthy and tedious analytic calculations were extensively cross-checked using numerics out of which we discuss and present here three important parts: the cancellation of mL terms iñ D 23 , the time dependence from the term C intBI−II (H.31), and match between the time dependence of D 23 and our analytic predictions (3.30). Beyond these, we also numerically verified other parts of the calculations, such as: numerically monitoring the validity of our manipulations performed in Appendix G (using sine-Gordon first breather form factors) and in Appendix I. From Appendix H, the calculations in (H.2) and (H.3) were also cross-checked, whereas the validity of the rest of the Appendix (i.e. cancellation of the 4th order singularities) is verified by the match between the predicted time evolution of D 23 and our analytical considerations, discussed below.

J.1 Cancellation of mL terms inD 23
During the evaluation ofD 23 in Appendix G and H we ignored discussing and showing the cancellation of O(mL) type terms. Here we verify their cancellation by considering the numerical values for D 23 for various time instants and mL system sizes. For simplicity, we use the Ising case S = −1 for our checks with analytic expression for the form factors from [77] and with a singular K function K Ising (ϑ) = −ig 2 sinh 2ϑ . (J.1) Note that the only way to obtain a quench with a zero momentum particle in the Ising model is when quenching from the ferromagnetic to the paramagnetic phase [25,26,27]. For such a quench calculations based on a form factor expansion presupposing a small post-quench density are not expected to give accurate results. Nevertheless, the cancellation of volume dependent terms is related to the order-by-order structure of the expansion independent of the eventual behaviour of the expansion itself. We evaluated numerically for time instants mt = 0, 2.5, 5, 10, 20 and volume sizes mL = 30, 40, 50, 60 performing discrete summation on the quantised rapidities. To make the summation finite we introduced a rapidity cut-off ϑ c = 4.834 for the particles involved. In Fig. (2) we show the numerical values for e imtD23 g/2 with g = 1. Real and imaginary parts of e imtD23 g/2 for time instants mt = 0, 2.5, 5, 10 and 20 for a quench in the Ising model. Results for system sizes mL = 30, 40, 50 and 60 are shown with blue, purple, orange and red symbols, but with a single exception these cannot really be distinguished due to their almost perfect overlap.
The source of the only observable deviation (for the case mL = 30 and mt = 20 ) is numerical inaccuracy resulting from the oscillatory nature of the terms of the sum. In Fig. (3) we plot the difference between the numerical values obtained for various system sizes at fixed times. These results demonstrate that terms growing with mL indeed cancel.
(J.8) Figure 4 shows the numerical results for (J.7) for system sizes mL = 50, 60 and for various time instants. The kernel Ker was determined by numerical integration and we subtracted the mL residue term   Figure 5 shows the real and imaginary factor multiplying the term g 2 F 1 e −imt in C intBI−II calculated by a discrete summation (J.7) and by the integration (J.8). For the imaginary parts, we also keep the F S + F C part from Ker stac (I.7) but drop the √ sinh 2 ϑ type function, f . We can conclude that for the real part multiplying g 2 F 1 e −imt using the real part of Ker stac in (J.8) gives a correct result. For the imaginary part, however, the imaginary part without the singular f from Ker stac alone is not able to reproduce the result of the discrete summation. On the other hand, the time dependence of the imaginary coefficient is only logarithmic, and therefore the calculation of this sub-leading time dependence is not addressed in this work. We also checked if the linear time dependence can be described by as predicted in Appendix I and if the logarithmic time dependence for the imaginary part is a correct assumption. In Fig. 6 we therefore fit the functions a+b mt and a+b ln mt to the real and imaginary parts of e imt C intBI−II F 1 g/2 − mL Ker ′′ g 4 32 calculated by the discrete summation with mL = 60 omitting the first 3 data points corresponding to the shortest times mt = 0.01, 1, 2.5. calculated by a discrete sum using Ker stac for various time instants for a quench in the Ising model. The red dots correspond to system sizes mL = 60 and the blue line to the fitted curves of type a + b mt and a + b ln mt.
From the linear regression, the coefficient of the linearly time-dependent term is 0.0799167 which is to be compared with 1 4π = 0.0795775 as g = 1. The agreement is excellent, and although the error of the fitted parameter is 7 × 10 −5 , the match is convincing.

J.3 Comparing the time dependence ofD 23 with the analytic results
In this subsection we study the time dependence of e imtD 23 g/2 for the Ising model with K = −i g 2 sinh 2ϑ for mL = 70. We compute this quantity using discrete finite volume summation for various time instants, and fit its real and imaginary parts with the appropriate function dependences a + b √ mt + c mt + d mt ln mt and a + b √ mt + c mt dictated by our analysis done in Appendix I, omitting the first 7 data points corresponding to short times. Concerning the real parts, the values obtained from the fit are d = −0.0825821 and c = 0.149943, which must be compared with − 1 4π = −0.0795775 and K + 3 4π = 0.131292 resulting from (3.30) and (3.31) with g = 1 and S = −1. The accuracy of the fit itself is around 10 −3 − 10 −4 . The difference between the fitted parameters and the analytic predictions is now a bit larger compared to the previous subsection. However, neglecting more data points for short times the fitted values move towards the analytic predictions. We note however, that omitting too many points leads to a deterioration of the fit quality, as for an accurate determination of a logarithmic term data points over several orders of magnitude should be used, which is not possible to extract due to the inaccuracy of the discrete sum for large times. As an additional test we also tried the fitting function a + b √ mt + c mt and noted that the fit residuals were two orders of magnitude larger than for the case including the term mt ln mt, which is a confirmation of the presence of the term mt ln mt in the time dependence.
For the imaginary part the parameter c was found to be c = 0.114809 which must be compared with 1 8 = 0.125 for g = 1 according to (3.30) and (3.31). The total estimated error of the fitted value is of the order 10 −3 . In principle we should have either included ln mt in the fitting function or subtracted the contribution of ℑm e imt C intBI−II F 1 g/2 − mL Ker ′′ g 4 32 discussed in the previous subsection from the data points, but since this correction is rather small it was simply discarded.