Initial states in integrable quantum field theory quenches from an integral equation hierarchy

We consider the problem of determining the initial state of integrable quantum field theory quenches in terms of the post-quench eigenstates. The corresponding overlaps are a fundamental input to most exact methods to treat integrable quantum quenches. We construct and examine an infinite integral equation hierarchy based on the form factor bootstrap, proposed earlier as a set of conditions deter- mining the overlaps. Using quenches of the mass and interaction in Sinh-Gordon theory as a concrete example, we present theoretical arguments that the state has the squeezed coherent form expected for integrable quenches, and supporting an Ansatz for the solution of the hierarchy. Moreover we also develop an iterative method to solve numerically the lowest equation of the hierarchy. The iterative solution along with extensive numerical checks performed using the next equation of the hierarchy provide a strong numerical evidence that the proposed Ansatz gives a very good approximation for the solution.


Introduction
The study of quantum dynamics in one-dimensional integrable systems has led to intriguing discoveries, like the experimental observation of lack of thermalization [1,2,3], the theoretical prediction [4] and experimental observation [5] of an unconventional statistical ensemble known as the Generalized Gibbs Ensemble (GGE), the discovery of unexpected effects on quantum transport [6,7] and of novel quasi-local conserved charges [8,10,9]. A typical protocol employed for the preparation of a closed quantum system in an out-of-equilibrium initial state is the instantaneous change of some parameter of its Hamiltonian, a process called quantum quench [11,12]. After a quantum quench, the initial state of the system, which is typically the ground state of the pre-quench Hamiltonian, evolves unitarily under a different post-quench Hamiltonian.
Obviously in order to derive the time evolution of the system we generally need to know the post-quench excitation amplitudes in the initial state, i.e. the overlaps of the initial state with any post-quench energy eigenstates. Determining the excitation content of the initial state is typically easy for quantum quenches in which both the pre-and the post-quench Hamiltonians are quadratic in terms of some suitable physical fields. In this class of problems, which correspond to non-interacting models or interacting models that can be mapped into non-interacting ones, the relation between pre-quench and post-quench excitations is typically described by a so-called Bogoliubov transformation and the initial state is then a squeezed coherent state, more precisely a squeezed vacuum state, of the post-quench Hamiltonian [13,14,15]. Such states consist of pairs of excitations with opposite momenta and have the characteristic exponential form where A † (k) are operators that create particles of momentum k and |0 is the corresponding vacuum state. However, the task of determining the excitation content of the initial state is far more difficult in the context of genuinely interacting integrable systems (except in special cases [16]). The reason is that the pre-quench and post-quench excitations are no longer related through a simple linear transformation, but instead a nonlinear one that in general corresponds to an infinite series [17]. Even though the initial state has been derived exactly in a number of special cases [18,19,20,21,22,23], a general and systematic method for its determination remains so far elusive. On the other hand, while in those earlier studies the initial state was derived by means of finite volume calculations based on the Bethe Ansatz, in fact we are mostly interested in the thermodynamic limit where all finite size effects vanish and integrable quantum field theoretic methods come into play [24,17,25,26,27]. In particular, it has been argued [28] that in the thermodynamic limit, i.e. in the limit of large system size L and particle number N with a fixed density of particles, it is sufficient to know only the extensive part of the logarithmic overlaps between the initial state |Ψ and the post-quench eigenstates |Φ , i.e. the quantity Indeed it was shown in [28] that a single post-quench eigenstate that is representative of the initial state, is sufficient for the description of the asymptotic values of local observables at times t → ∞. Such a representative state is completely determined by the quantity E Ψ [Φ] above, moreover, the full time evolution can also be derived from the same quantity. The argument is based on the fundamental idea of statistical physics that microstates can be classified according to their macroscopic properties into macrostates, which in integrable models are fully characterized by their Bethe Ansatz root density [29]. This approach, also known as Quench Action method, has successfully predicted the stationary values of observables at large times [30,31,32,33,34,35] as well as their full time evolution [36,37,38,39]. The thermodynamic quantity (1.2) for any specific quench initial state is generally given by an expression much simpler than the exact expression for the finite volume overlaps [19], suggesting that it may be easier to derive it directly in infinite volume by a field theoretic method. One of the earliest field theoretic approaches was given in [11,12] in which the quantum quench problem was reformulated as a boundary problem defined in an infinite strip geometry with boundary conditions given by the initial state. This way the authors exploited results already known from Boundary Renormalization Group (RG) theory to study a broad class of quantum quenches in Conformal Field Theories (CFT). In particular it is known that universal boundary behaviour at equilibrium can be classified within a small number of different types, each corresponding to a scale invariant boundary state. For example, the Dirichlet boundary state |D corresponds to a quantum quench with vanishing initial field fluctuations, as in the case of an infinite pre-quench particle mass m 0 → ∞.
Exploring the generalization of this approach to systems described by an Integrable Field Theory (IFT) with no off-diagonal scattering, the authors of [24,40] studied a special class of initial states of the form where Z † (θ) is an operator that creates a post-quench excitation of rapidity θ, |Ω is the ground state of the IFT, while K(θ) is the pair excitation amplitude. States of this form, which is a generalization of the squeezed vacuum state (1.1), include (but are not limited to) the so-called boundary integrable states [41], i.e. boundary states that respect the integrability of the bulk. This approach was later extended to include a slightly more general form of initial states [42]. For this type of states it was possible to analytically show that the system equilibrates and is described by a GGE. Furthermore the large time decay of observables can also be analytically calculated for quenches starting from an analogous class of initial states in the sine-Gordon model [36], using a form factor based approach. In addition, using the semiclassical approach proposed in [43] and extended to a quench situation in [44], the time evolution of correlation functions after a quench in the sine-Gordon model starting from a state of the form (1.3) was determined in [45]. Moreover, in recent investigations of the Loschmidt echo and statistics of work done during a quantum quench [46,47], an initial state of the form (1.3) was considered as a starting point. Therefore the knowledge of the amplitude K, which fully characterizes the state (1.3), is sufficient to determine the expectation values of local operators in the post-quench stationary state. In addition, the amplitude K serves as an input quantity to form factor based and semiclassical approaches to correlators, and to determinations of Loschmidt echo and the statistics of work. As a result, knowledge of K is of primary interest in the description of integrable quantum field theory quenches.
However, there is no a priori reason to believe that a quantum quench in an IFT would lead to an initial state of this form. Indeed, even though Dirichlet states are of this form and correspond to a quantum quench with an infinite pre-quench particle mass m 0 , these are ill-defined states: they exhibit ultraviolet divergent physical observables due to the fact that their excitation amplitude K D (θ) does not decay for large momenta. This is easily understood: initial state excitations are cut-off at large momenta by the mass scale m 0 which in this case is taken to be infinite. To regularize the boundary state in the CFT case, the authors of [11,12] used the concept of "extrapolation length" τ 0 , known from Boundary RG theory, where it expresses the difference of the actual boundary state from the idealized Dirichlet state. The parameter τ 0 plays the role of an exponential large momentum cut-off, dependent in general on the initial parameters and typically proportional to 1/m 0 in the CFT case. A generalization of this idea to IFT [24] (although not justified by RG theory since massive IFT are non-critical) would amount to replacing the Dirichlet state |D by the regularized state e −Hτ 0 |D where H is the post-quench Hamiltonian, or equivalently replacing K D (θ) by where E(θ) is the energy of an excitation with rapidity θ. In the IFT case τ 0 can be considered as a phenomenological parameter. Although it may be expected to be of order 1/m 0 as in the CFT case, its precise dependence on the quench parameters is not generally known. The relation of this phenomenological parameter to the quench parameters is crucial, since the values of physical observables depend explicitly on it. In [25] an estimation of τ 0 in the free massive case by comparison of the field fluctuations in the actual and approximate state demonstrated that this approach does not reproduce correctly the known exact results for the large time values of observables, as far as numerical factors are concerned. Moreover it was shown that the estimation of τ 0 is not independent of the choice of observable used to make the comparison between actual and approximate state. Different choices lead to different scaling for the m 0 dependence of τ 0 as m 0 → ∞. This discrepancy is an indication that the effect of large momentum excitations that are present in the initial state cannot be incorporated in a suitable and unique definition of the momentum cut-off. On the other hand, there is no fundamental reason preventing us from choosing τ 0 to be momentum dependent itself. In this way the actual K(θ) may not necessarily decay exponentially with the momentum and such alternative choices would modify the predicted values of physical observables. In fact, such a generalization is justified from the point of view of the same boundary formulation: RG theory teaches us that, in estimating the difference of the actual boundary state from the idealized Dirichlet state |D , any boundary irrelevant operator could be inserted as a boundary perturbation [48]. In an IFT such boundary irrelevant operators include (but are not limited to) all conserved charges of the bulk theory. Indeed adding such perturbations does not critically change the system's behaviour [49]. This means that, in the same way that the extrapolation length τ 0 is introduced essentially as a perturbative parameter associated to the Hamiltonian, one could in principle introduce a different parameter τ s for each conserved charge Q s =´dθq s (θ)Z † (θ)Z(θ). This generalizes our Ansatz for the regularized initial state from e −Hτ 0 |D to the more general e − s Qsτs |D , which is clearly equivalent to introducing a momentum-dependent "extrapolation length" τ (θ). Such a generalized regularization K D (θ) → K D (θ)e −2E(θ)τ (θ) could be any function of θ that fulfills the ultraviolet convergence condition. Other irrelevant perturbations that may be included as perturbations are not simply quadratic in Z(θ) and Z † (θ) (as are the charges) but of higher order instead [48]. These would lead to deviations from the form (1.3). Overall this conclusion brings us back to our starting question of how to determine correctly the pair excitation amplitude K(θ) of (1.3), or more generally all (possibly independent) excitation amplitudes of the initial state.
Note that, as mentioned above, this is a physically important problem in the context of quantum quenches, since the values of observables depend significantly on the initial state and even its regularization. In particular a modification of the amplitude K(θ) may crucially affect the subsequent evolution and the final equilibrium. In the CFT case, for example, considering only the Hamiltonian as perturbation leads to the conclusion that the system tends to an effectively thermal equilibrium with a temperature proportional to τ 0 [11,12]. On the contrary, including more perturbations results in the emergence of a generalized equilibrium (GGE) with as many different temperatures as the perturbative parameters [48]. General studies in interacting to free quenches show that in the case of massive evolution the only information that survives at large times is that of the initial two point correlation function and the equilibration is described by the GGE [50]. On the other hand, in the massless case this is not true and much more information about the initial state must be taken into account [51], unless the initial state is gaussian in which case the GGE is still valid [52,53,54,55].
In [25] a systematic approach was proposed to determine the expansion of the initial state in terms of post-quench excitations from first principles. The specific problem under investigation was a quantum quench of the particle mass m and the interaction g in the Sinh-Gordon model, a prototypical integrable model with a single type of particle, starting from a large but not infinite initial mass and zero interaction. It was shown that the condition that the initial state, being the ground state of the pre-quench Hamiltonian, is annihilated by all pre-quench annihilation operators, imposes an infinite system of equations that must be satisfied by the initial state excitation amplitudes. These are integral equations that involve infinite series of form factors of the physical field φ. The simplest of these equations for the pair excitation amplitude K(θ) was derived explicitly, assuming that the form (1.3) is valid. By truncating the form factor series and analyzing the integration kernels, a simple factorized Ansatz for the solution K(θ) was proposed and numerically verified. This justified the expression (1.3) for the initial state as the leading order result in a systematic expansion, rather than an ambiguous approximation. On the other hand it also showed that a consistent regularization of K D (θ) for large initial mass is decaying not exponentially as suggested by (1.4), but rather as a power of the momentum.
The Ansatz of [25] itself rested on a number of assumptions, for which only partial justification was given, and the numerical checks only treated a single member of an eventually infinite integral equation hierarchy, which the amplitude has to satisfy. In this paper we first derive the whole hierarchy of equations in explicit form at all orders of the form factor series. To achieve this we employ a simpler derivation method that works directly in the infinite volume limit and verify its agreement with our earlier method that is based on finite volume regularization. Second, we eliminate a large part of the assumptions present in [25] and give plausible arguments for the rest. In addition, we perform a thorough numerical analysis of the hierarchy. Unlike [25] where it was only checked that the Ansatz satisfies approximately a truncated version of the lowest order equation, here we compute a numerical solution of the equation by means of a newly developed iterative method without bias (i.e., without assuming the Ansatz). Moreover we perform numerical consistency checks that the derived solution also satisfies the next order equation and finally compare it with our Ansatz concluding that they are in perfect agreement. We also give independent numerical evidence for the correctness of our theoretical arguments. While we treat the same quantum quench in Sinh-Gordon model, as before, our theoretical considerations, supported by the numerics, also make it plausible that the expression (1.3) for the initial state, and the factorized form of the Ansatz is valid for a much larger class of models.
The paper is organized as follows. After setting up some necessary notions and notations in Section 2, we proceed to a general overview of the hierarchy of integral equations in Section 3, and consider the problem of existence and uniqueness of its solution. The properties of the solution are explored in Section 4, using both general principles of field theory and reasoning connected to integrability, leading to a strong plausibility argument for the Ansatz used in [25]. Then we proceed to the numerical investigation of the integral equation hierarchy in Section 5, and draw our conclusions in Section 6. The paper also contains two appendices, Appendix A containing concepts and results of finite volume regularization that are necessary for the main text, and Appendix B collecting some numerical data for illustration.

