Finite temperature one-point functions in non-diagonal integrable field theories: the sine-Gordon model

We study the finite-temperature expectation values of exponential fields in the sine-Gordon model. Using finite-volume regularization, we give a low-temperature expansion of such quantities in terms of the connected diagonal matrix elements, for which we provide explicit formulas. For special values of the exponent, computations by other methods are available and used to validate our findings. Our results can also be interpreted as a further support for a previous conjecture about the connection between finite- and infinite-volume form factors valid up to terms exponentially decaying in the volume.


Introduction
Correlation functions and expectation values of operators are important objects in quantum field theory, both from the theoretical and phenomenological point of view. Integrable quantum field theories have numerous applications to condensed matter systems; given that experiments are necessarily conducted at nonzero temperature the construction of finite temperature expectation values and correlation functions in integrable quantum field theories is an interesting problem. Almost fifteen years ago, LeClair and Mussardo [1] put forward a conjecture for both the one-point and the two-point functions of integrable models with diagonal scattering, expressed as a spectral series using exact form factors and the thermodynamic Bethe ansatz. In [2], another approach to finite temperature expectation values in the sine-Gordon model was proposed by Lukyanov, and more recently Negro and Smirnov provided a resummation of the spectral series of the one-point functions, again in the sinh-Gordon model [3].
For generic one-point functions, the LeClair-Mussardo proposal was eventually proven to be valid in [4], using the finite volume form factor formalism introduced in [5,6]. Concerning two-point functions, their proposal is more controversial [7] and probably not valid in its original form. However, the finite volume form factor formalism provides an alternative and systematic method to construct the two-point function. This approach solves the problem faced by earlier studies which could not resolve the issues related to a proper regularization of kinematical singularities of the form factors [8,9]. An early implementation of the finite volume approach for the two-point functions was used to describe finite temperature line shapes and dynamical correlations [10,11]. The full formalism itself was developed in [12,13]. We note that an alternative approach to thermal correlations was developed by Doyon [14,15], however, at present it seems to be confined to the Ising model.
The finite volume form factor methods were recently shown to yield results agreeing with other approaches in non-equilibrium steady state systems [16], and are also relevant in the context of quantum quenches [17,18,19].
Presently, the available results on form factor expansions of thermal correlators in integrable field theory are limited to the case of diagonal scattering. Conversely, much less is known about non-diagonal integrable field theories: this is partly due to the fact that the LeClair-Mussardo expansion, in its original formulation, requires the solution of the thermodynamic Bethe ansatz equations, which are considerably more difficult [20] when the theory is not diagonal. The finite volume formalism independently provides a way to extend the results to non-diagonal scattering, and recently finite volume form factors for nondiagonal scattering were constructed [21,22]. Albeit at present diagonal matrix elements of multi-soliton states are still not completely known, the available results permit the evaluation of the spectral series below the three-particle threshold. In this paper, we take the first step and consider finite-temperature expectation values in the sine-Gordon model, i.e. the one-point functions, which is performed in section 2. In section 3 we construct the connected diagonal matrix elements of the exponential operators, which allows the evaluation of the series for these observables. Exponential operators are particularly useful because they appear in many physical systems in one dimension, in connection with the characterization of lattice models at low temperatures, by passing to a continuous description through an effective bosonic action (see, e.g., [23,24,25]). In addition, exponential operators generate all the normal-ordered powers of the sine-Gordon field, provided it is possible to compute their expectation values with generic parameter in the exponent.
In section 4 we compare the spectral series for the case of the trace of the energy momentum tensor to results that follow from the NLIE approach [26,27,28,29] and find very good agreement. Unfortunately, for reasons discussed towards the end of section 4, we cannot perform an analogous numerical verification of our method for other operators at present. Nevertheless, our present results provide a nontrivial verification of the method and an analytic check of the form of the diagonal matrix elements conjectured in [22] (where this conjecture was tested numerically against TCSA).

