Boundary State in an Integrable Quantum Field Theory Out of Equilibrium

We study a quantum quench of the mass and the interaction in the Sinh-Gordon model starting from a large initial mass and zero initial coupling. Our focus is on the determination of the expansion of the initial state in terms of post-quench excitations. We argue that the large energy profile of the involved excitations can be relevant for the late time behaviour of the system and common regularization schemes are unreliable. We therefore proceed in determining the initial state by first principles expanding it in a systematic and controllable fashion on the basis of the asymptotic states. Our results show that, for the special limit of pre-quench parameters we consider, it assumes a squeezed state form that has been shown to evolve so as to exhibit the equilibrium behaviour predicted by the Generalized Gibbs Ensemble.

We study a quantum quench of the mass and the interaction in the Sinh-Gordon model starting from a large initial mass and zero initial coupling. Our focus is on the determination of the expansion of the initial state in terms of post-quench excitations. We argue that the large energy profile of the involved excitations can be relevant for the late time behaviour of the system and common regularization schemes are unreliable. We therefore proceed in determining the initial state by first principles expanding it in a systematic and controllable fashion on the basis of the asymptotic states. Our results show that, for the special limit of pre-quench parameters we consider, it assumes a squeezed state form that has been shown to evolve so as to exhibit the equilibrium behaviour predicted by the Generalized Gibbs Ensemble. Introduction. Research in non-equilibrium processes of Quantum Field Theory (QFT) and their statistical mechanical properties constitutes a fast developing and widely applicable area of theoretical physics. The correct understanding of phenomena out of equilibrium not only plays a crucial role for our knowledge about as diverse topics as cosmology, heavy-ion collision experiments and cold atom systems (see, for instance, [1][2][3] and references therein), but it also poses purely theoretical questions in the subject of QFT itself. This is particularly true in the case of (1+1) dimensional integrable QFTs, i.e. systems which have an infinite number of local integrals of motion [4]: in this case, connecting far-from-equilibrium dynamics at early times with the approach to equilibrium at late times may be a true challenge for the theory. In particular, the experimental evidence of lack of thermalization in a 1d system of bosons with pointlike interactions [5] (a system described by a special limit of the integrable QFT of the Sinh-Gordon model [6,7]) led to the conjecture that quantum integrable systems exhibit equilibration to a Generalized Gibbs Ensemble (GGE) rather than the usual Gibbs ensemble of thermal equilibrium [8].
The GGE is associated to a density matrix that involves all local integrals of motion Q i of the integrable model, including the Hamiltonian. The validity of the GGE has been verified with a variety of different approaches and settings in many systems which can be mapped to free boson or fermion systems, even though such mappings are often highly non-trivial [9][10][11][12][13][14][15][16][17][18][19]. For genuine interacting integrable QFT there have been so far only a few studies [20][21][22][23][24][25][26][27][28][29][30], and presently the most general result concerns the time-average of one-point functions of local operators, for which it was shown in [28] that their values can indeed be recovered by the GGE average.
In many non-equilibrium situations of a QFT, the future evolution of the system is entirely encoded into the specification of the initial state |B , also called boundary state. This is what happens, for instance, in the global Quantum Quench (QQ) process, where a parameter of the Hamiltonian is abruptly changed at t = t 0 and the role of the boundary state is played by the ground state of the pre-quench Hamiltonian.
The subject of this paper is the theoretical investigation of the boundary state closely related to Dirichlet boundary conditions of an interacting integrable QFT, in particular its determination according to basic principles. An important issue of our analysis will concern the proper treatment of the ultraviolet unbounded behaviour originally present in the expression of |B , a task that will lead us to a non-trivial set of equations involving the exact matrix elements of the field φ(x) on the asymptotic states. For convenience we will present our main results through the simplest representative of these interacting theories, i.e. the Sinh-Gordon model based on the bosonic field φ(x), the generalisation to more complicated integrable QFT being straightforward.
2. Boundary States. In global QQ, |B is the ground state of the pre-quench Hamiltonian and an important step for solving the subsequent out-of-equilibrium dynamics is to express this state in terms of the operators that create the particle excitations of the post-quench Hamiltonian. A familiar and simple example of this procedure, which will be important for our future considerations, is the quench process of the mass term m 0 → m of a free bosonic massive QFT [10]: in this case, introducing with E(p) = p 2 + m 2 and E 0 (p) = p 2 + m 2 0 , and denoting by (A 0 (p), A † 0 (p)) and (A(p), A † (p)) the pre/postquench annihilation and creation operators, these two sets are related by a Bogoliubov transformation In this example, the boundary state is identified by the condition A 0 (p)|B = 0, which can be expressed in terms of the post-quench operators as The solution of this equation provides the sought expression of the boundary state in terms of the post-quench operators where |Ω is the ground state of the post-quench Hamiltonian and From the point of view of the post-quench system, the boundary state is then an infinite superposition of pairs of equal and opposite momentum, each of them weighted with the amplitude K free (p) (see Figure 1). Figure 1: With respect to the post-quench Hamiltonian, the boundary state |B appears as a coherent superposition of an infinite number of pairs of particles with equal and opposite momentum.
In the case of interacting integrable QFT (for simplicity we are considering integrable QFT with only one type of particle excitations), a generalization of this class of boundary states consisting of an infinite number of pairs of equal and opposite momentum is given by the general where Z † (θ) are creation operators of the post-quench Hamiltonian, |Ω is its ground state and the variable θ conveniently parameterizes the dispersion relation of the particle excitations of mass m, given by E = m cosh θ and p = m sinh θ. The operators Z(θ) and Z † (θ) provide a complete basis of the Hilbert space of the post-quench integrable QFT and satisfy the Zamolodchikov-Faddeev algebra that involves the exact two-body S-matrix S(θ 1 − θ 2 ), function of the rapidity differences. Boundary states of the exponential form (5) have been considered in [20,32], where it was shown that they automatically lead to equilibration of one-point observables according to the GGE (in agreement with the general analysis done in [28]). Such a class of exponential states includes the important examples of integrable boundary states, i.e. boundary states that respect the integrability of the bulk [33], like the Dirichlet states that can be prepared by a special QQ where the mass parameter m 0 of the pre-quench Hamiltonian is sent to infinity. While it is presently not known whether a general QQ in an integrable QFT leads to an exponential state (5), we will argue that this may be the case for a natural class of quench processes in the Sinh-Gordon model. Before facing this aspect of the problem, let us discuss another important issue related to boundary states.
3. The problem of ultraviolet behaviour. The Dirichlet state |D , like all other known integrable boundary states, suffers from an ultraviolet unbounded behaviour. For example, taking literally the expression that comes out in the limit m 0 → ∞ of the free bosonic case (3), one realizes that the corresponding amplitude K D (p), associated to an idealized Dirichlet state, is a constant and does not decay sufficiently fast for large momenta. This is easily understood since excitations on the initial state are cut-off by m 0 which in this case is taken to be infinite. Similar unbounded behaviour is present in any other integrable QFT and, as a result, gives rise to divergent expectation values.
When the post-quench Hamiltonian is a Conformal Field Theory, a cure of this problem was proposed by Cardy and Calabrese [10,11], who assumed a large but not infinite m 0 and made use of the concept of 'extrapolation length' τ 0 , already known from boundary Renormalization Group (RG) theory, to account for the difference of the actual initial state from the idealized Dirichlet state. The parameter τ 0 is a 'small' regularization param-eter that plays the role of the inverse of an exponential cut-off and depends on the initial parameters: chosen to be of the order 1/m 0 , it turns out to give indeed a good approximation of the initial state when the post-quench system is critical.
In the case of massive post-quench theories, the obvious generalisation of the above idea is to replace K D (θ) by K D (θ)e −2E(θ)τ0 and to postulate that τ 0 is still of order 1/m 0 as in the conformal case. In order to check the validity of this assumption and estimate a suitable value for τ 0 , one may choose an observable and equate its expectation values in the exact and in the approximate initial state. If the τ 0 regularization were consistent, this estimate should be independent of the choice of observable under consideration. However this turns out not to be true, as shown in detail in the Supplementary Material. For instance, for a mass quench in free bosonic theory, choosing as observable the operator φ 2 (x), one arrives at the scaling relation where γ is the Euler-Mascheroni constant. Although this goes as in the conformal case, τ 0 ∼ 1/m 0 [10,11], the prefactor is however different. A similar scaling law, but with another prefactor, is obtained choosing as observable the Hamiltonian: in this case one arrives at Even though the two scaling laws (7) and (8) are very close numerically, the impossibility to arrive to a universal expression of τ 0 is nevertheless a flaw of the present regularization scheme. This discrepancy can be interpreted as an indication that the effect of higher energy excitations present in the initial state cannot be incorporated in an appropriate and unique definition of an energy cut-off: observables that weigh differently the effect of low and high energy excitations can then reveal different ultraviolet behaviour.
A way out of these difficulties is to assume τ 0 to be not a constant but a quantity that depends on p itself. In fact, such a dependence is perfectly justified from the point of view of boundary RG, according to which the actual boundary state may be constructed involving any boundary irrelevant operator, as recently discussed in [35,36]. In this approach, the introduction of the extrapolation length τ 0 amounts to a perturbation of the boundary state generated by the bulk Hamiltonian i.e. the state becomes e −Hτ0 |D . In addition to the latter, one must in general introduce a different τ 0 for each bulk conserved charge Q s , which are indeed boundary irrelevant operators. This would lead to a regularised initial state obtained by e − s Qsτ0,s |D , which turns out to be still of the form (5) but with a τ 0 that is a momentum-dependent function. This is because all charges Q s of an IFT can be put in the form dθe sθ Z † (θ)Z(θ) [4]. However the problem of how to determine the suitable function τ 0 (θ) or equivalently K(θ) remains.
In the following we will study a QQ in the Sinh-Gordon (shG) model in which we start from a large but not infinite mass m 0 and use the exponential form (5) as an Ansatz, providing a series of arguments for such a choice. We then derive, from first principles, a sequence of integral equations that must be satisfied by the function K(θ) and propose a solution based on an analytical approximation which we verify numerically with a high level of accuracy. This provides a posteriori a non-trivial check of the validity of our initial Ansatz (5).