The Sinh-Gordon model
The Sinh-Gordon theory is defined by the Hamiltonian where µ is the classical particle mass and g the coupling constant. It is the simplest example of Lagrangian integrable field theories, and is invariant under the Z 2 symmetry φ → −φ. The spectrum of the model is made up by multi-particle states of a single massive bosonic particle, whose exact mass at the quantum level is denoted by m. Parameterizing the dispersion relations in terms of the rapidity θ, the energy and momentum of a single particle excitation are E(θ) = m cosh θ and p(θ) = m sinh θ and the two-particle S-matrix is given by [56] S(θ, B) = tanh where θ = θ 1 − θ 2 is the relative rapidity of the particles, and B is the so-called renormalized coupling constant For real values of g the S-matrix has no poles in the physical sheet and hence there are no bound states. A complete basis of eigenstates of this Quantum Field Theory is provided by the n particle asymptotic states where the operator Z † (θ) creates a particle excitation with rapidity θ and |Ω is the vacuum state of the theory. The creation and annihilation operators Z † (θ) and Z(θ) satisfy the Zamolodchikov-Faddeev algebra The form factors F O n of the Sinh-Gordon model are matrix elements of a generic local operator O(0, 0) between the vacuum and a set of n particle asymptotic states Any other matrix element of the operator O(0, 0) can be obtained from the form factors above by exploiting the crossing symmetry of quantum field theory. The latter imposes that matrix elements of O(0, 0) between in-and out-states can be obtained by analytical continuation of the out-state rapidities, as discussed in detail in Smirnov's book [57]: where the sum is over all possible splittings of the sets Θ = {θ 1 , ..., θ n } and B = {β 1 , ..., β l } into nonoverlapping subsets Θ 1 , Θ 2 and B 1 , B 2 respectively (with Θ 2 non-empty), S(Θ, Θ 1 ) and S(B, B 1 ) the S-matrix products needed to reorder the rapidities and ∆(Θ 2 , B 2 ) = Θ 2 | B 2 . Note that using the translation operator U = e −i(Ht−P x) we can always shift any operator O(x, t) to the origin O(x, t) = U † O(0, 0)U . Based on general properties of a Quantum Field Theory (as unitarity, analyticity and locality) and on the relations (2.5), the form factor bootstrap approach leads to a system of linear and recursive equations for the matrix elements F O n [57] Form factor solutions for the Sinh-Gordon model were constructed in [58,59]; explicit expressions of the form factors F φ n (θ 1 , θ 2 , ..., θ n ) for the elementary field φ can be found in the Supplementary Material of [25].
Observable operators can be written as an expansion in terms of the Zamolodchikov-Faddeev operators [57] (for a more recent exposition in the framework of algebraic QFT cf. [60]) where the functions f can be expressed in terms of the form factors and Z † , Z are the Zamolodchikov-Faddeev creation and annihilation operators, which satisfy the algebra (2.5). It can be easily verified that the above expansion reproduces correctly the form factors of the local field O.
The above formula is equivalent to (2.7). To show that is sufficient to consider the general matrix element of the operator O, substitute (2.10), perform the contractions using the algebra (2.5) and compare the result to (2.7). The simplest case is the matrix elements with no particle in the out-state, as those appearing in the definition of form factors (2.6). We can easily see that only one term of the expansion (2.10) contributes to such a matrix element and therefore we obtain Notice that the ordering of contractions imposes that the order of rapidities gets inverted. Considering a general matrix element instead, there are also disconnected contributions from all lower order terms of the expansion. Therefore to determine the coefficients f O l,n one can start from the l = 0 case obtained above and work recursively in l, deriving each new one using all earlier derived coefficients. It turns out that the disconnected contributions match one by one the disconnected terms of Smirnov's formula, and we find that (2.11) holds.
From (2.8) the functions f O l,n satisfy the permutation relations We also introduce the Dirichlet boundary states, defined as the eigenstates of the field operator located at the temporal boundary φ(t = 0). Below we only need the Dirichlet state |D that corresponds to zero field value, i.e. defined by the condition (2.14) This is an integrable boundary state, i.e. its exact expansion on the eigenstate basis (2.4) is of the squeezed vacuum form [41] |D = exp with amplitude K D (θ) given by This formula can be obtained by analytically continuing the result for the first breather of the Sine-Gordon model, obtained in [61], to imaginary coupling. Under such continuation the Sine-Gordon model turns into the Sinh-Gordon model to all orders of perturbation theory, and the first Sine-Gordon breather is mapped into the Sinh-Gordon particle, both of which are physically identical to the elementary particle created by the field.
3 The integral equation hierarchy for the initial state