One-point functions at finite temperature
The classical action of the sine-Gordon (SG) field theory is: where λ and β are real parameters, of which β is dimensionless and λ determines the mass scale of the model. Classically λ has dimension of mass squared, but in the quantum theory it acquires an anomalous dimension The fundamental excitations of the model are known to be the soliton, with mass m and unit topological charge, and the antisoliton, with equal mass and opposite charge; the exact relation between λ and m has been derived by Zamolodchikov in [30]. In addition to solitons, the spectrum may also contain breathers which are bound states of a soliton and antisoliton; after quantization their spectrum becomes discrete and only a finite number of such states exists. Introducing the parameter ξ = β 2 8π−β 2 , it is possible to distinguish two regimes: a repulsive one ξ > 1, in which only the soliton and the antisoliton are present in the spectrum, and an attractive one ξ < 1, in which ⌊1/ξ⌋ different bound states (breathers), whose mass is are allowed. The scattering matrix between the elementary excitations of the system has been computed in [31]; the non-zero elements of the S-matrix in the solitonic sector are The S-matrix elements involving breathers are diagonal and can also be found in [31]. Continued to Euclidean time τ = −it, the action with periodic boundary conditions in τ describes the model at finite temperature T = 1/R. Note that by swapping the role of Euclidean time and coordinate, the finite temperature/infinite volume action can also be considered to be a zero temperature/finite volume action, and so the one-point functions constructed below also have this dual physical interpretation. The exponential fields e ikβφ is the most interesting class of operators to be studied, both because they serve as a generating function for all the normal-ordered powers of the SG field and in connection with one-dimensional lattice systems, where they commonly emerge as a counterpart of lattice operator via bosonization of the effective low-temperature field theory. For the case k = ±1 the expectation value of the exponential operator is identical to that of the perturbing operator cos βφ, which is in turn related to the trace of the stress-energy Θ tensor through [32] is the scaling dimension of the exponential field at the conformal point. The finite temperature one-point function of the exponential operators is defined by Gibbs average: is the Hamiltonian, and the summation runs over a complete set of energy eigenstates |n with energies E n . In infinite volume, the form factors ..an (θ 1 , . . . , θ N ) = 0|O|θ 1 , . . . , θ N a1...an (2.9) of local operators can be computed exactly using the form factor bootstrap [33,34,35], from which any multi-particle matrix element can be reconstructed using crossing symmetry. Here |θ 1 , . . . , θ N a1...an denotes a multi-particle state composed of particles with species a 1 , . . . , a N and rapidities θ 1 , . . . , θ N . The analytic structure of the form factors is fixed by a set of equations, which are built upon the factorized scattering of the model as input. Local operators of a given model can be defined as towers of solutions of the form factor bootstrap equations [36]. For many integrable models, including sine-Gordon theory, the exact solutions are known [36,37,38], therefore the spectral sum (2.7) can be evaluated in principle. However, due to the singularity structure which arises from the form factor axioms, the diagonal matrix elements of the fields are not well-defined, hence the spectral sum needs to be regularized. The regularization of form factors by using a finite volume setting is a useful technique for dealing with lowtemperature expansions of one-point and two-point functions [6,12,13]. To evaluate the one-point function one can apply the method outlined in [6]; however, a careful extension of the approach is necessary due to the presence of non-diagonal scattering, which can be performed using the recent results in [21,22] on finite volume form factors for non-diagonal scattering.
We recall that in finite volume the rapidities are quantized and the space of multi-particle states can be labeled by momentum quantum numbers I 1 , . . . , I N . We introduce the following notation for them: where the index r enumerates the eigenvectors of the n-particle transfer matrix, which can be written as where θ 1 , . . . , θ N are particle rapidities. Due to factorized scattering, the transfer matrix can be diagonalized simultaneously for all values of θ: We can assume that the wave function amplitudes Ψ (r) are normalized and form a complete basis: these eigenfunctions describe the possible polarizations of the N particle state with rapidities θ 1 , . . . , θ N in the state indexed by the species quantum numbers i 1 . . . i N . The transfer matrix can be diagonalized using the algebraic Bethe ansatz (cf. Appendix A of [39]), which enables one to compute the exact form of eigenvalues t (r) and eigenvectors Ψ (r) . The rapidities of the particles in the state (2.10) can be determined by solving the quantization conditions When considering N rapidities which solve these equations with given quantum numbers I 1 , . . . I N and a specific polarization state r, they will be written with a tilde asθ N . For states containing up to two particles, the only subspace in which the transfer matrix has to be diagonalized is the two-dimensional subspace of states containing one soliton and one antisoliton. The basis of the eigenstates is given by [21] Assuming that the finite volume energy eigenstates are chosen orthonormal, the partition function up to (and including) two-particle contribution expands as in which tildes denote rapidities which are quantized according to the Bethe-Yang equations in finite volume L, the index j = s, a is used to denote the elementary solitonic excitations, and the index b enumerates the breathers. The prime in the summations is a reminder that states with equal rapidities for the same kind of particle are not allowed solutions 1 of (2.14) and are thus excluded. Furthermore, the indexes +, − denote the symmetric (antisymmetric) combination of the neutral soliton-antisoliton states: The finite temperature expectation value can then be written as Following the derivation detailed in [6], this can be expanded as again up to (and including) two-particle contributions. We emphasize that in the above expression we have explicitly subtracted the terms in which two elementary excitations with the same topological charge or two breathers have the same rapidity.
Next we use the relation between the finite and infinite volume form factors (valid up to exponential terms) conjectured in [21,22]. The densities of states in rapidity space, corresponding to (2.14) are and the finite-volume matrix elements are given by [21,22] up to terms that vanish exponentially for large L. In the above formulas, the symmetric evaluation is defined as where the (±, ±)-polarized form factors F O ± are defined by In addition (in fact F k s = F k a due to charge conjugation invariance). Following [6], we can also express these results with the connected part of the diagonal matrix elements which is defined as follows. Consider the form factor in which kinematical (simple) poles appear as the regulators ǫ 1,2 → 0 independently; the connected form factor F k(c) ab (θ 1 , θ 2 ) is defined as the part which is nonsingular in both ǫ 1,2 . Using the same arguments as in [6], it can be easily checked that the symmetric and connected diagonal matrix elements are related by where the function G 0 is the logarithmic derivative of the soliton-soliton scattering phase (2.4): and On the other hand, for the states in which the scattering among the particles is diagonal, such as the breather-breather and the soliton-breather, the finite volume matrix elements are known from [6]: which can be expressed in terms of the connected form factors as with j = s, a for soliton/antisoliton, and b standing for the breather kind, where and the function G jb (θ), is defined analogously to (2.29), but starting from the soliton-breather scattering phase while in turn, the scattering phase between breathers b 1 and b 2 defines the function: Substituting into (2.35), and keeping terms up to two particles: which gives the low-temperature expansion of the expectation value of the local field O in the sine-Gordon theory, up to and including two-particle contributions.
3 Connected diagonal matrix elements of the exponential fields