The Sinh-Gordon Model. The shG Hamiltonian is
is a real scalar field, π(x) its conjugate momentum, µ the mass and g the coupling constant. In this integrable field theory there is only one type of particle with physical renormalized mass m given by m 2 = µ 2 sin απ/απ where α is the dimensionless renormalized coupling constant α = g 2 /(8π + g 2 ). Particle scattering is fully determined by the two-particle S-matrix given by [5,37] where θ is the rapidity difference between the particles.
Let us consider a QQ in the shG model starting from a large initial mass m 0 and, for reasons that become clear soon, zero interaction α 0 = 0: the final quantities are finite values of m and α. Such a quench may be regarded as made of a sequence of processes: an initial quench of the mass in free bosonic theory, m 0 → m, swiftly followed by a switching on of the coupling, α 0 → α.
To determine the boundary state |B for this QQ in terms of the post-quench Hamiltonian, let us use the condition that |B is annihilated by the pre-quench annihilation operator Z 0 (p) [39]. The choice of the initial coupling value α 0 = 0 is particularly convenient because in this case Z 0 (p) is just the annihilation operator of the free bosonic theory, easily expressible in terms of the physical field operator φ(x) and its conjugate momentum π(x) as Since in a QFT we have π =φ = −i[φ, H], we arrive to the following equation To make progress in the solution of this equation, let us first expand the state in the post-quench basis in the most general way where is the post-quench eigenstate containing s particles with rapidities θ 1 , ..., θ s and p(θ r ) = m sinh θ r is the momentum corresponding to rapidity θ r .
Additional constraints on |B come by exploiting the symmetries of the quench process and the boundary state. Since this is the ground state of the pre-quench free Hamiltonian, it is invariant under parity and translation invariance. Moreover both symmetries are preserved by the quench process: hence, for parity reason, the sum runs over even integer numbers of particles only, while, for translation invariance, each term in the sum has zero total momentum, as ensured by the δ-function.
Applying suitable test states on the left of (12), we can derive integral equations satisfied by the amplitudesK s of the excitations present in |B . However our investigation drastically simplifies if we assume that the state is of the exponential form (5). If we apply first the assumption that the state consists only of pairs of particles with opposite rapidities, we have where, due to the algebra (6), the amplitudes K s satisfy the properties Assuming further that the state is of the more special form (5), the amplitudes K s are related to each other by The plausibility of such an Ansatz comes from a series of reasons: first of all, from the vanishing of the expectation values on the state |B of all infinite conserved charges Q − a (a = 1, 3, 5, . . .) of the Sinh-Gordon model which are odd under parity transformation. Indeed, if P is the parity operator which is conserved in the quench process, then P Q − a P = −Q − a and since |B is an even state, P |B = |B . Therefore Since on the asymptotic states such charges act as Q − a |θ 1 , ..., θ n = n k=1 sinh(aθ k ) |θ 1 , ..., θ n a pair-wise structure of the boundary state automatically guarantees the vanishing of the expectation values on the state |B . Secondly, imagine to realize the overall quench in terms of a sequence of quenches, the first QQ 1 in which we change only the mass (at α 0 = 0), the second QQ 2 in which we switch on the coupling. After QQ 1 , the resulting boundary state is |B free given in (3), which is made of pairs of equal and opposite particles created by the free operators Z † 0 (p) with mass m. After QQ 2 , where we have switched on the coupling constant α, the infinite number of pairs present in |B free start interacting between them. However the interaction provided by the Sinh-Gordon model cannot create or destroy particles since it is integrable and when the particles cross each other, they just experience a time-delay dictated by the elastic S-matrix given in (10). It is therefore conceivable that the only effect of interaction is to "dress" both the free particle amplitude K free (θ) → K(θ) and the free creation operators Z † 0 (p) → Z † (p), preserving though the pair-wise structure of the boundary state.
Assuming the validity of the pair-wise structure of the initial state and the exponentiation of the amplitudes, i.e. assuming the form (5), let us start our analysis from the limit m 0 → ∞ which corresponds to the Dirichlet state |D satisfying the conditioñ Such a boundary state is known to be of the exponential form (5) with amplitude K D (θ) given by [40] K D (θ) = i tanh (θ/2) Such a known case provides a non-trivial check of the approach we are going to propose. Indeed, if we now take as test state an arbitrary 1-particle excitation θ| applied to the left of (16), substitute (5) and expand, we obtain in this way the following integral equation that must be satisfied by the amplitude where F n ({θ j }) are the matrix elements (the so-called Form Factors) of the field φ defined by F n ({θ j }) ≡ 0|φ|{θ j } . In the derivation of (18) we have exploited both the crossing symmetry and the analytical properties of QFT [3,4] which have allowed us to write the matrix elements θ|φ|{−θ i , +θ i } in terms of the Form Factors above. The exact expressions of the Form Factors of the Sinh-Gordon model were computed in [5,37] (for convenience, their exact expressions can be found in the Supplementary Material). Moreover since the numerical value of the F n 's decreases with the order n, the series (18) shows a fast convergent behaviour and can be approximated to the desired order of accuracy simply restricting to the lowest terms. There is however a technical issue to take care of: since the Form Factors have poles whenever an in-and an outrapidity coincide, one needs to choose a suitable prescription on how to pass around the poles, in order that the equation above makes sense. This prescription is encoded in the integration contours C i which can be determined by means of a finite volume regularization [6,7,42,43] (details are discussed in the Supplementary Material). Once this prescription is implemented, the first few terms of the series give as a result the equation In Fig. 2 we plot the numerical solution of (19) when we keep its first three and five terms, along with the analytical result (17). As shown in the Figure, the agreement is quite satisfactory even when the series is truncated up to the first three terms and it improves significantly once the next two terms are included. Using multi-particle test states applied on the left of (16), one obtains a series of equations that must all be satisfied by the amplitude K D (θ). More details of such computations will be presented elsewhere [46].
Supported by this positive check for the case m 0 → ∞, let us now address the problem of determining the amplitude K(θ) in the case of large but finite m 0 . The equation that defines the initial state is now (12). Considering a 1-particle test state as before and substituting (5), we find, after some algebra, that the new equation is which, after a truncation of the series to the same order as before, becomes One way to obtain the solution of this equation is to notice that, for a smooth function K, the integral of the first line is dominated by the contribution of the kinematical poles of the Form Factor at θ = ±θ . At these poles, the prefactor (E 0 (θ) − E(θ) + 2E(θ ))/(E 0 (θ) + E(θ)) of the integration kernel becomes equal to unit. This suggests the approximate solution since then the first line of (21) becomes approximately the same as the first line of (19). The second and third lines contribute only small corrections to the solution. The correctness of our approximate solution can be verified numerically. Fig. 3 shows a typical plot of a numerical solution of (21) truncated to the 3rd or 5th terms, along with the Ansatz (22), for some values of the ratio m 0 /m and interaction α. The agreement is quite satisfactory, even when we include the contribution of the second and third lines. Further comparative plots for a wide range of parameter values will be presented elsewhere [46].
Obviously, the proposed solution (22) is expected to be more accurate when the parameters m 0 , m and α are such that the domination of the poles is more prominent and the higher order terms of the series give smaller contributions. The second condition is satisfied for example when both masses m 0 and m are large or when the interaction α is small.
In analogy to the Dirichlet state case, using multiparticle test states we can derive a series of equations that must be satisfied by the amplitudes K s of the initial state in the form (14). Based on the same combination of arguments used above (truncation of the form factor series and pole dominance of the integrals), it is possible to show that also these equations reduce approximately to the ones corresponding to the Dirichlet case when the K s are chosen to be K s (θ 1 , ..., θ s ) = 1/s! s r=1 K 1 (θ r ) with K 1 (θ) ≈ K D (θ)K free (θ). In this way the exponential form of the Dirichlet state leads also to approximate exponentiation of the QQ initial state.