The infinite hierarchy
We consider a quench in the Sinh-Gordon theory from an arbitrary initial mass m 0 and zero initial interaction to arbitrary (renormalized) mass m and interaction coupling g. The initial state |B is the ground state of the pre-quench Hamiltonian i.e. it is defined through the condition that it is annihilated by the pre-quench annihilation operators for all momenta p. It is convenient to express the annihilation operators A(p) in terms of the elementary field φ and the canonical conjugate field π, which by continuity of the time evolution through the quench satisfy the initial conditions This can be done easily since the pre-quench annihilation operators A(p) correspond to a free field theory and are therefore given in terms of φ and π by the standard relation is the Fourier transform of φ(x) and is the energy of a pre-quench excitation with momentum p. Moreover, from the relation π where H is the post-quench Hamiltonian, we obtain the equation This equation was first derived in [25] and, as we will now see, it has the advantage that all operators are straightforwardly expressible in the framework of the post-quench field theory, i.e. the Sinh-Gordon theory.
Expanding the initial state on a complete set of states we have where |β 1 , . . . , β r are the multi-particle eigenstates of the post-quench Hamiltonian with eigenvalues Without loss of generality, the functions K r (β 1 , . . . , β r ) are chosen to satisfy the symmetry relation by virtue of the commutation rules of the Zamolodchikov-Faddeev algebra (2.5). Projecting (3.6) to all the linearly independent states of the post-quench Hilbert space, we can write an infinite hierarchy of integral equations for the functions K r by taking the inner product with all possible post-quench eigenstates |θ 1 , . . . , θ N . We expect that this set of equations determines all the functions K r uniquely (up to overall normalization). Note that due to translational invariance of the global quench, the initial state |B is in the subspace of zero total momentum, i.e all amplitudes K r contain a delta function δ r j=1 m sinh β j . As a result, the nontrivial equations are only obtained for i.e. when the momentum p is opposite to the total momentum of the test state. Also note that the number N of particles in the test state must be chosen odd, otherwise the equations are trivially satisfied. This can be seen easily from (3.10) since the operatorφ(p) + φ (p), H /E 0 (p) is odd i.e. antisymmetric under the transformation φ → −φ, while the state |B is even since the pre-quench Hamiltonian is symmetric under the same transformation. This means in particular that it contains only excitations with even number of particles, i.e. all odd amplitudes K r in (3.7) must vanish.