Form factors of exponential fields
Multi-soliton form factors of exponential operators in the sine-Gordon model were obtained by Lukyanov in [37] exploiting the bosonized form of the Zamolodchikov-Faddeev operators [40]: in which : : denotes the appropriate normal ordering while φ,φ are generalized free fields defined in [40] and the contours Γ ± are specified below. These operators satisfy the algebra of asymptotic soliton/antisoliton creation operators with the scattering matrix (2.3). The contractions of the vertex operators are defined as follows: where ; analytic continuations valid for a wider range of rapidities can be found in [21,41].
Using the above definitions, Lukyanov constructed the multi-soliton form factors in the form is the vacuum expectation value of the field and the integration contours are specified as follows. Calling "principal" poles the singularities of the function W (γ) located at γ = −iπ/2, the contours run above the "principal" singularities of the W functions arising from contraction of a given operator e −iφ with all fields on its right and below the principal poles originating from the contraction with the ones on its left. Accordingly, in the definition (3.2) Γ + (Γ − ) denote the contours which pass above (below) the pole at The tilde in (3.10) refers to the fact that our conventions for the form factors differ from those of Lukyanov's by the relation as noted in [21]. This is due to a difference in the conventions of the form factor equation, which corresponds to a redefinition of relative phases of the multi-particle states [21]. Note that all the form factors need to be normalized by the vacuum expectation value of the exponential field, for which a formula has been derived in [42]. It is useful to write that formula in a way which can be used for any value of k. Given two integers M, and defining ∆ = ∆(β) = β 2 /8π, one has: which can be obtained by exploiting the integral representation of the logarithm of the Euler's gamma function. Finally, it was noted in [39] that to make agreement with TCSA results, a sign had to be inserted in the diagonal matrix elements with an odd number of solitons or antisolitons. Therefore, our diagonal matrix elements can be obtained in the following way: where F k denotes diagonal matrix elements in Lukyanov's conventions [37].

