LeClair-Mussardo series for two-point functions in Integrable QFT

We develop a well-defined spectral representation for two-point functions in relativistic Integrable QFT in finite density situations, valid for space-like separations. The resulting integral series is based on the infinite volume, zero density form factors of the theory, and certain statistical functions related to the distribution of Bethe roots in the finite density background. Our final formulas are checked by comparing them to previous partial results obtained in a low-temperature expansion. It is also show that in the limit of large separations the new integral series factorizes into the product of two LeClair-Mussardo series for one-point functions, thereby satisfying the clustering requirement for the two-point function.


Introduction
One dimensional integrable models are special interacting many body systems, where the eigenstates and eigenenergies can be computed with exact methods. Various forms of the method called Bethe Ansatz have been developed that apply to a wide variety of models including spin chains, continuum models, and Quantum Field Theories [1,2,3]. Interest in integrable models has been sparked by recent experimental advances [4]: it has become possible to measure physical quantities of such systems both at equilibrium and in far-fromequilibrium situations. In order to compare to experimental data it is essential to calculate the correlation functions.
However, computing exact correlations in Bethe Ansatz solvable models is a notoriously difficult problem. Depending on the models a number of methods have been developed; they include the Algebraic Bethe Ansatz [1], methods based purely on symmetry arguments (for example the vertex operator approach to spin chains [5]), or the so-called form factor approach. Here we do not attempt to review the vast literature, instead we focus on the form factor method.
The main idea of the form factor approach is to evaluate two-point functions as a spectral sum over intermediate states. By definition, the form factors are the matrix elements of local operators on infinite volume scattering states, and they are naturally related to finite volume on-shell matrix elements [6,7]. Traditionally there are two ways to obtain the form factors: either by solving a set of functional relations that follow from factorized scattering and relativistic invariance [8,9,10,11,1,7], or by explicitly embedding the local operators into the Algebraic Bethe Ansatz framework [12,13,14,15]. In many cases these methods lead to explicit and compact representations of the form factors. On the other hand, the summation of the spectral series is typically a very challenging problem, and its treatment depends on the specifics of the physical situation.
First of all, the form factor approach was applied in massive Integrable QFT in order to obtain two-point functions in the physical vacuum [10]. In these cases only the so-called elementary form factors (matrix elements between the vacuum and a multi-particle state) are needed, and the resulting integral series has good convergence properties. Typically it can not be summed up analytically, but a numerical treatment gives highly accurate results [16].
A completely different situation arises when there is a finite density of excitations in the system. Examples include the ground states of certain non-relativistic models such as the anti-ferromagnetic spin chains or the Lieb-Liniger model, or finite temperature situations. Quantum quenches also belong to this class of problems: the long-time limit of local correlation functions can be evaluated on a finite density background given by the so-called Generalized Gibbs Ensemble [17,18,19]. In these cases the form factors and correlation functions display different types of singular behaviour, depending on how the finite density state is constructed. Starting with the infinite volume form factors with a fixed number of particles one encounters the so-called kinematical poles, whose treatment requires special care. On the other hand, starting with a finite volume and increasing the number of particles proportionally with the volume can lead to a non-trivial scaling behaviour for the transition matrix elements [20]. In the XXZ spin chain and related models this approach was used to compute the long distance behaviour of correlations (see [21] and references therein), by studying the asymptotics of explicit determinant representations of the form factors, and performing the relevant summations.
In the context of integrable QFT (iQFT) a framework was proposed in [22] to deal with the kinematical poles of the form factors in finite temperature situations. Integral series were derived for one-point and two-point functions. The resulting series (today known as the LeClair-Mussardo series) are built on the basic form factors with a finite number of particles. Additionally, they involve a thermodynamic function that describes the distribution of Bethe roots, and in the case of the two-point function a thermodynamic dressing of the energy and momenta of the intermediate particles. Arguments and counter-examples against the LM series for two-point functions were presented in [23,24,25,26,27], whereas the result for the one-point function was still believed to be true [23].
In [28] it was shown that the LM series for one-point functions follows from a finite volume expansion of mean values, which uses the so-called connected limit of the infinite volume form factors. The expansion itself was conjectured in full generality in [29] (see also [23]), and for the non-relativistic Lieb-Liniger model it was proven using Algebraic Bethe Ansatz in [28]. Finally, the expansion was proven also for iQFT in the recent work [30], which thus completed the proof of the LM series for one-point functions.
On the other hand, the problem of the finite temperature or finite density two-point functions has remained unresolved. Multiple works performed a low-temperature expansion and obtained explicit formulas for the first few terms [31,32,33,34,35,36]. These results were free of any singularities, and they could be understood as the first few terms in an expansion of a hypothetical LM-type series, but it was not clear what the general structure of this integral series should be.
As an alternative approach it was suggested in [25,26,27] that finite temperature correlations should be computed using form factors that take into account the dressing due to the finite density background. This is in contrast to the logic of the LM series, which uses the zero-density form factors calculated over the vacuum. Even though there is clear physical motivation for this approach, the proposed program has only been applied to free theories. Quite interestingly, a similar picture emerged in the recent work [37], which considered correlation functions in generic inhomogeneous non-equilibrium situations within the framework of Generalized Hydrodynamics. An integral series was derived for the large time and long distance limit of correlations, which would apply as a special case also to static correlations in an arbitrary finite density background. The results of [37] only concern the large time limit, nevertheless it is remarkable that the integral series takes the same form as suggested by the works [25,26,27].
It is somewhat overlooked in the literature, that the formalism of the LeClair-Mussardo series was already developed much earlier in the context of the infinite volume, non-relativistic Quantum Inverse Scattering Method [38,12,39,40,13,41]. In the particular case of the 1D Bose gas explicit formulas were obtained for the form factors of the field operators [12] and the particle current [39] from the so-called Quantum Gelfand-Levitan method. In [40] these were shown to satisfy a set of functional relations that are known today as "form factor axioms" in iQFT. Furthermore, an integral series for the two-point function was also derived in [13,40], that has the same structure as the LM series for one-point functions. In these works the two-point function is treated as a composite object, and the resulting series is built on the form factors of the bi-local operator. In this approach there is no need for an insertion of intermediate states, because the underlying method (the Quantum Gelfand-Levitan equation) allows for an explicit representation of the bi-local product of field operators in terms of the Faddeev-Zamolodchikov creation/annihilation operators.
This program remained confined to the 1D Bose gas, essentially due to the fact that in iQFT the form factors are determined from the solution of the form factor axioms, their structure is considerably more complicated (for example the operators don't preserve particle number), and there is no efficient method to treat the matrix elements of bi-local operators. Moreover, even in the 1D Bose gas alternative approaches (such as the finite volume Algebraic Bethe Ansatz) became dominant, because they lead to intermediate formulas that are more convenient for subsequent analytic or numerical analysis.
In a completely independent line of research the paper [42] developed an alternative framework to deal with finite density states in integrable models. Here the goal was to provide a field theoretical derivation of the TBA equations, by computing the mean value of the Hamiltonian density using its form factor series. The main idea of this work is to consider smeared states such that one does not hit the singularities of the form factors directly. Although this work only considered the Hamiltonian density, its derivations only rely on the general properties of the form factors, therefore all of the intermediate results (before specifying the operator through its form factors) are valid for general one-point functions. This would imply an independent proof of the LeClair-Mussardo series for onepoint functions. We believe that this connection has not yet been noticed in the literature.
An important lesson of the works [12,39,40,13] is that the bi-local operators satisfy the same kinematical pole equation as the local ones. In the 1D Bose gas this was established using explicit form factors for the bi-local operators. Therefore, any type of regularization procedure that treats the kinematical singularities has to work equally well for the one-point and two-point functions. In the present work we build on these ideas: we present arguments for the validity of the kinematical pole equation for the bi-local operators even in integrable QFT, and derive a well-defined LeClair-Mussardo formula for the two-point function.
The article is composed as follows. In Section 2 we review previous approaches towards the LM series, both for the one-point and two-point functions. In Section 2.3 we also formulate our main result using the form factors of the bi-local operators. These are computed in Section 3 by inserting a complete set of states between the two operators. The properties of the bi-local form factors are studied in Section 4. Two different compact representations for the LM series are derived in Section 5, where we explicitly prove the clustering property of the integral series for the two-point function. Our results are compared in Section 6 to earlier calculations in a low-temperature expansion. Finally, Section 7 includes our conclusions.

The LeClair-Mussardo series
In this work we consider one-point and two-point functions in massive, relativistic, integrable QFT. We limit ourselves to theories with one particle species with mass m. The scattering phase shift will be denoted by S(θ), where θ is the rapidity variable. We put forward that the generalization of our results to more particles with diagonal scattering is straightforward, but theories with non-diagonal scattering pose additional technical challenges, which are not considered here.
Our goal is to derive integral representations for the mean values of the one-point and two-point functions in finite density situations. A finite density state can be characterized by a density of rapidities ρ r (θ) such that in a finite volume L the number of particles between θ and θ + ∆θ is ∆N = 2πρ(θ)∆θ. As usually we also define the density of holes ρ h (θ), which satisfies the integral equation is the scattering kernel and p ′ (θ) = m cosh θ is the derivative of the one-particle momentum. We also define the filling fraction .
The physical applications include the standard Gibbs and the Generalized Gibbs (GGE) ensembles [17]. In the first case the thermal average is defined as where β = 1/T . Similarly, for the GGE Tr e − j βj Qj . (2.6) In both cases the ensemble average can be simplified to a single mean value on a representative state |Ω , whose root density ρ r (θ) is determined by the Thermodynamic Bethe Ansatz (TBA) equations [44]. Defining the pseudo-energy as ε(θ) = log ρ h (θ) ρr(θ) , the standard TBA equation reads with e(θ) = m cosh θ is the one-particle energy. Together with (2.3) this determines the Bethe root distributions. In the case of the GGE the eigenvalue functions of the higher charges also enter the source terms. Our main goal is to develop integral series for the objects using their infinite volume, zero-density form factors, for arbitrary root distributions. We allow for two different operators O 1 and O 2 , even though in practice they are typically the same or simply just adjoints of each other. We will restrict ourselves to space-like separations x 2 − t 2 > 0, for reasons to be discussed below.
In evaluating the two-point function there are two main difficulties that need to be solved. First of all, the finite density state |Ω is not a well defined object if one starts from the infinite volume directly. Instead, different regularization schemes need to be applied that increase the number of particles gradually. Second, one has to deal with the singularities of the form factors, including the disconnected pieces and the kinematical poles. Even after the subtraction of the singular parts it is highly non-trivial to find the remaining finite contributions.
In 2.1 we review previous approaches to the one-point functions that lead to the corresponding LeClair-Mussardo series. Later in 2.2 we review previous attempts for the two-point function and in 2.3 we formulate the main results. However, before turning to the LM series we list here the main analytic properties of the form factors of local operators, and discuss the diagonal limit with a finite number of particles.
A generic form factor can be expressed with the so-called elementary form factors by applying the so-called crossing relation (2.9) Here the second line includes the disconnected terms. The elementary form factors satisfy the relations [43,45,46,47,48] I. Lorentz transformation: where s O is the Lorentz-spin of the operator. II. Exchange: F n (θ 1 , . . . , θ k , θ k+1 , . . . , θ n ) = S(θ k − θ k+1 )F n (θ 1 , . . . , θ k+1 , θ k , . . . , θ n ). (2.11) III. Cyclic permutation: F n (θ 1 + 2iπ, θ 2 , . . . , θ n ) = F n (θ 2 , . . . , θ n , θ 1 ). (2.12) IV. Kinematical singularity: There is a further relation related to the bound state structure of the theory, but it will not be used in the present work. On the other hand, we will assume that the form factors show the asymptotic factorization property, when a subset of the rapidities is boosted to infinity: First observed by Smirnov [43] and later proven in [49], this relation holds for relevant scaling operators, and it is used to identify solutions to the form factor axioms with concrete operators. For future use it is useful to display the kinematical pole relation in the form (2.15) The diagonal limit is reached by setting n = m and letting ϑ j → θ j . In this limit the form factor has an apparent n-fold pole, but a straightforward calculation shows that the residue is actually zero. It can be shown that around this point the form factor behaves as F n,n (ϑ n , . . . , ϑ 2 |θ 2 , . . . , where the coefficients A i1i2...in are symmetric in the n indices. There are two natural ways to define a regularized diagonal form factor. The so-called connected form factor is defined as F n,c (θ 1 , . . . , θ n ) ≡ F.P. F 2n (θ 1 + ε 1 , . . . , θ n + ε n |θ n , . . . , θ 1 ) , (2.17) where F.P. stands for finite part, i.e. the terms which are free of any singularities of the form ε j /ε k . According to the expansion (2.16) this coincides with n!A 12...n . The second possibility is to define the symmetric limit as F n,s (θ 1 , . . . , θ n ) ≡ lim ε→0 F 2n (θ 1 + ε, . . . , θ n + ε|θ n , . . . , θ 1 ).
Both diagonal form factors are symmetric in their variables. The linear relations between them can be found using (2.15); they were studied in detail in [29]. Further analytic properties of the diagonal form factors were studied in [50].

One-point functions
In the seminal work [22] the following result was proposed for the finite temperature one-point functions: where ε(θ) is the solution of the TBA equations (2.7) and F n,c are the connected diagonal form factors defined in (2.17).
The main idea of [22] was to consider the finite-temperature problem in finite volume L and in the low-temperature limit T ≪ m. In this case the volume parameter can be chosen to satisfy Le −m/T ≪ 1 such that the partition functions in (2.5) are dominated by states with few particles. In this case it is possible to perform an expansion in the small parameter e −m/T such that the disconnected terms in the numerator are canceled by the Boltzmann-sums in the denominator. A formal calculation gives disconnected terms proportional to δ(0), which were interpreted in [22] as diverging terms proportional to the volume. It was argued that each order in the cluster expansion becomes finite after the cancellation of all Dirac-deltas. However, this procedure only kept to most divergent pieces in L, and sub-leading singularities can also affect the remaining finite answer.
A more rigorous approach was initiated in [29] which aimed to evaluate the average (2.5) in the small temperature limit by keeping all diverging pieces polynomial in L. This was achieved by developing a precise description of the finite volume diagonal matrix elements. In the following we briefly review the results of [29].
Let us denote finite volume states by Here it is understood that the rapidities solve the Bethe equations where p j = p(θ j ) = m sinh(θ j ). For transition matrix elements it was found in [6] 20) where ρ N and ρ M are Gaudin determinants: They can be interpreted as the density of states in rapidity space, and in non-relativistic models as the exact norm of the Bethe Ansatz wave function. The relation (2.20) simply states that (apart from the physically motivated normalization) the form factors are the same in finite and infinite volume. For a coordinate Bethe Ansatz interpretation of this statement see [7]. For finite volume mean values the following result was found in [29]: where J + is the sub-matrix of J corresponding to the particles in the set {θ + }. There is also an alternative representation built on the symmetric diagonal form factors (2.18), but it will not be used here. A rigorous proof of (2.22) was given for local operators of the Lieb-Liniger model in [28], and for relativistic QFT it was finally proven in [30]. It is important that [30] used the kinematical pole relation and the result (2.20) for off-diagonal matrix elements.
In [29] the expansion (2.22) was used to perform a rigorous low-temperature expansion of the Gibbs average, and the LM series was confirmed up to third order. These calculations were performed in the regime Le −m/T ≪ 1, such that it was enough to consider states with low particle numbers.
An alternative, all-orders proof was later given in [28]. Here the idea was to consider a representative state at some temperature (or with some fixed root distribution) and to perform the thermodynamic limit directly on the formula (2.22). In this approach the physical amplitudes F O 2n,c {θ − } of (2.22) enter the integrals of the LM series, and the ratios of the determinants produce the weight functions 1/(1 + e ε(θ) ). More generally it was shown that for an arbitrary state we have where f (θ) is the filling fraction (2.4). It was also shown that the LM series can be expressed using the so-called symmetric diagonal form factors defined in (2.18) as where The two formulas (2.24) and (2.25) can be considered to be partial re-summations of each other. Results similar to (2.25) (including the weight function ω(θ)) had been obtained for the Lieb-Liniger model using Algebraic Bethe Ansatz [51,1]. It is remarkable, that results of the form (2.24) were obtained much earlier in the context of the 1D Bose gas in [38,12,39,40,13]. These works considered both one-point and twopoint functions, for example the field-field correlation Ψ † (x)Ψ(0), which in the x → 0 limit becomes the particle density operator. An integral series with the same structure as the LM series was derived for Ψ † (x)Ψ(0), and in the x → 0 limit the Yang-Yang thermodynamics was obtained completely independently from the TBA arguments [38]. The main idea of these works is to introduce a regularization involving a Galilean boost operator [38], which renders the various singular terms finite. A very important result of this approach, presented in [40], is that the form factors of the two point function Ψ † (x)Ψ(0) satisfy the same kinematical pole equation (2.13) as those of the local operators, see for example (3.11) of [40].
An alternative infinite volume regularization scheme for one-point functions was developed in [42], which considered smeared states to avoid the singularities of the form factors. The main goal of [42] was to provide a field theoretical derivation of the TBA equations, by calculating the mean value of the energy density operator in a finite density state. Although the regularization in [42] is different from that of [38], the treatment of the form factors is the same and the calculation relies only on the kinematical pole property. The work [42] only considered the energy current operator, but its methods could be adapted in a straightforward way to arbitrary local operators, and this would give an independent rigorous proof of the LM series (2.19). We believe that this connection has not been noticed in the literature before.

Earlier results for two-point functions in iQFT
In [22] the following series was proposed for the finite temperature two-point functions: can be interpreted as the dressed momentum and it is given by where ρ 1 (θ) is the solution of the integral equation The form factors appearing in (2.27) are defined by It is an important feature that the x and t dependent phase factors involve the dressed energies and momenta of the particles, where the dressing is due to the finite density background. On the other hand, the form factors are the bare quantities. An explicit counterexample to (2.27) was found in [23]; this counter-example involved a chemical potential, and it was not clear whether (2.27) could still hold for the µ = 0 case for which it was originally derived. In [24] the temperature dependent two-point function of the stress-energy operator T (x) was evaluated in the scaling Lee-Yang model, in the massless limit. The results were compared to benchmark calculations from CFT. The work [24] only considered terms of the LM series up to N = 2, but after an investigation of the convergence properties of the series it was concluded that (2.27) can not be correct. In [25,26,27] it was shown that the LM series is not correct in the Ising model for most operators; it only works for operators with at most two free field factors.
Also, there is a central problem with the proposal (2.27), which is independent from the previous counter-examples: the higher terms with N ≥ 3 involve ill-defined integrals and it is not specified how to subtract the singular pieces. For N ≥ 3 all the terms for which the σ variables are not equal have poles in the θ-variables. These poles are squared and there is no prescription given for the subtraction of these singularities. As far as we know, this problem has not been emphasized in the literature. The calculation of [24] did not encounter this problem, because they only included terms up to N = 2, for which the pre-factor of the double pole is zero.
Later works attempted to perform a well-defined low-temperature expansion of the Gibbs average for the two-point function. In [31,32,33,34] form factor series were developed for spin chains and field theoretical models. The paper [34] includes both finite and infinite volume regularization, involving also a constant shift in the rapidity parameters to treat the singularities of the form factors; this corresponds to adding a Lorentz-boost operator into the definition of the partition function. Adding a boost is essentially the same technique that was used in [38,12,39,40,13], where the Galilean boost operator was used due to the non-relativistic kinematics of the Lieb-Liniger model.
The finite volume regularization was applied later in [35,36], where the main goal was to derive a regularized form factor expansion irrespective of the details of the model or the operator. The summation over Bethe states in the Gibbs average was transformed into contour integrals such that the infinite limit volume could be taken in a straightforward way. However, a shortcoming of this method was that it required a term by term analysis, and it was not evident how to express the general higher order terms.
All of these approaches use the bare form factors, that is, the infinite volume zero density form factors of the theory. In [25,26,27] it was suggested that finite temperature correlations should be described by dressed form factors that already take into account certain effects of the background. New form factor axioms were set up for the finite temperature case, and they were solved for problems involving Majorana fermions. However, the approach has not yet been worked out for interacting theories. On the other hand, an exact result was computed in [37] for the large scale correlations in generic inhomogeneous, non-equilibrium situations, which also include static backgrounds as special cases. Here we do not discuss the results of [37] in detail, as they apply to the large time limit; however, in the Conclusions we give remarks about the possible relations to our work.

LeClair-Mussardo series for two-point functions
In the following we derive a new form factor series for the two-point function in finite density situations. The structure of our result is essentially the same as that of the formula (2.19), therefore it can be called the LM-series for the two-point function. The central idea is to treat the product of two local operators as a composite object, and to investigate its matrix elements.
Let us define the form factors of the bi-local operator It is important that we restrict ourselves to space-like separations x 2 − t 2 > 0, in which case the two operators commute with each other. This restriction is essential for some of our arguments to be presented below. The form factors are always Lorentz-invariant, therefore we could transform any |t| < |x| to t = 0 with a proper boost; however, the finite density state |Ω is typically not Lorentz-invariant, therefore we reserve the possibility of a finite time displacement. In the following we present three independent arguments that show that the bi-local form factors satisfy the same kinematical pole property (2.15) as the usual form factors.
Our first argument is based on the analytic properties of the Operator Product Expansion (OPE) 1 Writing the bi-local operator as the OPE each term on the r.h.s. satisfies the analytic properties, therefore their sum does too. This argument relies on the existence and the absolute convergence of the OPE. Whereas this has not been rigorously proven, it is generically believed to be true in QFT for space-like separations [52,53,54,55,56].
Our second argument relies on the coordinate Bethe Ansatz wave function. It was shown in [7] that in non-relativistic cases the kinematical pole is easily proven by investigating the real space integrals in the form factor. For this argument we require t = 0, but here we deal with only a finite number of rapidity variables (independent of the background |Ω ), therefore we can always set t to zero by an appropriate Lorentz-transformation. We argue that the bi-local form factor satisfies the relation which has the same form as (2.15). The singularity can be understood easily by real space calculation of the form factors [7]. The two terms in the pole represent divergent real space integrals when the particles with rapidities θ 1 and ϑ 1 are before or behind the operator and all the other particles, and the pre-factor reflects the change of the phase of the wave functions as the particles with θ 1 and ϑ 1 are moved from the first position to the last. The kinematical pole arises from infinite x → ±∞ real space integrations, therefore it is not sensitive to the precise locality properties of the operators, the only requirement being the product should have a finite support in real space. Similarly, relativistic effects that modify the Bethe Ansatz wave function on small distances (comparable to the Compton wavelength of the particles) do not play a role either, because the non-relativistic derivation of [7] only uses the long distance behaviour of the wave function, which is described by the Bethe Ansatz even in QFT. Therefore, this argument also supports the kinematical pole property even for the bi-local operators.
The third argument is based on explicit representations of the bi-local form factors in terms of the form factors of the individual operators; such integral formulas will be presented in Section 3. The kinematical pole is then evaluated explicitly in Section 4.3. This derivation is mathematically rigorous, but it leaves the physical meaning of the kinematical residue somewhat obscure. We believe that it is our first two arguments which provide the physical understanding.
In 2.1 we reviewed three different approaches that lead to the LM series for one-point functions: • The infinite volume regularization of [38,12,39,40,13] that uses a boost operator.
• The infinite volume regularization of [42] that uses smeared states.
• The finite volume regularization of [29] (supplied with the proof of [30] for the diagonal matrix elements), and the presentation of [28] regarding the thermodynamic limit.
All three approaches use only the kinematical pole property, and we have argued above that this property holds for the bi-local operators as well. It follows that finite density two-point functions are given by the LeClair-Mussardo series where G x,t n,c is the connected diagonal form factor of the bi-local operator, which is defined as the finite part of G x,t n,n (θ 1 + ε 1 , . . . , θ n + ε n |θ n , . . . , θ 1 ), which is free of any singularities of the form ε j /ε k . Alternatively, the LM series can be expressed in the same form as (2.25): where G x,t n,s (θ 1 , . . . , θ n ) = lim ε→0 G x,t n,n (θ 1 + ε, . . . , θ n + ε|θ n , . . . , θ 1 ), (2.35) and ω(θ) is given by (2.26). The equivalence between (2.33) and (2.34) is guaranteed by the theorems for the connected and symmetric diagonal form factors presented in [29]; these theorems only use the kinematical pole property, therefore they apply to the two-point function as well. Expressions (2.33)-(2.34) are implicit: they do not specify how to compute the form factors of the bi-local operator, they only describe how to deal with the singularities of this object. In the next section we also show how to compute the bi-local form factors using those of the local operators. This is achieved by inserting a complete set of (infinite volume) states between the two operators, which makes the LM series completely explicit. It is very important that the composite form factors depend on the position (x, t), and the diagonal limit has to be taken by considering a full dependence on the rapidities, including the kinematical factors involving (x, t). Examples in Sec. 6 show that this can lead to secular terms.
We remark again, that in the non-relativistic case an LM series of the form (2.33) has been already established in [40], see for example eq. (3.18) in that work, together with (3.14), which agrees with our definition for G x,t n,c . An important difference between the two situations is that in the 1D Bose gas the field operators change particle number by one, and the bilocal form factors could be obtained explicitly as sums of algebraic expressions. On the other hand, in the field theoretical case the operators can have matrix elements between any Nand M -particle states, and the bi-local form factors are obtained as an infinite integral series, to be presented in the next section.

Form factors of the bi-local operators
In this Section we derive explicit integral representations for the generic n-m form factor G x,t n,m (ϑ n , . . . , ϑ 1 |θ 1 , . . . , θ m ) = ϑ n , . . . , by inserting a complete set of states between the two local observables. The idea and the methods of this section are essentially the same as in the original works by Smirnov [43]. The analytic properties of the resulting series are investigated in Section 4. A naive insertion of states would lead to an encounter with the kinematical poles: this happens when intermediate rapidities approach one of the ϑ or θ variables. In order to avoid these singularities we employ a well-defined expansion of the local operators in terms of the Faddeev-Zamolodchikov (FZ) creation/annihilation operators Z † (θ), Z(θ) [57,58].
The FZ operators satisfy the commutation relations Local operators are represented in terms of the FZ operators as [43,59,60] O where where the functions f can be expressed in terms of the form factors f O k,l (θ 1 , . . . , θ k |η l , . . . , η 1 ) = F O l+n (θ k + iπ + i0, . . . , θ 1 + iπ + i0, η 1 − i0, . . . , η l − i0) (3.5) and K x,t is the phase factor due to the displacement of the operator: where E(θ) = m cosh(θ) and P (θ) = m sinh(θ) are the energy end momentum of the oneparticle states. An important ingredient in the expansion above is the presence of the infinitesimal shifts of ±i0. They are irrelevant for the off-diagonal form factors of O, but they are necessary to obtain a well-defined integral representation for the bi-local form factor (3.1).
The bi-local form factor is obtained by inserting two instances of the expansion (3.3)-(3.4) into (3.1) and computing the contractions of the FZ operators using the algebra (3.2). It is easy to see that only those terms contribute, where there are at least n FZ creation and m annihilation operators. Also, there can not be more than m annihilation (or n creation) operators immediately acting on the ket (or bra) states. These constraints are satisfied by the triple sum In H n,m k,l,p the indices k and l show the numbers of the disconnected ϑ and θ rapidities. We introduce the notation for the ordered set of rapidities {θ j1 , θ j2 , . . . , θ jn } = {θ} J< and {ϑ jn , ϑ jn−1 , . . . , ϑ j1 } = {ϑ} J> , where the elements of J are ordered as j i < j i+1 . I n,< denotes the ordered set {1, 2, . . . , n} and I n,> denotes {n, n − 1, . . . , 1}.
For each term in H n,m k,l,p we need to evaluate a contraction of the form {ϑ} Im,> = 0|Z(ϑ m ) . . . Z(ϑ 1 ). (3.10) The contractions lead to various Dirac-deltas and phase factors, which can be treated using straightforward calculations, leading to G x,t n,m (ϑ In,> |θ Im, where the two inner summations run over bipartite partitions of the index sets I n and I m , and we omitted the particle number subscript of the form factors. The phase factors S < , S > and S ↔ above result from the exchanges of the FZ operators such that the set of particles {θ A − } can be connected to the operator O 1 , and the set {ϑ B − } to O 2 . They are defined as follows. For a given partitioning }} the ordered set of rapidities that is obtained by concatenating the two ordered sets {θ A + < } and {θ A − < }. The phase factor S < follows from the rearrangement of rapidities as (3.12) Explicitly it is given by where a + i , i = 1, . . . , |A + | and a − j , j = 1, . . . , |A − | are the indices in the sets A + and A − . Similarly, S > is defined as the phase factor coming from the rearrangement with the explicit expression Finally, the mutual phase factor S ↔ is given by The intuitive interpretation of the formula (3.11) is the following: The incoming and outgoing particles are connected to one of the local operators, whereas they are disconnected from the other one; this choice is given by the partitioning of the index sets. On the other hand, the integration variables µ j represent additional intermediate particles between the two operators. The phase factors in (3.11) arise from the interchange of positions of the particles. We stress that the ±i0 shifts for the rapidities are essential to obtain a regular expression for the bi-local form factor. A graphical representation of (3.11) is given in Fig.  1.
It is useful to derive an alternative representation of the bi-local form factor, where the integration contour over the intermediate particles is well separated from the incoming and outgoing rapidities. To this order we shift the integration contour in (3.11) to R + iα with α being a small real number. The advantage of this step is that afterwards the kinematical poles at ϑ j → θ k can be studied directly, without paying attention to the µ-integrals. On the other hand, if some of the µ-integrals would run between the rapidities ϑ j and θ k , then the contour could be pinched, and such cases would complicate the analysis of the kinematical singularities.
First we assume that α > 0, then the contour shift picks up poles of the form factors of O 1 . After a straightforward, but tedious calculation we obtain the formula G x,t n,m (ϑ In,> |θ Im,< ) = ∞ p=0 G x,t n,m,p (ϑ In,> |θ Im,< ), where G x,t n,m,p (ϑ In,> |θ Im, where we omitted the ±i0 shifts, because they are not relevant anymore. Details of this calculation are presented in Appendix A; here we just note that the changes in the phase factors (including the reordering of the rapidities in the form factor F O2 ) at a given p result from adding pole contributions from terms with p ′ > p. We can also shift the contours to the negative imaginary direction, which results in the alternative formula G x,t n,m,p (ϑ In,> |θ Im,< ) = . (3.18) The differences between the integrands in (3.11) and those of (3.17) and (3.18) lie in the phase factors, including the ordering of the rapidities in the form factors functions. The convergence of the µ-integrals in the above expressions depends on the separation between the operators. The µ-dependence for large |µ| is given solely by the K-factors, given that the asymptotic factorization property (2.14) holds. It follows that the integrations towards µ → ±∞ are convergent in (3.17) if x < 0, and in (3.18) if x > 0. If the two operators O 1,2 are not identical then the two-point function might not have the x → −x symmetry, and only one of the representations can be used, depending on the sign of x. On the other hand, the µ-integrals would never be convergent for time-like separations: they would diverge either for µ → +∞ or µ → −∞.
It is important that the integral series is effectively a large-distance expansion for the bi-local form factor: each µ-integral carries a weight of e −ms with s 2 = x 2 − t 2 . This is analogous to the case of the vacuum two-point function.
We note that essentially the same integral representations as presented in this section were used by Smirnov to prove the local commutativity theorem [43]. This theorem states that if there are two operators such that the form factors satisfy the axioms as given above, then the two operators commute for space-like separations. We do not replicate the proof here. Instead, we just point that it also uses a shift of the integration contours in the complex plain, which is only possible for space-like displacement between the operators.

Form factor properties of the bi-local form factors
In this Section we investigate the analytic properties of the bi-local form factors based on the integral representations (3.17)- (3.18). Most importantly we confirm explicitly that the bi-local form factor satisfies the kinematical pole relation (2.32), which is the basis for the derivation of the LeClair-Mussardo series. In 4.4 we also discuss the implications of the form factor axioms for the bi-local operators.

Lorentz transformation
After a Lorentz transformation every rapidity is shifted by Λ. The S-matrix factors are invariant under the shift, since they only depend on rapidity differences. Each elementary form factor satisfies the Lorentz transformation axiom (2.10). The phase factor transforms as

Exchange property
It is easy to see that the form factors of the bi-local operators satisfy the exchange property for rapidities both sides of the operators

(4.3)
For those terms where the rapidities to be exchanged are in the same sets A ± or B ± the property follows from the property of the elementary form factors (2.11). In those cases, where the rapidities are in different sets, the extra S-matrix factor exactly cancels out or introduces the terms in the ) phase factors to get the appropriate ordering.