The hierarchy in terms of form factors
To derive explicit expressions for equations of the hierarchy (3.10) we substitute the expansion (2.10) for the operatorφ(p) + φ (p), H /E 0 (p) and the general expansion (3.7) of the state |B into (3.10). Then the general equation with an N particle test state is the sum of all possible contractions between the N particles of the test state with the operatorφ(p) + φ (p), H /E 0 (p) and the state |B . From (3.4) and using the translation operator e iP x (where P is the momentum operator) to shift the field φ(x) to the origin, we can easily find that with expansion coefficients f φ l,n (θ 1 , . . . , θ l |η 1 , . . . η n ) given by the form factors of the elementary Sinh-Gordon field φ according to (2.11). We also need the relations and that give the action of the Hamiltonian operator H on the bra and ket states. We can now substitute the above relations to (3.10) ∞ l,n,r=0 and use the Zamolodchikov-Faddeev algebra (2.5) to perform the contractions. From the latter it is easy to see that the matrix element is only non-zero when all η j are contracted with some of the β k and all ζ i are contracted with some of the θ s , while the remaining β k are contracted with the remaining θ s . This means in particular that l and n are constrained to vary between the values 0 ≤ l ≤ N , 0 ≤ n ≤ r and to satisfy the equation N − l + n − r = 0. These conditions restrict the sums in the above equation. Moreover, according to earlier comments, l + n must be an odd positive integer, r must be even and N odd. There are r!/(r − n)! ways to contract η j and β k , N !/(N − l)! ways to contract ζ i and θ s and (r − n)! ways to contract the remaining β k and θ s , which are all equivalent up to S-matrix factors due to permutations of the rapidities. Note that S-matrix factors due to permutations of integrated rapidities can always be absorbed by re-ordering the rapidities to a fixed ordering using (2.13) and (3.9) and by renaming them since they are dummy integration variables. However there still remain S-matrix factors corresponding to the permutations of the test particle rapidities θ s that are necessary in order to choose those that are contracted with ζ i . That is, there are N l distinct choices that give terms multiplied by different S-matrix products. Overall the number of combinations of contractions that give equal contributions is After all contractions have been performed, (3.16) becomes where "perm.s" denotes the other choices of splitting the test particle rapidities θ s into those contracted with ζ i and those contracted with β k , which contain S-matrix products as explained above. Taking into account that the amplitudes K r contain a factor δ r j=1 m sinh β j expressing the translational invariance of |B , the δ-function in the above equation can be replaced by δ p + N s=1 m sinh θ s which can be taken out of the integral and sum. This overall factor means that, as mentioned earlier, nontrivial equations are only those for which p is opposite to the total momentum of the test state (3.11). Provided that this condition is fulfilled, the final form of the N test particle equation is In particular, for a single test particle, there are only two possible values for n, either zero or one, and so n = r − 1 or n = r respectively. Therefore the N = 1 equation is  with p(θ) = m sinh θ and therefore E 0 (p(θ)) = m 2 sinh 2 θ + m 2 0 . In the limit of large initial mass m 0 → ∞, the factors in round brackets in the above equations are dominated by E 0 (p(θ)) ∼ m 0 and therefore they can be replaced by m 0 , which, being an overall factor multiplying the equation, can be omitted. As already explained, in this limit the equations must be satisfied by the Dirichlet state given by (2.15) and (2.16).
The above equations can be represented diagrammatically as in Figure 3.1. The terms (3.22) in the expansion of the operatorφ(p) + φ (p), H /E 0 (p) are represented as square boxes with l legs on the left and n legs on the right, which should be contracted respectively with the N rapidities of the test state, represented as external legs on the left side of the graph, and with the r rapidities of K r , represented as legs emerging from the rectangular box on the right. The remaining rapidities of K r should be contracted with the remaining rapidities of the test state. The sum is over all possible orders n, r and all possible combinations of contractions. Note that in principle the ordering of rapidities in such graphs matters and that exchange of two consecutive rapidities results in multiplication with the corresponding S-matrix factor.
The above derivation was based directly on the infinite volume formulation of Integrable Field Theories through the crucial use of the Zamolodchikov-Faddeev operators and the related operator and state expansions. In order to verify the validity of the equations from the finite volume formulation of integrable models, in Appendix A we present an alternative derivation of the above equations (up to a certain order) based on the Bethe-Yang equations. We compare the two systems of equations and confirm that they are indeed identical.