The diagonal matrix elements
As a warm-up and a demonstration of how the analytic continuation is performed, it is useful to write down the two-particle form factor in the repulsive regime: in which the contour Γ +(−) passes above (below) the pole at γ = i π 2 (γ = −i π 2 ) and has been deformed to the real axis in the second line. We also introduced the notation A = −2k − ξ −1 .
Consider now the integral part, which is divergent for k ≥ 1/2. It can be written as the analytic continuation of a Fourier transform to imaginary values z = −iA: where the definition (3.6) has been used. Now the integral is convergent, but we still need to continue the integral to the value z = i (2k + 1/ξ) by adding the poles that are crossed in the contour deformation. The result is then and is now convergent for all real values of k. Note that this part is O(ǫ 0 ). The connected part of the matrix element (3.14) is the total O(ǫ 0 ) contribution, which can be collected as This will be compared with the exact formula below in section 4.2. Note that the functionŴ defined in (3.8) is regular in the point −iπ/2 and its derivative can be computed straightforwardly from the definition. Let us now write explicitly the regularized diagonal four particle form factor: where we again denote A = −2k − ξ −1 . The contraction in the last line contains one factor which is independent of the integrated variables and reads: with the usual notation ǫ 21 = ǫ 2 − ǫ 1 ; on the other hand, the integral parts are as follows with the contours depicted in figure 3.1. Moreover, we have Again, the contraction in the last line contains one factor which is independent of the integrated variables and reads: on the other hand, the integral parts are as follows with the contours depicted in figures 3.3 and 3.2. The above formulas for the regularized diagonal form factor contain a divergent integral whenever k ≥ 1/2. It is therefore necessary to find a suitable and meaningful definition for these quantities. We refer the interested reader to Appendix B, where explicit formulas are provided.
4 The trace of the stress-energy tensor

The NLIE prediction
As shown by Zamolodchikov [32], one can compute the expectation value of the trace of the stress energy tensor for any volume from the exact ground state energy. In order to compute the latter, for the sine-Gordon theory, we can make use of the nonlinear integral equation (NLIE) [26,27,28]: where η < π min(1, ξ) and the ground state energy of the theory in finite volume R can be calculated using It satisfies E(R) → 0 as R → ∞ (4.4) In the repulsive regime, one can continue analytically to η = π. Introducing the functions [43] ǫ(θ) = −iZ(θ + iπ) the ground state energy can be written as: The functions ǫ,ǭ are analogous to the pseudoenergies of the thermodynamic Bethe ansatz approach and the two terms can be thought of as resulting from the soliton-antisoliton doublet. However, in contrast to TBA pseudoenergies, they are complex valued functions. From (4.1) it can be deduced that they satisfy the equations where Note that the function G 1 has a pole at θ = 0; on the other hand, in the calculation of all physical quantities, this is of no consequence. The two quantities ε andε are related by complex conjugation, as well as their derivatives below. In fact, one can derive with respect to the volume and obtain where the complex quantities f (θ) = 1 1+e ε(θ) andf (θ) = 1 1+eε (θ) have been used. On the other hand, the derivative of the first of (4.7) with respect to the rapidity, denoted by a prime, is obtained as follows It is then a simple matter to differentiate (4.6) and obtain The subtraction of −1 is related to the asymptotic property (4.4), which means that the ground state energy computed from the NLIE has the bulk term subtracted. For the vacuum expectation value derived from it, this entails that we obtain the finite size corrections to the infinite volume value (2.6,3.12). For later convenience, we also normalized the expectation value by its infinite volume limit.