Kinematical poles
The kinematical pole property (2.13) is the essential ingredient for the proof of the LeClair-Mussardo formula. Here we prove that it holds for the form factors of the bi-local operators too. We use the representation (3.17) in which there are no singularities associated with the intermediate particles. However, the elementary form factors in the summation (3.17) have direct kinematical poles for ϑ 1 ∼ θ 1 , and these need to be summed up to obtain the total residue of the bi-local form factor.

(4.5)
Adding up the two kinds of residues we arrive to the kinematical pole property of the bi-local form factor This is identical to the crossed version (2.32) of the form factor kinematical pole axiom (2.13). We remark that an analogous calculation can be performed also for the alternative representation (3.18).

Other properties and discussion
We have checked that the crossing and periodicity properties (2.9) and (2.12) also hold for the bi-local form factors. However, these relations are not relevant for the LM series, therefore we refrain from presenting the proof. We also remark that the asymptotic factorization property (2.14) does not hold in the bi-local case, which can be understood simply from a physical point of view, or from the integral representations.
So far we have argued that the bi-local form factors satisfy all axioms that were originally derived for the local ones. Here we discuss the implications of this statement.
First of all, it is an interesting idea to apply the local commutativity theorem to the composite objects. This theorem by Smirnov states that if the form factors of two operators satisfy the set of axioms presented in Section 2, then they commute at space-like separations [43]. On the other hand, our results suggest an extension of this theorem to bi-local objects. Such a theorem would state that for any three operators with the required form factor properties we have given that all three separations are space-like. The trivial way of proving (4.7) is to apply Smirnov's theorem twice to commute the two objects. However, a second proof can be given by using the form factor properties of the bi-local operator, and then repeating Smirnov's calculation for the composite object. We refrain from presenting a rigorous proof of this extended theorem, as it is not relevant to the present work. Nevertheless, we stress that the local commutativity theorem can not be applied backwards: it does not imply that only the local operators can satisfy the form factor axioms. Finally we comment on the possibility of a form factor bootstrap for the bi-local objects. In the case of the local operators the form factor axioms (together with the asymptotic factorization property) include enough constraints to fix the form factors completely, both in relativistic and non-relativistic situations. The reason for this is the following: The form factors can always be written as a product of a "minimal" part (which is two-particle irreducible) and a physical amplitude, which is a symmetric polynomial in the appropriate variables. Then the kinematical pole axiom is used to fix this polynomial. In the bi-local case the polynomial would depend also on the displacement between the two operators, and the same number of constraints can not fix it. Therefore, the standard bootstrap procedure can not be applied to the bi-local form factors.