The hierarchy for free theory: Bogoliubov transformation
To argue for the statement that the hierarchy uniquely determines the initial state up to normalization, consider first the case when the interaction is zero, corresponding to the free bosonic theory with the Hamiltonian where the canonical momentum is π = ∂ t φ. We can define creation and annihilation operators as follows 1 where E(p) = m 2 + p 2 and the canonical commutation relations are Now consider a global quantum quench when at time t = 0 the mass parameter is abruptly changed in the theory m 0 → m. As the time evolution of the fields φ and π is continuous, taking their spatial Fourier transform at time t = 0 and defining the creation and annihilation operators with them, the following equations must hold where E 0 (p) = m 2 0 + p 2 , E(p) = m 2 + p 2 and A 0 (p) and A(p) are the pre-and post-quench mode operators. This is nothing but the familiar Bogoliubov transformation, which allows one to express |B using the fact that it is the pre-quench ground state, i.e. A 0 (p)|B = 0. We now demonstrate that the usual expression for the initial state can also be obtained from the integral equation hierarchy. For better transparency of the subsequent manipulations we start directly with the operator form of the hierarchy (3.6) and go through essentially the same steps that yielded (3.20) from (3.6).
Using (3.24), the operator form of the hierarchy can be written as The fact that A(p) and A † (−p) occur together means that the equations only link components which differ by pairs of particles of opposite momenta. The lowest component is the post-quench vacuum, so |B can be expressed in terms of states composed entirely of pairs of particles with opposite momenta. This allows us to write the amplitude K f ree 2n (−k 1 , k 1 · · · − k n , k n ) as a function of only n variables K f ree 2n (k 1 . . . k n ), hence .3: Diagrammatic expansion of the integral equations for a free theory. The first line represents the one-particle equation, while the second line shows the general multi-particle equation. As can be easily seen in the third line, a factorized solution K 2n+2 ∼ K 2n × K automatically satisfies the multi-particle equation which then reduces to the one-particle equation for the function K. The latter is therefore equal to the pair amplitude K = K 2 . This explains the exponential form of the solution (3.35).
In order to find K 2n , eqn. (3.27) can be further specified by applying a given test state on the left side of eqn. (3.27), which gives the individual equations of the hierarchy. To find K 2 , we apply a one-particle test state p 1 | = 0|A(p 1 ) to eqn. (3.27) and substitute the expansion of (3.28) to obtain where the matrix elements can be evaluated using the canonical commutation relations (or, equivalently, Wick's theorem) which leads to from which follows.
To determine the functional form of K f ree 2n (k) for a general n, test states with higher number of particles are applied on the left side of (3.27); only test states with odd numbers of particles give a nonzero result. Let there be 2n − 1 particles in the test state; then only two terms from (3.28) give a nonzero contribution, which can be attributed toφ(p) containing only one creation and one annihilation operator in its expansion (3.13) in the free theory, or equivalently toφ(p) having non-zero form factors f φ l,n only when l = 1, n = 0 or l = 0, n = 1. These considerations lead to Due to the symmetry properties of K f ree 2n , and using (3.31) This is exactly the squeezed state form that can also be obtained via a direct application of the Bogoliubov transformation (3.26) to the condition A 0 (p)|B = 0 (see e.g. [62] for a derivation in the quantum quench context).

Uniqueness of the solution for the interacting case
As we have seen, the hierarchy (3.10) has a unique solution (up to normalization) for a mass quench in a free field theory, which coincides with the well-known squeezed state resulting from the Bogoliubov transformation. We now argue that for the interacting case (Sinh-Gordon theory) we expect the same uniqueness property. The interacting systems is much more involved than the free one since each of the equations (3.16) involve all of the K r excitation amplitudes and therefore it does not have the chain-like organization of the equations of the free case, where only two terms were present for each equation, allowing us to calculate each amplitude of arbitrary order one-by-one from its predecessors. However, since the zero-coupling limit gives back the free equation, and at least for a small enough value of the coupling the dynamics and form factors of Sinh-Gordon theory are known to be well-described by perturbation theory, the introduction of a small coupling does not spoil the uniqueness of the solution. If anything, we expect the solution to be slightly deformed from the free solution, but to preserve much of its properties. Therefore if a given Ansatz can give a solution, or at least a very good approximation of a solution, we can argue it is close to the unique solution of the hierarchy.
We note that a remnant of the chain structure is still present in the interacting case. Denoting the number of particles in the bra state by n L and in the ket state by n R , for a free field theory the only nonzero matrix elements are the ones with n L − n R = ±1. From perturbation theory, it is clear that these terms dominate for weak coupling. However, in a two-dimensional field theory this argument goes even further: as a simple consideration of available phase space shows, the Feynman graphs corresponding to form factors with larger differences in the number of incoming and outgoing particles are suppressed. The reason for this is that the difference in number of particles can be interpreted as particle creation by the operator inserted; however, in two space time dimensions the available phase space eventually decreases with increasing particle multiplicity [63], therefore processes involving a change in particle number are suppressed; the larger the change, the stronger the suppression. As a result, we expect that the terms are hierarchically organized by the value of ∆n = n L − n R , the ones with ∆n = ±1 being the largest, followed by the terms with ∆n = ±3 and so on 2 . This fact is important for our ability to treat the infinitely many integral equations, each composed of an infinite number of terms, making up the hierarchy, as it implies that they can be well-approximated by equations that are truncated to a finite number of terms, and more terms can be gradually included to improve the approximation.
4 General properties of the initial state

Pair structure and exponentiation of the initial state
In the introduction we presented an argument for the exponentiation of the initial state, based on RG theory. More specifically we argued that the initial state after a quench in an IFT can be built by inclusion of the IFT charges Q s as perturbations to the Dirichlet state: e − s Qsτs |D . Such perturbations preserve the exponential form and merely modify the pair excitation amplitude K(θ). Below we will show an independent argument based on two general properties of the initial state: the extensivity of the charges in the initial state and its pair structure. The first property is due to the fact that the charges are local and the quench is a change of a global parameter. The second property is expected to hold from perturbation theory considerations. Combination of the two properties leads to exponentiation as described by (1.3) and simplifies our search for a solution to the hierarchy equations.

Requirements from extensivity of local charges
We consider a general state |Ψ that satisfies translational invariance and therefore contains only excitations with zero total momentum. Without loss of generality we can write such a state in the form of a cumulant expansion where the amplitudes K Ψ r (θ 1 , ..., θ r ) contain a factor δ( i p(θ i )) which expresses the conservation of total momentum due to the translation invariance. We further assume that the state |Ψ is Z 2 invariant, i.e. symmetric under the transformation φ → −φ, therefore it must contain only excitations with even number of particles or, in other words, all odd amplitudes K Ψ 2n+1 (θ 1 , ..., θ 2n+1 ) must vanish. We will show that the expectation values of the local charges Q s =ˆdθ e sθ Z † (θ)Z(θ) (4.2) in the above state are extensive quantities (i.e. they increase linearly with the system size L) only if the amplitudes K 2n do not contain any other δ-function factor except of the one that accounts for the translation invariance, δ( i p(θ i )). The initial state after a quantum quench, by definition, must satisfy the requirement of extensive local charges. Indeed, the charges Q s , being spatial integrals of local operators, are extensive thermodynamic quantities for all states that satisfy the cluster decomposition principle, as all ground states of physical Hamiltonians do.
The expectation values of Q s in the state |Ψ are and correspond to the sum of all possible ways to contract left and right excitations in the expansion of (4.1) with the charge operator and with each other, in such a way that no part is disconnected from the rest. Diagrammatically this is represented by fully connected graphs as in Fig.4.1. In order to determine the scaling of thermodynamic quantities with L we should simply count the number of redundant δ-function factors of momentum variables in the expectation values of the corresponding operator. As mentioned above, the amplitudes K Ψ 2n contain a factor δ( i p(θ i )) due to the translation invariance. Taking into account all of those δ-function factors involved in each connected graph and assuming that the K Ψ 2n do not contain any other δ-function factors, we can easily see that each graph has exactly one left over factor of δ(0), which is nothing but a factor equal to the system size L. This means that the contribution of such graphs is linear in L i.e. extensive. If however some K Ψ 2n contains two δ-function singularities in the  momentum variables, then the corresponding graphs contain one extra δ(0) ∼ L factor i.e. scale more than linearly with L and are not extensive. For any single charge Q s , extensivity may still be restored by cancellation among the contributions scaling with a higher than linear power of L. Therefore we must assume that the available charges Q s are functionally independent and complete; then the above argument means that all of the graphs must be extensive in order for the charges to be such and this is true only if there is no amplitude K Ψ 2n with more than one δ-function. The same condition ensures the extensivity of the generating function log Z = log Ψ|Ψ (analogous to the free energy) from which other thermodynamic quantities can be derived.
Note that the argument requires the functional completeness of the charges, which may require taking into consideration quasi-local charges recently introduced in [8,10,9]; similar extension of the class of charges was also advocated for field theories in [64]. However, as long as these additional charges are extensive for large system sizes, such as the case e.g. for the charges recently used in completing the GGE for the XXZ spin chain [9], the argument is left unchanged.

Pair structure from integrable dressing
As already mentioned, due to the symmetry of the pre-quench Hamiltonian under the Z 2 transformation φ → −φ, the initial state consists only of excitations with even number of particles. However in the present case it is expected to satisfy an even stronger condition: to consist solely of pairs of particles with opposite momenta. An argument for this can be formulated by considering how switching on the integrable interaction dresses the initial state. To start with, any quench within the free case g = 0 corresponds to an initial state consisting of pairs. Now a quench to an interacting point g = 0 can be considered perturbatively and the issue now is how the state is dressed by turning on an integrable interaction. Note that such a dressing must map the non-interacting vacuum state |0 to the interacting vacuum state |Ω and the free one-particle states to the asymptotic one-particle states of the interacting theory, which is rather nontrivial. However, there is at least one case in which we know the result of such a dressing: integrable boundary states. Considering only the parity invariant case, let us start from the free field theory with Robin boundary condition for which the boundary state has the form is the Robin reflection factor. When switching to Sinh-Gordon theory (with an integrable boundary condition), K R gets dressed up into [65] where E parametrizes the boundary interaction, and can be considered as a dressed version of λ. We see that in the above case turning on an integrable interaction dresses the state so that the pair structure is preserved. It is plausible that a perturbative proof can be given, similar in spirit to the mechanism of how particle number changing amplitudes cancel at each order of perturbation theory in the Sinh-Gordon model [66]. We expect the pair structure of the dressed state to be the consequence of the pair structure of the starting state and the integrability of the dressing interaction, and to hold in general. In any case, by analogy we expect the initial state to have this pair structure for any quench from some mass and zero interaction to any other mass and any interaction.
We remark that an approximate pair structure is expected to hold for small quenches irrespective of integrability. A quench is considered small when the post-quench energy density is small compared to the natural scale m 2 where m is the mass gap in the post-quench system. In this case the density of particles created after a quench is small, and generally 3 their creation is dominated by pairs separated by a distance larger than the correlation length ξ = m −1 . One expects that the pair creation amplitude K 2 (θ) is small, and furthermore, the amplitude for creating any higher number of particles is well-approximated by the product of amplitudes corresponding to independent pairs, leading to an initial state (4.1) having the form where we shortened K 2 to K.

Exponentiation from extensivity
In the cumulant expansion of a state with such pair structure, all functions K Ψ 2n would contain n δfunctions in order to account for the pairing of momenta. According to the above, such a state would not satisfy the extensivity requirement unless all K Ψ 2n with n > 1 vanish. We therefore conclude that the pair structure and the extensivity requirement constrain the initial state to be of the squeezed coherent form. We note that the above argument is analogous to the one used to show asymptotic exponentiation in [42], although the actual statement and the detailed reasoning are different. Our argument is an application of the formalism developed in [67], extended to the integrable case.
In this special class of initial states, the one test particle equation (3.19) becomes as we can see by substituting into (3.19). Note that based on the Zamolodchikov-Faddeev algebra (2.5) the pair amplitude K satisfies the relation K(θ) = K(−θ)S(2θ).

An Ansatz for the solution of the hierarchy
Let us now consider the initial state on the basis of our heuristic argument for the pair structure. We attributed the preservation of the pair structure to the fact that when switching on the coupling, the state gets dressed by an integrable interaction. In the free case, the state can be expanded as where index 0 of the two-particle state refers to g = 0, and |0 is the free boson vacuum. In the interacting case the initial state takes the form with |Ω denoting the vacuum state of the interacting theory. We remark that the relative sign between (4.14) and (4.15) is consistent with our earlier conventions (1.3) and (3.28). Note that the interaction gives a non-trivial renormalization of the vacuum state and also of the particles. Motivated by the dressing argument of Subsection 4.1.2, let us look for the pair amplitude K in the form where k = m sinh θ and the dressing factor D(k) only depends on the parameters of the post-quench Hamiltonian, i.e. it is independent of the pre-quench mass m 0 . Now consider the limit for the pre-quench mass m 0 → ∞. In this case the free amplitude tends to 1, while the amplitude K is expected to become identical to (2.16) for the integrable Dirichlet boundary state, which means that D(k) = K D (θ). This leads to the proposal of the following solution 4 for the initial state, where This is exactly the Ansatz proposed in our previous work [25] on a slightly different basis, and it was demonstrated to solve the one-particle equation of the hierarchy to a good approximation; note that it has no free parameters to play with. In fact our arguments suggest that its functional form is exact; however, our reasoning does not preclude a coupling constant dependent renormalization of the parameter m 0 /m on which the K f ree factor effectively depends, similarly to λ in the case of the integrable reflection factor (4.9). However, in view of the excellent agreement of our Ansatz with the numerical solution discussed in the next Section, such a renormalization seems unlikely.

Numerical solution of the hierarchy
Following the arguments of the previous section we can now assume that the initial state has the form As we will see below, for a numerical solution for the function K one needs to consider analytic continuation of the equations (4.12) to complex rapidities. These can be derived by deforming the integration contours and taking into account the residues of kinematical poles given by (2.8). An alternative way to obtain them is provided by the finite volume formalism briefly reviewed in Appendix A. To be certain that the equations are correctly computed, we performed the finite volume derivation and then cross-checked the result against (4.12).

Keeping the vacuum and two-particle terms
Keeping the first few terms in the expansion (5.1) and applying a one-particle test state θ|, the integral equation can be written for the Dirichlet case m 0 = ∞, when the operator equation (3.6) simplifies toφ(p)|D = 0 [25]. For the iteration procedure we need to generalize (5.2) to complex test rapidities. As long as Im θ < ε, (5.2) turns out to be valid, but for Im θ > ε, has to be used. This can be derived from (5.2) by analytical continuation using the analytical properties of the form factors. In the case of finite mass quenches (5.4) holds for Im θ < ε and × F φ 5 (θ + iπ, −θ 1 , θ 1 , −θ 2 , θ 2 )K(θ 1 )K(θ 2 ) + ... . for Im θ > ε [25]. The structure of the integral equation is shown in Fig 5.1.  Numerical calculations immediately show that the vacuum and two-particle terms (those that are zeroth and first order in K) in eqns. (5.4)-(5.3) are much larger than the rest, which is expected from the considerations of Subsection 3.4 as they correspond to ∆n = ±1. As a consequence, it makes sense to construct an iterative method based on the truncated version of (5.4)-(5.3) that includes only the first three terms, while the rest are omitted: for the Dirichlet, and for the finite mass quench.

Dirichlet problem
Discussing first the Dirichlet problem, our aim is to numerically calculate the K D (θ) function for real test rapidities, however, the equations contain this function for complex rapidities θ +iε as well. Consequently, to obtain a closed iterative scheme, we also have to plug in complex test rapidities with imaginary part larger than the shift of the contour to the equations, that is, in each iteration we have to calculate two iterative functions. Instead of calculating the iterative functions at real and shifted rapidities, for practical reasons, both of them are calculated at shifted rapidities 0 < ε 1 < ε 2 . At first, is calculated based on (5.6) as the integration contour is shifted with ε 2 , and then K(θ) is calculated based on (5.7) as the integration contour is now shifted with ε 1 . The equations of this iterative scheme derived from (5.6) and (5.7) read (5.10) As can be seen from (5.10), averaging with the previous iterative function is also performed in the scheme.The solution along the real axis can be obtained by the following equation (5.11) The above iteration scheme will be denoted by S2D. The iterations under this scheme rapidly converge, either the functions with shifted rapidities in their arguments or the real rapidity solutions obtained with (5.11) are concerned. This is shown by calculating the function difference of the consecutive iterative functions with the integralˆd for the real and complex rapidity solutions as well. We took the real part of the functions, as the imaginary parts were much smaller and slowly varying. For a growing number of iterations all the three function differences decreased monotonously and rapidly, and we stopped the iterations when the largest of the three function differences crossed a threshold value from above. In Fig. (5.2) the Dirichlet solution can be seen together with the 6th and 7th iterative functions as the iteration was stopped after the 7th run, the largest function difference being 0.0033. The iterative solution is always very close to the expected exact result (2.16) over all the fundamental range 0 ≤ B ≤ 1 of the coupling strength. The small deviation between the numerical and the exact solution can be attributed to the truncation of the infinite integral equation to the vacuum and twoparticle terms. This is supported by the results presented in Subsection 5.2, where it is shown that the iterative solution of the integral equation is much closer to the exact result K D when higher particle number terms are taken into account.

Finite mass quench
For the finite mass problem (5.10) is modified as yielding scheme S2F. The solution along the real axis is obtained by (5.14) For the finite mass quench problem fast convergence is witnessed again, furthermore, the iterative solution appears to be close to the proposed solution (4.19) (Fig. 5.3). Just like in the Dirichlet case, these observations hold for a large regime of the coupling strength B and the quench parameter m 0 . Note that the deviation between the iterative solution and the Ansatz (4.19) is of the same magnitude as in the Dirichlet case, therefore it can safely be attributed to the truncation of the form factor series, which will be strongly confirmed after taking into account the four-particle contribution in Subsection (5.2).

Adding the O(K 2 ) terms
To construct an iteration scheme including the four-particle contributions of (5.4)-(5.3), we can simply add them (5.6)-(5.9). For a closed iteration, however, we now need also an iterative function that is defined for real rapidities. Although it is possible to construct a scheme working with only two iterative functions, the one presented below uses three of them: for K(θ) (k+1) and K(θ) the integration contour is shifted with ε 2 , and for K(θ) (k+1) ε 2 the contour is shifted with ε 1 . Unlike in the equation for K(θ) where the denominator is 1 + S(−2(θ + iε 1 )), in the equation for K(θ) (k+1) , the denominator is 1 + The equations for this scheme S4D read

(5.17)
We omit the equations of the iteration scheme for the finite mass case, as the corresponding finite mass scheme S4F is easily obtained from S4D by plugging the extra E 0 −E E 0 +E type factors. Similarly to schemes S2D and S2F, averaging with the previous iterative functions is present in each iterative step.
This scheme was chosen from other possibilities by observing that it always performed better than all other schemes we tried. Unlike S2D and S2F and just like other schemes including the four-particle terms, S4D and S4F iterations are unstable. To overcome this issue we stopped the iteration, when the difference between two consecutive iterative functions is the smallest. The function difference was measured by the (5.12) as, similarly to the previous schemes, the difference between the imaginary parts tends to be much smaller and slowly varying.

The three-particle condition
For a numerical verification of the exponential form (1.3) of the initial state, it is necessary to consider higher members of the hierarchy. This rapidly becomes impractical due to the computational cost of the integrals involving very high particle form factors. In addition, constructing the form of the equation with the integration contours arranged conveniently for numerical evaluation is also quite tedious.
For this reason, we restrict ourselves here to the case of three-particle test states. The three-particle condition is a sum 0 = T 0 + T 2 + T 4 + T 6 + . . . (5.18) with T n denoting the n-particle contribution from initial state (1.3). Since the iterative solution of the one-particle condition is very well described by the Ansatz (4.19), the latter can be used to perform the evaluations. We have computed the first three terms explicitly with the result given in (A.15). For the contribution T 6 , deriving the integral form proved to be so tedious that we decided to estimate its contribution directly evaluating the corresponding term in the finite volume form (A.10) for a large value of the volume (mL = 250).
To verify that (5.18) holds, we computed the sum for several values of the Sinh-Gordon coupling in the fundamental range 0 ≤ B ≤ 1, for different values of the test rapidities at each point, and for several different quenches, parametrized by the mass ratio m 0 /m.
To make certain that the integral form (A.15) of the condition was derived correctly, we supplemented the calculation of T 2 and T 4 by a direct evaluation of the finite volume sum for mL = 250. The finite volume form was always found to agree with the integral form within the numerical precision of the former. Note that for practical evaluation the finite volume sum must be truncated to states below some upper energy cutoff. For T 2 and T 4 it was possible to keep this truncation high enough to achieve better than 3 digits accuracy; however, for T 6 this was not always the case, as we discuss below.
A sample of the resulting data is collected in Appendix B. As a benchmark, we always quote the Dirichlet case m 0 /m = ∞, which can be obtained by omitting the energy factors in square brackets from each term. For the Dirichlet case, the equation must hold exactly when terms T n are included for all n, since then (1.3) corresponds to a boundary state known exactly from reflection factor bootstrap [41,61]. We can see that the terms T 2 and T 4 (corresponding to ∆n = ±1) are always dominant, and typically cancel each other within a few percent. To go further, it is necessary to include T 0 and T 6 , corresponding to ∆n = ±2. In most cases, this improves the cancellation to better than a percent. There are two exceptions to this, however. First, when the sum T 0 + T 2 + T 4 is relatively small, to verify the improvement T 6 would need to be evaluated to a very high precision which is not possible using the finite volume summation. Second, when some of the test rapidities are relatively large, the cut-off necessary to evaluate the finite volume sum for T 6 prevents the evaluation of all the dominant contributions, as these come from regions which cannot be explored within reasonable computer time. However, in all these cases we also see the same deviations for the m 0 /m = ∞ case, for which we known the full equation should hold exactly. Taken altogether, these facts show that the deviations can be explained by the approximations made during numerical evaluation. On the other hand, we have data for a much larger number of couplings and rapidities than shown in Appendix B, and all of them fall in the same pattern.
Summing up, the numerical evaluation of the three-particle member of the hierarchy strongly confirms the exponential form of the initial state.

Conclusions
The present paper continues a program started in [25], which is aimed at constructing the initial state of quantum quenches in integrable quantum field theories in terms of the post-quench eigenstates. Such a construction amounts to the determination of the overlaps of the initial state with the eigenstates of the post-quench Hamiltonian, which is an essential input for computing the temporal evolution and the steady state expectation values of physical observables.
The particular class of quenches considered in this paper was within the Sinh-Gordon model, starting from the ground state of mass m 0 with zero coupling, to a post-quench system with a mass m and a nonzero value of g. The important characteristics of this quench is that the initial state can be specified via an operator condition. This allows the derivation of an infinite hierarchy of integral equations based on the form factor bootstrap for the overlap functions K 2n . Each integral equation can be written as a form factor expansion, consisting of an infinite number of terms.
We then examined the hierarchy, and presented a number of arguments concerning the nature of its solutions. Some of these arguments are likely valid for more general models and initial states, described by similar integral hierarchies. The results we expect to be generally valid are the following: 1. For a free field, the unique solution of the hierarchy is the usual squeezed state obtained from the Bogoliubov transformation, and perturbation theory considerations imply the existence and uniqueness of its solution for the interacting case.
2. The terms in the hierarchy are ordered in magnitude by the difference between the bra and ket state particle numbers ∆n; the larger ∆n, the smaller is the corresponding term. This means that each equation in the hierarchy can be well approximated by a finite truncation, and that higher equations corresponding to test states with more particles probe overlap amplitudes with states containing more particles. For the elementary field φ considered in this paper it is the∆n = ±1 terms that dominate, and a rather good approximation can be obtained by truncating to only two terms.
3. The iterative solution method that was developed on the basis of the aforementioned truncatability property.
In addition, for the particular class of quenches we presented a heuristic argument, dubbed "integrable dressing", that the state only contains pair states. This argument is heuristic in many respects. First, it supposes the construction of a mapping of the free theory eigenstates to the interacting eigenstates; there are general results such as Haag's theorem [68] showing that such a mapping does not exist in mathematically strict sense. This can be avoided by introducing an ultraviolet cut-off and removing it after renormalization, but then it must be shown that the resulting relations really confirm the argument. An equivalent approach is to consider fast but smooth quantum quenches [69] in which the time interval during which the quench is performed tends to zero but kept larger than the inverse energy cut-off. Either way, this is definitely out of the scope of the present paper; however, completing the argument or finding an alternative is definitely worthwhile to consider, especially in view of extending the result to other integrable quenches. The heuristic nature of the pairing argument makes it very important to supplement the analytical considerations by numerical ones. First, the iterative solution confirms to high precision that the proposed factorized Ansatz (4.19) indeed solves the one-particle test state condition, which is the lowest member of the hierarchy. Second, an independent test of the squeezed state form, which includes the pair structure and the exponentiation, was provided by checking the next member of the hierarchy, i.e. the three-particle condition. Note that since the Ansatz contains no free parameters, there is no way to adjust it: it either fails or passes. In the end, the Ansatz passed all the tests we could pose; while these are limited by computing power, they still impose very stringent constraints. This shows that the Ansatz is at least a very good approximation to the solution, and raises the possibility that it is eventually exact; in this regard it could be of interest to work out the integrable dressing argument (cf. Subsection 4.1.2) in more detail for the quench situation.
Our results show that an exponential cut-off regularisation of the state of the form (1.4) cannot provide a consistent description of the initial state (or its asymptotic behaviour at large times) in the limit of large initial mass m 0 . This is because the original amplitude K(θ) of (4.19) decays algebraically as a function of momentum p for large momenta (it behaves as 1/p 2 ) and therefore physical observables with different dependence on momentum scales would exhibit different scaling behaviour in the Ansatz and in the exponentially cut-off state as m 0 → ∞. In other words, fixing the correspondence between τ 0 and m 0 for some test observable would give inconsistent results for other observables. A demonstration of this discrepancy was presented in [25]. The same argument would hold for any other choice of regularisation, different from the exponential cut-off of (1.4), unless it decays exactly as the Ansatz itself.
Interesting open issues include applying the present approach to quenches in other field theoretical models, such as sine-Gordon theory or the O(3) sigma model, that are more relevant for applications in condensed matter.
Another challenge is to provide a less heuristic argument for the pair structure, or to prove the exponentiation of the initial state directly from the hierarchy equations. One possible approach to tackle this problem is to first focus on the Dirichlet case, for which both the exponential form and the amplitude K D are fully determined, in a completely different way, from the functional constraints imposed by the boundary bootstrap [41]. Since the pole structure and asymptotic behaviour of K D are linked to the above constraints, which are sufficient to determine it, one expects that they must also be sufficient (along with the properties of the form factors) in order to verify analytically that it satisfies the hierarchy equations. On the other hand, exponentiation may also be explained by a reduction of the N -particle equation to the one-particle one in a way analogous to the free case, though definitely more elaborate.
It would also be worthwhile to find more efficient numerical methods to solve the hierarchy determining the overlap amplitudes.

A.1 The finite volume regularization
In a finite volume representation the initial state takes the form [70] |B L = |Ω L + I >0 where the summations run over all positive integer quantum numbers for the particles, and the rapidities are related to the quantum numbers by the Bethe-Yang equations: for the two-particle, and for the four-particle term. These are quantization conditions for states that are constrained to have a pair structure, and similar equations can be written for higher particle numbers. We can also define the unconstrained multi-particle states |θ 1 , . . . , θ n L , which (up to corrections that vanish exponentially with the volume) satisfy the Bethe-Yang quantization equations Q (n) The density of unconstrained states in rapidity space is given by the Jacobi determinant [71] ρ n (θ 1 , . . . , θ n ) L = det ∂Q In these notations, the finite volume normalization factors are [70] N 2n (θ 1 , . . . , θ n ) L = ρ 2n (−θ 1 , θ 1 , · · · − θ n , θ n ) L ρ 2n (θ 1 , . . . , θ n ) L .
Let us introduce the short-hand notation In finite volume, we can write the following form for the n-particle equation in the hierarchy are the infinite volume multi-particle form factors of O(p). They can be expressed using the form factors of the elementary field φ(x) as Using these formulas, one can derive the equations satisfied by the state |B in the infinite volume limit by applying the methods developed in [73]. The trick is to rewrite the summation over the quantum number using multi-dimensional residue integrals, and open the contours which necessitates compensating for some additional singularities by subtracting their residues. In the Supplementary Material of the previous work [25] we have shown an explicit example of this sort of calculation; therefore here we refrain from going into more details.
For one-particle test states, the resulting equation up to (and including) four-particle terms from the initial state can be written as (5.4) valid as long as Im θ < ε.
For Im θ > ε, the equation is changed by kinematic poles of the form factors F φ 3 and F φ 5 crossing the contours of integrations; the appropriate analytic continuation of the equation reads (5.5). In the limit m 0 → ∞, the energy terms E(θ) can be dropped and the equation divided through by E 0 (p), leading to the simplified form (5.2) for Im θ < ε, and (5.3) for Im θ > ε. These express the condition φ(x)|D = 0 (A.14) satisfied by the Dirichlet state.
For three-particle test states, the result is with T n denoting the n-particle contribution from the state |B : , (A.18) and p = m sinh θ 1 + m sinh θ 2 + m sinh θ 3 , valid as long as Im θ i < ε. For other complex values of the test rapidities it can be continued analytically similarly to the one-particle equation; the condition satisfied by the Dirichlet function K D can be obtained by dropping the terms containing combinations of E 0 and E in the square brackets.

A.2 Comparing to the infinite volume formalism
The equation hierarchy can be obtained directly from substituting the infinite volume matrix element (2.11) into (3.16). Considering the case of a one-particle test state, from (4.12) one obtains In the finite volume formula it is necessary to take the form (5.5) valid for Imθ > ε to have the same ordering of the imaginary parts between the unprimed and primed rapidity variables as in (A. 19) above. Shifting back the contours to the real axis, and absorbing S-matrix factors by reordering the rapidity variables in the corresponding form factor gives × F φ 5 (θ + iπ + iε, −θ 1 , θ 1 , −θ 2 , θ 2 )K(θ 1 )K(θ 2 ) + ... .
Due to the +i0 shifts in the unprimed rapidities, the −i0 shifts in (A.19) can be eliminated, making the two equations identical. Similar identity can be demonstrated for the three-test particle condition; as it contains no essential novelty compared to the one-particle case, for the sake of brevity we omit the details here.
B Tables for the numerical verification of the three-particle condition The tables in this Section contain a sample of numerical data obtained from numerical evaluation of the three-particle condition (A.15). The first three terms T 0 , T 2 and T 4 were evaluated using the integral formulae. The first line labeled "sum" gives the sum of these three terms, and verifies how precisely the condition is satisfied in this truncation. For the evaluation of T 6 it proved practical to use the finite volume sum form (A.11). The second line labeled "sum" gives the value of the three-particle condition once the computed result for T 6 is included as well.