Connected form factors for the trace of the stress-energy tensor
Diagonal matrix elements of the trace of the stress-energy tensor can be computed as outlined in [1], which is reviewed in Appendix B. In particular, the diagonal one-particle form factor is given by Comparing this to the result (3.17) gives a first check of the regularization method, which is shown in Figure 4.1. Moreover, the connected diagonal two-particle form factors are also explicitly known (see Appendix A):  Figure 4.1: One-particle diagonal matrix element F Θ s = F Θ a of the trace of the stress-energy tensor: comparison between the regularization procedure outlined in Section 3.2 (dots) against the analytic formula obtained from Lukyanov's form factor (solid line). The normalization through the operator vacuum expectation value is not included in this plot.
We checked (for values 1/2 < ξ < 2, where we need to take into account only principal poles) that the connected form factors as obtained in the form (B.20,B.23) from Lukyanov's expressions (as given in Appendix B) agree with these (see figure 4.2). In addition to that, the two can be compared with the results from the regularization procedure of [41], again resulting in agreement among the results, thus providing a threefold check. For future reference, we also write here the soliton-breather and the two-breather connected diagonal matrix element of the trace of the stress-energy tensor:

Comparing the NLIE to the series
Let us start with an analytical comparison in the repulsive regime ξ > 1. It is easy to see that expanding the exact result (4.11)to second order in e −mR of produces exactly the terms of series (2.35) up to this order.
First we need to expand the expression (4.7) and the relative (complex) filling factor, which leads to: while the one forε andf are obtained by the substitution G 1 →Ḡ 1 . The expansion of the NLIE result gives: Using expressions (4.12,4.13) for the connected form factors of the trace of the energy-momentum tensor, one can easily recognize that this is indeed identical to the solitonic terms in (2.35), once the normalization of the two-particle form factor according to (4.12) is taken into account. Note that this also validates, using this operator as an example, the conjecture stated in [22], for which only heuristic argument and numerical support was obtained in the original paper. Convergence of the series is illustrated in Figure  4.3, in which we compare the expansion for the trace of the stress-energy tensor to the NLIE data obtained by a recursive solution of (4.1).    The analytic continuation of the counting function which underlies formula (4.11) to imaginary values of its argument has to take into account that the first determination strip of the G 1 function has a width of π min(1, ξ). This implies that the so-called second determination must be used in the attractive regime; moreover the series must also be partially summed to obtain the expected exponential decay with exponent determined from the breather mass. To avoid these complications, we resort to numerical calculations, and in figure 4.4 it is shown that that the exact expectation value of Θ is well reproduced by the conjectured expansion. We can also provide the finite temperature corrections to the vacuum expectation value of all the operator e ikβφ R for integer k. Their soliton-antisoliton form factors are finite in the diagonal limit and can be straightforwardly derived from [37] for k = 1, 2 . . . , while our computation provides the two-particle connected diagonal matrix elements in (B.20,B.23). However, at present we can not make a useful comparison to an independent calculation. The NLIE only provides access to the cases k = ±1, i.e. the trace of the stress energy tensor. For higher k, one could resort to numerical determination using truncated conformal space approach (TCSA), originally invented by Yurov and Zamolodchikov [44] and extended to perturbations of the free boson theory in [45]. However, the evaluation of matrix elements is plagued by ultraviolet differences. Using the results in [46], it can be easily seen that the matrix elements of e ikβφ are only finite in TCSA for ξ < 1/(2|k| − 1). Even for k = ±2 this is deep in the attractive regime, with numerous light breather states. For these couplings, before getting to the interesting novel part of our expansion which involves non-diagonal scattering (i.e. the two-soliton states), a lot of corrections coming from states with diagonal scattering must be summed over. Besides this being a very tedious task, the part we would really wish to put to the test would be so tiny as to escape any useful comparison.
An alternative possibility is to renormalize TCSA to obtain the expectation values at a less attractive, or even repulsive coupling. Such renormalization has already been performed for energy levels [47,48]. By extending the methods of [46] it should be possible to extend the procedure to expectation values. However, a preliminary investigation that we performed indicates that this is a very nontrivial task, and is clearly out of the scope of the present work. We hope to return to this issue in the future.