Compact representations of the LeClair-Mussardo series
Formulas (3.17)- (3.18) give an integral representation for G x,t n,m in terms of the elementary form factors. The LM series is then given by (2.33) and (2.34), and this gives a well-defined way to evaluate the two-point function.
The final integral formulas (2.33) and (2.34) require the evaluation of the connected and symmetric diagonal form factors. It follows from the proof of the kinematical pole property, that these operations can be performed on the objects G x,t n,n,p for each n and p separately. On the other hand, computing the diagonal limit term by term in (3.17)-(3.18) leads to singularities, and it is only the sum over the partitions in (3.17)-(3.18) that has the desired residue structure G x,t n,n (ϑ In,> |θ In, One way towards singularity free expressions is to develop an integral formula which automatically produces the diagonal limit of the form factors. This can be achieved by introducing auxiliary integration variables. It follows from (5.1) that where C j are small contours around θ j . These contours can be opened to encircle the whole real line. Opening the contours we do not encounter additional n-fold poles, therefore we can write where C is a narrow contour around the real line and the factor of 1/n! has been inserted due to the possible permutations of the set {ϑ}. It is important that C is narrow enough so that it does not hit the integration contours R + iα or R − iα for the intermediate particles in (3.17)-(3.18). It follows that the LeClair-Mussardo series can be written finally as .

(5.4)
An alternative compact formula can be given by using the symmetric diagonal limit in (2.34). The idea is to exchange the ε → 0 limit with the summations. This is justified (at least in the thermal situation) due to the exponentially decaying factors e −mR cosh(θ) and e −mx cosh(θ) associated with the numbers n and p, and the limits on the growth of the form factors following from the asymptotic factorization property (2.14). We then obtain  G x,t n,n (θ 1 + ε, . . . , θ n + ε|θ n , . . . , θ 1 )   .

(5.6)
We stress that it is important to keep the iε shift in the kinematical factors K x,t , because they can combine with the poles the form factors to produce terms proportional to x and t. Examples for this are given in Sec. 6. On the other hand, the shifts could be omitted from the phase factors, because these factors only depend on the rapidity differences within a given set. It is also important that all ε variables have to be identical, because any other choice with some ε j /ε k = 1 would lead to different finite terms. Representation (5.6) is perhaps the most transparent from a physical point of view: the rapidities {θ} stand for particles that are present in the representative state |Ω and interact with the two-point function, whereas the {µ} are additional intermediate particles between the two operators. A graphical interpretation is given in Figure 2. We also point out the striking similarity with the original proposal (2.27), however, there are crucial differences. In (2.27) the energy and momenta of the intermediate particles is dressed, and it is not specified how to deal with the kinematical poles in the higher order terms. In contrast, here the energy and momenta are not dressed, and all terms are regular due to the well-defined shifts iα and ε. The only "dressed" quantities in (5.6) are the weight functions f (θ) and ω(θ), and the derivation in [28] shows that these are statistical weights that are determined by the Bethe Ansatz description of the state |Ω .
To conclude this section we discuss the clustering property of the two-point function. In the limit of large separations it is expected that This identity is motivated by physical requirements about the pure state |Ω , but its explicit confirmation is a highly non-trivial test of the various integral formulas [35,36]. Here we prove that the LM series satisfies the clustering property; the simplest way is to use the representation (5.6).
The kinematical factors K x,t include multipliers of the form e ±ip(θ)x , which result in terms of O(e −m|x| ) after integration in θ. Therefore, in the large distance limit only those terms survive where all K-factors are trivial. This happens when there are no intermediate particles, and the partitions have to satisfy A + = B − and A − = B + . In these cases the incoming and outgoing θ j rapidities are always attached to the same operator, and we obtain .
The phase factors cancel each other by definition, and then the integrals completely factorize according to the partitions. By permutation symmetry we obtain A graphical interpretation of this identity is shown in Fig. 3.

Comparison to our earlier work
Here we compare the formula (5.4) to our previous results for the finite temperature case [35,36]. These articles approached the evaluation of finite temperature two-point functions through finite volume regularization, leading to a double series   .18), the terms can not be matched one to one. This is an effect of the differences in the regularization schemes, and it can be seen explicitly in the example presented below.
We compare the LM series to (6.1) up to first order in e −mR cosh(θ) . Due to the infinite number of possible intermediate particles between the two operators this is already a strong independent confirmation of (2.33).
In the first order approximation the (2.7) gives simply f (θ) = e −mR cosh(θ) and we expect to match the zeroth and first order contributions as These two terms are evaluated in the following two subsections.
We have simply which is nothing more that the spectral series for the zero-temperature two-point function.
In this case there is no singularity for the µ-variables, and the integration contour could be pulled back to the real line. The same result was given for M D 0M in [35].
We calculate the connected part of G 1 , which is defined as the ε independent part of the bi-local form factor G 1 (θ + iε|θ). Two-particle form factors don't have kinematical poles, therefore in this case the connected part coincides simply with the ε → 0 limit. From (3.17) we have the integral representation (6.5) The first and fourth terms are finite in the ε → 0 limit, but the second and third terms have apparent first order poles at ε = 0. These singularities cancel each other, and care needs to be taken to obtain the finite left-over pieces.
We introduce an expansion of the form factors to make the singularities apparent With the help of the expansion we can express the ε independent part of the integrands in the second and third lines of (6.9). Following [35,36] these terms will also be called the "connected parts". They are given by In the (6.8) the tilde denotes, that the connected form factor incorporates the contribution form the energy and momentum phase factor as well.
With the regular pieces introduced above the connected bi-local form factor is written as The D 1M contributions to (6.1) were calculated in [35,36] as