Observables at large times.
According to the analysis done above, a first order approximation of the initial state |B for the QQ under consideration is given by the exponential form (5) with amplitude K(θ) given by (22). This decays for large momenta as a power law (∼ p −2 ) and ensures a smooth ultraviolet behaviour through a momentum dependent τ 0 -regularization. An initial state of this form belongs to the class studied in [20,28,32] and therefore at least the one-point functions of local observables equilibrate according to the GGE: their long time values are given by θ n , ..., θ 1 |O(x)|θ 1 , ..., θ n c (23) whereK is given by the solution of the generalised Ther-modynamic Bethe Ansatz (TBA) equation 24) where ϕ(θ) ≡ −i d(log S(θ))/dθ, as explained in [20]. From the above equations we can calculate numerically, for instance, the GGE prediction for the operator exp(kφ) from which one can also derive all field fluctuations φ 2n by differentiation with respect to k at k = 0. In Fig. 4, the three curves represent three successive partial sums of the series (23). The convergence of this series is particularly fast near the point k = 0, even though to compute the higher moments φ 2n with sufficient accuracy one needs to employ more terms of the series.  6. Conclusions. In this paper we have studied a QQ of the mass and coupling constant in the Sinh-Gordon model, in the special case of a large initial mass and zero initial interaction. We have seen that the ultraviolet regularization of this state is a non-trivial and physically relevant problem. This has led us to develop a systematic method to determine the expansion of the boundary state in the post-quench basis: this consists in solving integral equations for the excitation amplitudes K s which involve the finite-volume prescription of the exact Form Factors of the elementary field. Assuming that the boundary state is of the exponential form (5), we have obtained a first but quite accurate approximation of the solution by truncating the related Form Factor series. The proposed solution is used to derive the large time behaviour of observables.
The fact that the large energy behaviour of excitation amplitudes in the initial state is relevant for the calculation of physical observables at large times, means that models that are effectively equivalent as far as their ground state or thermal equilibrium properties are concerned, may not be equivalent out-of-equilibrium. More generally, we conclude that RG methods and concepts that are valid at equilibrium cannot always be applied directly to out-of-equilibrium problems.
It would be quite interesting to extend the analysis done in this paper to other QQ protocols and determine the relevant amplitudes K s of the corresponding boundary state from first principles.
Acknowledgements. Spyros Sotiriadis acknowledges financial support by SISSA -International School for Advanced Studies under the "Young SISSA Scientists Research Projects" scheme 2011-2012 and by the ERC under Starting Grant 279391 EDEQS. GT was partially supported by a Hungarian Academy of Sciences "Momentum" grant LP2012-50/2012. The work of GM is supported by the IRSES grants QICFT. This work was also supported by the CNR-MTA joint project "Nonperturbative field theory and strongly correlated systems". We would also like to thank the Max Planck Institute for the Physics of Complex Systems (Dresden, Germany) for their hospitality during the international workshop QSOE'13.