Conclusions
In this work we have studied the low-temperature expansion for one-point functions in sine-Gordon model, which can be considered as a paradigmatic example of integrable field theory with non-diagonal scattering. Following the ideas of Pozsgay and Takács [6], we proposed a series expansion which enables to compute the vacuum expectation value of exponential operators. The formula is expressed in terms of connected components of the diagonal matrix elements of the operator and can be considered as a generalization of the LeClair-Mussardo expansion [1] to non-diagonal scattering. However, in contrast with the latter and with the established formalism in diagonal theories, we could not write our result in terms of single particle energies and momenta dressed by the thermodynamic Bethe ansatz. In this respect, there remains a marked difference between diagonal and non-diagonal scattering theories. Our results were verified for the case of the trace of the stress-energy tensor by comparing against the NLIE approach.
We have also developed a way to evaluate the connected diagonal matrix elements analytically, starting from the integral expressions of the form factors obtained by Lukyanov [37]. For a special value of the exponent, the form factors were compared to those of the trace of the stress-energy tensor which can be obtained by different routes, providing a nontrivial validation of the procedure.
Finally, we have provided analytic support to the relation between finite and infinite-volume form factors conjectured in [22] which had previously been supported only by intuitive reasoning and numerical evidence.
An interesting open problem is to extend our results to states containing more than two solitons/antisolitons, which would then enable the explicit evaluation of the higher terms in the proposed series expansion.

A The diagonal matrix elements of the trace of the stress-energy tensor
Here we summarize the calculation, introduced in [1,49], of the diagonal matrix elements of the operator Θ = T µ µ . First we recall that for a given local operator O, the form factor dependence under space and time translations can be written as in which the energy and momentum of the j-th particle, having mass m j , have been parametrized as e j = m j cosh θ j and p j = m j sinh θ j , respectively. Conservation of the stress tensor implies that: for some scalar field A. Knowing the form factors of this field, as well as the property (A.1), allows to compute those of Θ by the use of: . . , θ n + iπ + ǫ n , θ n , . . . , θ 1 ) where the overall normalization N is left undetermined.
Following the procedure explained in [49], one can determine the two-particle form factor from the expectation value of the Hamiltonian (the T 00 component of the stress energy tensor) by comparing it with the single particle energy The above formula (A.1) implies that the behavior of the two particle form factor is in the diagonal limit ǫ → 0. In order to match Lukyanov's normalization [37] of the exponential operator, an overall normalization has to be left undetermined. This normalization can be determined by comparison with the exact formula (4.12). Analogous reasoning and the expression (2.2) allows one to compute the breather diagonal matrix elements.
Once the proper normalization factor is fixed, higher form factors are uniquely determined and can be computed recursively the by repeated use of the kinematical pole equation, which encodes the singular part of the function when θ m → θ n + iπ: where C is the charge conjugation matrix, which in the case of sine-Gordon is the Pauli matrix σ x in the soliton-antisoliton sector, while it is the identity in the breather sector. The mutual locality factor ω OΨ encodes the braiding properties of the operator O with the field Ψ that interpolates particles [40]. One needs to select all the contributions which diverge as O 1 ǫjǫ k , which will give a finite contribution when inserted into (A.3).