Conclusions
In this work we derived a LeClair-Mussardo formula for the two-point function at spacelike separations. The general structure of the series is given by the two implicit representations (2.33) and (2.34). These formulas involve the form factors of the bi-local operators, for which the integral series are given by (3.17) or (3.18), depending on the sign of x. More explicit representations are given by (5.4), and finally (5.6). We believe that the latter form is the physically most transparent; its graphical interpretation is plotted on Fig. 2. The result can be considered as a large distance expansion, because each intermediate particle (each µ-integral) carries a weight of ∼ e −ms with s 2 = x 2 − t 2 .
It is an important property of the final formulas that each term is explicit, with a well defined prescription to deal with the kinematical singularities. This was missing in the previous works: the original proposal of [22] had singularities in the higher order terms, and the regularization of [35,36] was only performed on a case by case basis.
Our results are expressed in terms of the bare form factors of the theory, and the momenta and energy of the intermediate particles between the operators are also the bare quantities. This follows from our approach to treat the bi-local operator as a composite object, and to develop the LM series based on the bi-local form factors. The only dressed quantities are the statistical filling fractions f (θ) and the weight function ω(θ). These two are derived from the Bethe Ansatz description of the background |Ω , and they are the same for all local and bi-local operators.
The original proposal of [22] used the bare form factors of the theory, but it suggested a dressing for the intermediate particles. This might seem intuitive, but our derivation shows that it is inconsistent. We expressed the two point function with the bare quantities, and if there is a partial re-summation of the series, then it should lead to a dressing for both the form factors and the kinematical factors. Examples for such dressing can be found in the context of non-relativistic models [1], and a similar structure was also implied by the results of [37], at least in the large time limit.
The methods and results of this paper also apply to non-relativistic models, such as the 1D Bose gas. In fact, most of these results were already derived for the Bose gas in [38,12,39,40,13], and these papers served as a motivation for the present work. The main addition of our work is the realization that the kinematical pole structure of the bilocal operator is the same even in the relativistic case, and that this property is enough to establish the LM series. We presented three arguments for the kinematical pole; the most rigorous is the explicit check presented in 4.3.
The form factors of the bi-local operator were obtained by inserting a complete set of states between the two operators. This is the same technique that was used by Smirnov [43], which lead to the proof of the local commutativity theorem. There is one essential step in both calculations: a certain shift of integration contours, which is possible only for space-like separations.
As tests of our results for the LM series we evaluated the first order corrections in the low-temperature expansion in Section 6 and compared them to existing results; complete agreement was found for an arbitrary number of intermediate particles. As an additional test we performed the large distance limit, and confirmed that the integral series factorizes into the product of two LeClair-Mussardo series, thus fulfilling the clustering property. This is a highly nontrivial test of the final formula (5.6).
The most important open problem is to find methods for the actual implementation of the final formulas. As shown by the examples in Sec. 6, this can be a cumbersome task, both numerically and analytically.
It is also important to consider the case of time-like separations. Our first two arguments for the validity kinematical pole property (those based on the OPE and the Bethe wave function) do not apply in this situation, but preliminary calculations show that the integral series for the bi-local form factor can be modified such that the kinematical pole property can be proven nevertheless. However, the convergence properties of the resulting integral series can be drastically different, and we don't have a decisive answer whether a LM series can be established for this case too.
Coming back to space-like separations, it would be interesting to consider the first corrections in a large distance expansion. Generally it is expected that where ξ is a correlation length. In the thermal situation Euclidean invariance implies, that the above object is equal to a two-point function with time-like separation, evaluated in the ground state of a finite volume QFT with volume R = 1/T . In this case the correlation length is 1/ξ = E 1 − E 0 , where E 0 and E 1 are the exact ground state and first excited state energies in the given volume. It is an interesting question whether these quantities (or possibly even the exact finite volume transition matrix elements that determine the pre-factors of the correction terms in (7.1)) could be extracted from the LM series.
Finally, it would be interesting to establish a connection to the results of [37]. Equation (3.38) in that work is reminiscent of our formulas, even though it was derived for time-like separations in the large time limit. The concrete relations between the two integral series have yet to be investigated.
We hope to return to these questions in future research.
Note added: After this work was finished we became aware of the work [61] which derives the so-called first Lüscher's corrections to finite volume form factors in iQFT. The methods and results there are related to our work. In [61] the new results were shown to be consistent with the previous works [35,36], thus they are consistent with our present results too.
where C i is contour encircling ϑb i with positive orientation. We omitted the ±i0 shifts in the formula, since there are no more poles on the integration contours, the form factor is off-diagonal (ϑ i = θ j ) and the shifts are not necessary anymore. We note, that in the case |A − | = 0, p = |B + | and p = r, the residue is missing, since the two-particle form factor is regular. With (2.11) and (2.13) the residue evaluates to ∞ p=r 1 (p − r)! R+iα .
The next step is to expand the product in the formula. We separate the setB into two disjoint subsetsB I ∪B S =B. For indices inB I we pick the 1-term in the product, for the indices inB S we pick the product of the S-matrices. Using the exchange axiom (2.11) we simplify the expression to where for the set B − ∪B S we imply the rearrangement of the elements to increasing order, and we relabeled the integration variables.
It is important to note that apart from the sign of the expression, it only depends on the setsB,B I and B − ∪B S and p. Since we have to sum for all disjoint B + and B − that sums up to B, there is a chance for cancellation of terms. For fixedB,B I and B − ∪B S and p the sum of the sign pre-factors is