FORM FACTORS OF THE SINH-GORDON MODEL
The Sinh-Gordon theory is defined by the action It is the simplest example of Lagrangian integrable field theory and it is invariant under the Z 2 symmetry φ → −φ.
Parameterising the dispersion relations in terms of the rapidity β, p 0 = m cosh β , p 1 = m sinh β , its two-particle S-matrix is given by [1] where β = β 1 − β 2 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.
The form factors (FF) F O n of the Sinh-Gordon model are matrix elements of a generic local operator O between the vacuum and a set of n particle asymptotic states Based on general properties of a QFT (as unitarity, analyticity and locality), the form factor bootstrap approach leads to a system of linear and recursive equations for the matrix elements F O n [2,3] F n (β 1 , . . . , β i , β i+1 , . . . , β n ) = F n (β 1 , . . . , β i+1 , β i , . . . , β n ) S(β i − β i+1 ) , F n (β 1 + 2πi, . . . , β n−1 , β n ) = F n (β 2 , . . . , β n−1 , β n , β 1 ) , For an operator O(x) of spin s, relativistic invariance implies In the Sinh-Gordon model a convenient parameterization of the n-particle FF, which takes into account their kinematical poles, is given by [4] F n (β 1 , . . . , β n ) = H n Q n (x 1 , . . . , where x i ≡ e βi and β ij = β i − β j . F min (β) is an analytic function given by that satisfies the functional equations It has a simple zero at the threshold β = 0 (since S(0, B) = −1) and no poles in the physical strip 0 ≤ Im β ≤ π, with an asymptotic behaviour given by H n are normalization constants, which can be conveniently chosen as with and H 1 , H 2 are two independent parameters. Finally, the functions Q n (x 1 , . . . , x n ) are symmetric polynomials in the variables x i , which have to be fixed by the recursion equations satisfied by the form factors. For FF of spinless operators, their total degree is equal to n(n − 1)/2. On the other hand, the partial degree of Q n in each variable x i is fixed by the nature and the asymptotic behaviour of the operator O which is considered.
Exploiting the parameterization (31) together with the functional equations (29) and (33), the polynomials Q n (x 1 , . . . , x n ) have to satisfy the recursive equations [4] with D n (x, x 1 , . . . , x n ) = n k=1 k m=1,odd We have introduced the symbol [n] defined by (39) and the elementary symmetric polynomials σ (n) k (x 1 , . . . , x n ), given by the generating function A general class of solution has been identified by Koubek and Mussardo [5] and consists of the following expression where M ij (k) is an (n − 1) × (n − 1) matrix with entries These polynomials, which we call elementary solutions, depend on an arbitrary integer k and satisfy One can prove [5] that the whole set of FF of the elementary field φ(x) are given by Q n (0) while those of the trace of the energy-momentum tensor Θ(x) are given by the even polynomials Q 2n (1). One can also prove that these solutions of the FF equations correspond to the matrix elements of the exponential operators e kgφ .

FINITE VOLUME FORMALISM
A formalism that gives the exact quantum form factors to all orders in L −1 was introduced in [6,7]. The finite volume multi-particle states can be denoted where the I k are momentum quantum numbers, ordered as I 1 ≥ · · · ≥ I n by convention. The corresponding energy levels are determined by the Bethe-Yang equations We defining the two-particle phase shift δ(β) as and its derivative will be denoted by Due to unitarity, δ is an odd and ϕ is an even function. We can write Q k (β 1 , . . . ,β n ) = mL sinhβ k + l =k δ(β k −β l ) = 2πI k , k = 1, . . . , n where the quantum numbers I k take integer values. Eqns. (48) must be solved with respect to the particle rapidities β k , where the energy (relative to the finite volume vacuum state) can be computed as n k=1 m coshβ k (49) up to corrections which decay exponentially with L. The density of n-particle states in rapidity space can be calculated as The finite volume behavior of local matrix elements can be given as [6] {I 1 , . . . , I m }|O(0, 0)|{I 1 , . . . , I n } L = F O m+n (β m + iπ, . . . ,β 1 + iπ,β 1 , . . . ,β n ) ρ(β 1 , . . . ,β n )ρ(β 1 , . . . ,β m ) whereβ k (β k ) are the solutions of the Bethe-Yang equations (48) corresponding to the state with the specified quantum numbers I 1 , . . . , I n (I 1 , . . . , I n ) at the given volume L. The above relation is valid provided there are no disconnected terms i.e. the left and the right states do not contain particles with the same rapidity, i.e. the sets β 1 , . . . ,β n and β 1 , . . . ,β m are disjoint. In the course of work we do not need the relations containing disconnected contributions; these were first obtained in [7].
In addition, the proper definition of a finite volume boundary state was derived in [8]. For a state of the form it can be written as where the rapidities {β 1 , . . . ,β k } of the 2k-particle term satisfy the appropriate finite volume quantization conditions The density of these states is given bȳ and the finite volume normalization coefficients are given by

DERIVING THE CONDITIONS SATISFIED BY THE BOUNDARY STATE
The formalism described above is only valid up to finite size corrections that decay exponentially with the volume. However, these have no effect in the L → ∞ limit and so we omit them below.
We demonstrate the method on the simple case of a Dirichlet boundary state, which satisfies Substituting and taking the matrix element with a one-particle state, we obtain Using (53) this can be written as with the quantization conditions Neglecting the higher-order terms, we can then use eqns. (51) and (53) to write where the state densities are Using (56), we can write where we also extended the summation to negative values of rapidity. Now one can trade the sum for a contour integral where the contour C I run around the solution of (62) counter-clockwise, and so we obtain Open the contour into two lines running above and the below the real axis to obtain where the last two integrals are evaluated on contours surrounding β and −β, in counter-clockwise direction and We can easily calculateQ as L → ∞ and Using the third form factor equation from (29), we obtain F 3 (β + iπ, −β , +β ) = (1 − S(−2β )) iF 1 β + β + regular terms and so the residue terms are From (61,62) one can evaluate and so the pole contribution is where we also used the relations Our final result for the integral equation to this order is This calculation can be extended to higher terms in the expansion of the boundary state |B using the multidimensional residue method introduced in [9]. One can also use any other multi-particle test state in place of the one-particle state applied above. Computation of these terms is lengthy and tedious, but presents no further difficulties.

ESTIMATION OF τ0 IN A FREE BOSONIC MODEL
We consider a QQ of the mass in a relativistic free bosonic model. For this problem, the actual initial state is known exactly thus allowing the estimation of τ 0 . The initial state is where E p = p 2 + m 2 is the energy dispersion relation (with E 0p corresponding to mass m 0 ) and A † (p) is the post-quench bosonic creation operator. For m 0 → ∞ this state tends to the Dirichlet state which therefore is |D ∼ exp − ∞ 0 dp 2π A † (−p)A † (p) |0 and the approximation we want to use is In order to find an appropriate value for τ 0 we should compare the expectation value of an observable O(φ) in the exact and in the approximate state and require that they be equal. One choice of such an observable is the field fluctuations φ 2 . We therefore demand the condition Ψ|φ 2 |Ψ = Ψ |φ 2 |Ψ which reads leading for large m 0 (small τ 0 ) to the following scaling relation where γ is the Euler-Mascheroni constant. Notice that the last expression is consistent with the conformal limit m → 0 [10,11] where we know that τ 0 ∼ 1/m 0 , however the numerical prefactor is not the same.
It is easy to see that the estimation of τ 0 is not independent of the choice of observable under consideration, not even in the large m 0 limit where τ 0 should be uniquely defined if it is a good phenomenological parameter. Indeed let us repeat the estimation of τ 0 using a different observable, in order to test if we obtain the same value. Since the states we consider are gaussian, all higher moments of the field fluctuations φ 2n are simply powers of the variance φ 2n = (2n − 1)!! φ 2 n and therefore by fixing τ 0 from the variance as above, we fix all higher moments too. Therefore using any of these higher moments to obtain τ 0 would automatically result in the same estimate. Another observable of interest is the energy density. By requiring that Ψ|H|Ψ = Ψ |H|Ψ we have which finally gives This estimate differs from (82) in the numerical factor by about 3%. Fig. 5 shows τ 0 m 0 as a function of m/m 0 as estimated based on the field fluctuations φ 2 and the energy density H by numerical evaluation of the corresponding equations. The dependence of τ 0 on the observable is rather expected and can be explained by that fact that a QQ creates excitations of arbitrarily high energy and therefore cannot be analysed by means of the low-energy physics only. Perhaps a more serious flaw comes up when we compare the exact with the approximate large time asymptotic value of an observable, as derived based on the value τ 0 estimated from this same observable. For the field fluctuations our approximation turns out to give lim t→∞ Ψ |φ 2 (t)|Ψ (R) = π/8mτ 0 (the subscript (R) stands for "renormalised", i.e. after mass renormalisation has been applied) whereas the exact value is lim t→∞ Ψ|φ 2 (t)|Ψ (R) = πm 0 /8m, always in the limit of large m 0 , small τ 0 . The two results would be equal to each other if τ 0 was simply equal to m −1 0 instead of (82). Therefore our estimate (82) based on the comparison of the initial field fluctuations does not reproduce the exact results consistently as far as numerical factors are concerned.
As suggested in the introduction, in order to overcome these failures of the τ 0 regularisation, one may assume a momentum dependent τ 0 so that (79) reproduces exactly (80), i.e. by defining Notice that for large m 0 , the above equation gives τ 0 ∼ 1/E 0p ∼ 1/m 0 as expected.

NUMERICAL SOLUTION OF THE INTEGRAL EQUATIONS (19) AND (21) GIVEN IN THE TEXT
In order to numerically solve the integral equations (19) and (21) we employ the iterative method. Truncated after the first 5 terms, the equations are of the form K(θ) = K f (θ) + dθ G(θ, θ )K(θ ) + dθ 1 dθ 2 G(θ, θ 1 , θ 2 )K(θ 1 )K(θ 2 ) which are non-linear integral equations. If we further omit the 4th and 5th terms, the equations reduce to linear inhomogeneous Fredholm integral equations of the 2nd kind. In the iterative method we start by using a guess or approximation of the solution K (0) (θ) (for example K (0) (θ) = K f (θ) or K (0) (θ) = K D (θ) in the case of (19)) and calculate the first correction We proceed similarly to calculate the second correction K (2) (θ) from K (1) (θ) and continue iterating this step until the new correction is sufficiently close to the previous one. The convergence criterion can be chosen to be that the difference of two successive corrections is everywhere smaller than some chosen tolerance. To ensure or make convergence faster we can mix every new correction with the previous one using suitable weights.
A complication arises in the above general iterative procedure by the fact that the integration should be performed not along the real θ-axis but along a parallel line shifted by i . This means that an iteration step does not give K (n) along the same line where K (n−1) is known but a shifted one and therefore the iteration breaks. There are several ways to solve this problem. One way is to set to a very small value so that we can assume that the obtained function K (n) (θ) is approximately equal to K (n) (θ + i ) which is needed in the next iteration step. The value of has to be sufficiently small so that the cumulative error after all iterations needed in order to reach convergence is small. An alternative way to apply this approximation is to consider the positions of the kinematical poles shifted according to the -prescription and bring the integration contour along the real axis, so that the iterative method can be applied without problem. In this approximation the integration kernels can be simplified using the analytical properties of the form factors and can be expressed in terms of only one numerically evaluated function: the function F min (β + iπ) for real β which is smooth and easily calculated numerically. Another approach is to consider two copies of the integral equations for two opposite values of so that at each iteration we calculate K (n) (θ − i ) from K (n−1) (θ + i ) and K (n) (θ + i ) from K (n−1) (θ − i ), therefore obtaining all the information needed for the next iteration step. The advantage of this approach is that the value of does not have to be very small; simply smaller than the distance of the closest singularity of the integrands to the real axis. The solution can be found along the real axis after the last step.
Another practical point is that, since the 5th term of the equation requires more computing time due to the double rapidity integral that it contains but is smaller in size than the previous terms, one may choose a less dense discretisation of rapidity to reduce the computing time. In addition one may omit this term from the equation at first, solve it and then improve the solution by including the previously omitted term and solving again. This two-stage iterative method obviously requires less overall computing time.
In Fig. 2 and 3, we have used the first approach for the numerical implementation of the -prescription, i.e. we used a very small value = 0.0001 and ignored the line shifts. As a starting guess we used K D (θ) for equation (19) and our ansatz (22) for equation (21). We iterate the equation first omitting the 4th and 5th terms which we include later. Convergence of the iterative procedure is achieved with a tolerance of the order of 0.05 after only one step almost everywhere; except the neighbourhood of the point θ = 0 where pathological behaviour should be expected due to the fact that the truncation of the form factor series becomes worse at this point.