B The connected diagonal soliton form factors of the exponential field
We focus here on the four-particle diagonal matrix element. We need to collect the parts which stay finite when ǫ 1,2 → 0. According to the considerations in section 3.1, all the form factors in the soliton sector share the same structure: there is an overall factor, which depends only on the rapidities but not on the particle species, and a double integral. The double integral is generally time-consuming to evaluate numerically. It has the following form: two functions g and h, both of which depend only on one of the two integration variables, multiplyingḠ, which instead depends on the difference of the two. Therefore, the double integral can be written as a sum of products of two simple integrals by a simple trick [41]: dγ 1 e (α1+α2/ξ)γ1 g (γ 1 ) using the definition (3.7). Note that the integrals are evaluated along the contours Γ (1,2) , which will be deformed to either the real axis or the lines ℑmγ 1 = π (as for the one in figure 3.2). In the course of this deformation some poles are encountered; in the regime 1/2 < ξ < 2, one only needs to treat principal poles of the W functions (3.5). For each integral, there is a contribution of order 1 ǫ1 , which we write as 1 ǫ1 P (j) 1 (θ, ǫ 1 , ǫ 2 ), and another one of order 1 ǫ2 , written as 1 ǫ2 P (j) 2 (θ, ǫ 1 , ǫ 2 ). The index j = 1, 2 represents the integral from which the given contribution originates. We postpone for the moment the analysis of the functions P , by just remarking that they depend on the difference of rapidities θ = θ 2 − θ 1 only.
From this, one has one family of form factors, in the form: Each of the two integrals can now be written as a real integral and is finite whenever ǫ 1,2 → 0. On the other hand, they can have O(ǫ 1,2 ) contributions. We separate each term in the sum and the various orders in ǫ 1,2 by writing where α 1,2 = ± label the terms in the sum (B.1), while σ = ± labels the contour index in the sum (3.2). The prefactor A can be subjected to the same analysis. Its finite part and O(ǫ 1,2 ) contributions can be isolated as The same reasoning can be applied to the form factors in which particles of opposite charge are adjacent. In this case we have: 1;σ2α1α2 (θ, ǫ 1 , ǫ 2 ) ǫ 1 + P 2;σ2α1α2 (θ, ǫ 1 , ǫ 2 ) where we introduced the parametrization and A sasa (θ 1 + iπ + ǫ 1 , θ 2 + iπ, θ 2 + ǫ 2 , θ 1 ) = A sasa + ǫ 1 A sasa,1 + ǫ 2 A sasa,2 + ǫ 1 ǫ 2 A sasa,12 (B.6) Note that in this expression, due to the presence of an integration in the (3.2), it is more convenient to expand around the antisoliton rapidity. The only singularities are coming from the residues picked up from the principal poles of the W functions, during the process of deforming the contours, while the remaining integrals are regular. Note that because of the presence of two antisolitons, each of the contours generates a singularity. However, the contraction of the vertexes e −iφ(ϑ1) e −iφ(ϑ2) is zero for coinciding rapidities, hence the poles at ǫ 1 → 0 and ǫ 2 → 0 are simple.

B.1 The integral parts
To regularize the integral representation, we borrow a method from the original paper [40] and note that the integrals associated to antisolitons can be interpreted as the analytic continuation of Fourier transforms to imaginary arguments. Since the diagonal part of the form factor written in terms of integrals over hyperbolic functions only, one is able to compute the Fourier transformŝ explicitly, and then the resulting convolution can be evaluated numerically. Let us define the functions: which we can use to write Their Fourier transforms areĈ with the obvious symmetrieŝ On the other hand, the Fourier transform of the logarithmic derivative of the W function is easily obtained asL sinh π(ξ−1)z 2 sinh πz sinh ξπz 2 e izθ/2 (B.13) and also, using the definitions (3.6) (B.14) Using this procedure, one obtains the integral part of the diagonal F sa form factors in the following form: For the F ss , we have instead:

B.3 Collecting the finite contributions
Finally, we write the factor A sasa from the contraction (3.23) as in (B.6): A sasa =Ḡ(θ 21 − iπ) Substituting into (B.5) and selecting the order O(ǫ 0 1,2 ) results in 2;σ1α1α2 J (2,2) σ1α1α2 + J (1,1) σ1α1α2 P 1;σ2α1α2 + J (1,2) σ1α1α2 P and the functions P (j) k from section B.2. Using this result we have checked explicitly that F as can be obtained by complex conjugation, as it is, possible to put each term in (B.20) in correspondence with the terms in F as .