Exact finite volume expectation values of $\bar{\Psi} \Psi$ in the Massive Thirring model from light-cone lattice correlators

In this paper, using the light-cone lattice regularization, we compute the finite volume expectation values of the composite operator $\bar{\Psi}\Psi$ between pure fermion states in the Massive Thirring Model. In the light-cone regularized picture, this expectation value is related to 2-point functions of lattice spin operators being located at neighboring sites of the lattice. The operator $\bar{\Psi}\Psi$ is proportional to the trace of the stress-energy tensor. This is why the continuum finite volume expectation values can be computed also from the set of non-linear integral equations (NLIE) governing the finite volume spectrum of the theory. Our results for the expectation values coming from the computation of lattice correlators agree with those of the NLIE computations. Previous conjectures for the LeClair-Mussardo-type series representation of the expectation values are also checked.


Introduction
Finite volume matrix elements of local operators play an important role in several applications of integrable quantum field theories. Namely, they are fundamental building blocks of the form factor perturbation theory [1] and their determination is indispensable for the computation of the string field theory vertex [2] and of the heavy-heavy-light 3-point functions [3] in the planar AdS 5 /CF T 4 correspondence.
In the past decade a remarkable progress has been made in the computation of finite volume form factors in integrable quantum field theories [4][5][6][7][8][9][10][11]. Most of the methods use the infinite volume form factors [12] as a starting point and the finite volume form factors are to be determined in the form of a systematic large volume series. As a first step the large volume corrections, that decay with a power of the volume were determined [4,5] and not much later a method was proposed [6] for computing some special type exponential in volume corrections. These investigations shed light on the fact that the computation of diagonal matrix elements is a much simpler task than that of the non-diagonal ones. Recently, in [13] it has been shown, that the leading term in the large volume series representation of the diagonal form-factors in [5,15] can be derived from the formulas for the non-diagonal form-factors of [4], by taking the diagonal limit appropriately.
Though the structure of exponentially small in volume corrections for the non-diagonal matrix elements is still unknown, inspired by [14] for the diagonal matrix elements, a nice series representation was proposed in [10,11]. However the proposal is valid only to purely elastic scattering theories and its extension to non-diagonally scattering theories is still unknown in general.
Recently, in the Massive Thirring (sine-Gordon) model a similar series representation was proposed to describe the finite volume diagonal form factors of the theory [17]. The conjecture was based on the computation of the U (1) current from the light-cone lattice regularization [20] of the theory.
The purpose of this paper is two-fold. On the one hand we would like to demonstrate that the light-cone lattice approach admits an appropriate framework for computing the finite volume form factors of the Massive Thirring (sine-Gordon) model and on the other hand we would like to give further justification for the validity of the LeClair-Mussardo type series representation conjectured in [17].
To do so we compute the diagonal form factors of the composite operatorΨΨ from the light-cone lattice approach. There are several advantages of the choice of this operator. First of all, this operator is proportional to the trace of the stress-energy tensor. Thus the results of reference [18] imply, that up to a constant factor these expectation values can be computed simply from the non-linear integral equations (NLIE) governing the finite volume spectrum of the model [22]- [30]. This makes possible to check the results coming from the lattice computations against a result coming from a completely different method. Second of all, the operatorΨΨ is still simple enough not to mix with other operators under renormalization. Nevertheless, contrary to the case of the U (1)-current [17], in this case an infinite renormalization constant arises in accordance with field theoretical computations [19].
In the present paper, using the framework of Quantum Inverse Scattering Method [33]- [57], we compute such special spin-spin 2-point functions on the lattice, in which the spin operators are located at neighboring sites of the lattice. A straightforward computation shows, that the discretized version of the continuum operatorΨΨ corresponds to the lattice operator: σ + n σ − n+1 + σ − n σ + n+1 . We compute the expectation values of these operators between those Bethe eigenstates which correspond to the pure fermion (soliton) states in the continuum theory. Taking the continuum limit, we show that this quantity, as expected, is proportional to the result coming from the pure NLIE computation of the expectation values of the trace of the stress-energy tensor. As a by product, our method accounts for the lattice artifacts, as well. Our results show, that the leading order divergences are exactly of the same form as those expected from the renormalization group analysis of the Massive Thirring (sine-Gordon) model. Finally, we also checked that the all order conjecture [17] for the systematic large volume series representation of the diagonal fermionic (solitonic) form-factors of the Massive Thirring (sine-Gordon) model is also valid for this operator.
The outline of the paper is as follows: In section 2. we recall the most important properties of the Massive Thirring and sine-Gordon models and their light-cone lattice regularizations. This section contains the pure NLIE computation of the fermionic (solitonic) expectation values of the trace of the stress-energy tensor. In section 3. we summarize the Quantum Inverse Scattering Method framework and the lattice part of the computation of the special spin-spin 2-point functions of interest. The continuum limit procedure is described in section 4. In section 5. we rephrase our results in the form of a systematic large volume series and check the validity of the conjecture of [17]. Our summary and outlook closes the body of the paper in section 6. The paper contains three appendices, as well. In appendix A we rewrite the sums entering the lattice formulas for the two-point functions into integral expressions. In appendix B we describe how to compute the lattice cutoff tend to zero limit within these integral expressions. Finally, appendix C contains the large argument series representations of the convolution integrals necessary for the computations.
2 Light-cone lattice approach to the Massive-Thirring and sine-Gordon models The Massive Thirring (MT) model is defined by the Lagrangian: where m 0 and g denotes the bare mass and the coupling constant of the theory, respectively. As usual, γ µ s stand for the γ-matrices. They satisfy the algebraic relations: {γ µ , γ ν } = 2η µν with η µν = diag(1, −1). Throughout the paper we use the chiral representation for the fermions as follows: It is well known [31], that this fermion model can be mapped to the sine-Gordon (SG) model: provided the coupling constants of the two theories are related by the formula: A more detailed investigation of this equivalence [32] pointed out, that the two models are identical only in the even topological charge sector of their Hilbert-spaces and they differ in the odd topological charge sector.
The operator we study in this paper is the fermion bilinearΨΨ in the MT model. According to the equivalence [31] it is proportional to the potential of the sine-Gordon model [19]:Ψ Ψ ↔ 1 πa cos(βΦ), (2.5) with a being a cutoff in coordinate space. The perturbing operator cos(βΦ) of the SG model is related to the trace of the stress-energy tensor as follows 1 : From (2.5) and (2.6) the fermion bilinear can be expressed in terms of the trace of the stressenergy tensor as follows:Ψ Due to renormalization effects α 0 scales with the coordinate space cutoff a as α 0 ∼ a −β 2 /4π [19], thusΨ The minimal length a can be thought of as a lattice constant, as well. From (2.8) it can be seen that the matrix elements ofΨΨ are divergent in the attractive regime (β 2 < 4π) and the coefficient of the leading order divergence is proportional to the trace of the stress-energy tensor. In this paper we show that our light-cone lattice computations account for the scaling behavior (2.8) and up to a constant factor, allow one to compute the expectation values of T µ µ .

The light-cone lattice regularization
The light-cone lattice regularization scheme [20] admits an appropriate lattice approach to the even topological charge sector of the MT model. In this description the space-time is discretized along the light-cone directions: x ± = x ± t with an even number of lattice sites in the spatial direction. The sites of the light-cone lattice correspond to the discretized points of space-time. The left-and right-mover fermion fields live on the left-and right-oriented edges of the lattice. In this manner a left-and a right-mover fermion field can be associated to each site of the lattice (See figure 1.). Lattice Fermi operators satisfy the discretized version of the usual anti-commutation relations: As figure 1. indicates, the chirality of the Fermi operators is related to the parity of the lattice-site index. Namely, left-mover fields live on the odd-and right-mover fields live on the even-edges of the lattice, respectively: For later purposes it is worth to rewrite the lattice Fermi operators in terms of the spinoperators of the lattice. This can be achieved by a Jordan-Wigner transformation: where σ ± = 1 2 (σ x ± i σ y ) with σ x,y,z being the Pauli-matrices. The dynamics of the regularized model is given by light-cone evaluation operators: U L and U R . They are given by transfer matrices of an inhomogeneous 6-vertex model: where T is the trace of the monodromy matrix over the auxiliary space V 0 ≃ C 2 , The monodromy matrix is given by the R-matrix of the 6-vertex model in the usual way [36], , (2.14) such that ξ n s denote the inhomogeneities of the model. The entries of the 2 × 2 monodromy matrix act on the quantum space of the model H = ⊗ N i=1 V i with V i ≃ C 2 and they play crucial role in the algebraic Bethe-Ansatz solution of the model. In (2.12) H, P and a denote the Hamiltonian, the momentum and the lattice constant of the model, respectively. In order to get a massive interacting quantum field theory as the continuum limit of this lattice construction, the inhomogeneities of the vertex-model must be chosen as follows [20]: where the parameter ρ 0 is a-dependent and plays crucial role in the regularization scheme: Here M denotes the physical mass of fermions (solitons), L stands for the volume and N is the number 2 of lattice sites of the 6-vertex model. The parameters of the regularized lattice model are the inhomogeneities, the number of lattice sites and the anisotropy parameter γ. In (2.16) and (2.17) we described how to choose the inhomogeneities to obtain a massive interacting integrable quantum field theory in the continuum limit. The large N, ρ 0 = fix solution of the 6-vertex model shows [21] that this massive continuum quantum field theory is nothing but the MT or SG model, provided the following relation holds between the anisotropy parameter of the vertex model and the coupling constants of the Lagrangians (2.1) and (2.3): For later purpose it is worth to introduce a new parametrization for the anisotropy parameter: γ = π p+1 , then: We note, that the regimes 0 < p < 1 and 1 < p correspond to the attractive and repulsive regimes of the model, respectively. The definition (2.12) embeds the light-cone evolution operators of our model into the hierarchy of mutually commuting set of transfer matrices of the 6-vertex model. This implies that the Hamiltonian and the momentum of the model can be diagonalized via the Algebraic Bethe Ansatz method [33].

Algebraic Bethe Ansatz
In the framework of algebraic Bethe Ansatz, the eigenvectors of the transfer matrix (2.13) are constructed by successive application of creation operators on the bare vacuum of the model. The bare vacuum or reference state |0 is the completely ferromagnetic state with all spins up. The role of creation operators are played by the 12-matrix element of the monodromy matrix (2.14): T 12 (λ) = B(λ) which decreases the S z quantum number of a state by 1. A state constructed in this manner: is an is an eigenvector of the transfer matrix, provided the spectral parameters in the arguments of the creation operators satisfy the Bethe equations: In the Algebraic Bethe Ansatz approach the solutions of the Bethe equations play central role, since all physical quantities, such as the eigenvalues of transfer matrices can be expressed in terms of these roots: For the cases, when the number of Bethe-roots is large, it is more convenient to reformulate the Bethe-equations (2.21) in their logarithmic form. The central object of this formulation is the so-called counting-function given by the formula [26]: where φ ν (λ) is an odd function on the whole complex plane with all discontinuities running parallel to the real axis, such that in its fundamental domain, it is given by the analytic formula: The counting-function allows one to reformulate the Bethe-equations (2.21) in the form as follows: In this formulation, depending on the value of δ, an integer or half-integer quantum number I a can be assigned to each Bethe-root. When one considers states formed by only real Betheroots, then all these quantum numbers are different 3 and they characterize the state uniquely. The true vacuum corresponding to the ground state of the quantum field theory, is the S z = 0, anti-ferromagnetic vacuum with δ = 0. This state is formed by real Bethe-roots such that the quantum numbers of the Bethe-roots fill completely the whole allowed range [Z λ (−∞)/2π, Z λ (∞)/2π] . The excitations above this sea of real roots are characterized by complex Bethe-roots and holes. In this paper we will consider only hole excitations, since they correspond to fermion or soliton excitations of the continuum quantum field theory. The holes are such special real solutions of (2.25), which are not Bethe-roots 4 . Holes can be interpreted as missing Bethe-roots and the quantum numbers of the missing Bethe-roots can be assigned to them: where h k denotes the positions of the holes and their number is denoted by m H .

NLIE for the finite volume spectrum
When one has to deal with a large number of Bethe-roots, it is worth to rephrase the Betheequations (2.21) in a form of a set of nonlinear-integral equations (NLIE) [23]- [29].
Here we present the equations only for the pure hole sector of the theory [25] and here we will use the rapidity convention for the equations. This means a simple rescaling of the spectral parameter: θ = π γ λ. In the pure hole sector the counting-function in rapidity variable Z N (θ) = Z λ ( γ π θ) satisfy the nonlinear-integral equations as follows: where χ(θ) is the soliton-soliton scattering phase and G(θ) denotes its derivative: , (2.28) 0 < η < min(pπ, π) is an arbitrary positive contour-integral parameter, Θ = ln 2 N M L is the inhomogeneity parameter of the vertex-model and H k = π γ h k denote the positions of the holes in the rapidity convention. They are subjected to the quantization equations: The nonlinearity of the equations is encoded into the form of L The number of holes is not independent of the S z quantum number of the state. The connection between these two quantum numbers is given by the counting-equation 5 [26]: where here [...] stands for integer part. This equation immediately implies that on a lattice with even number of sites, only states with even number of holes can exist.
The main advantage of formulating the spectral problem in terms of the counting function is that it has a well-defined continuum limit. If one keeps the hole quantum numbers fixed, it is just the N → ∞ limit of the lattice counting-function [23,24]: The continuum counting-function satisfy the nonlinear-integral equations equations as follows [25]- [29]: where ℓ = M L with L being the volume and M is the fermion (soliton) mass. In the continuum limit the energy and momentum of the pure hole states read as: The counting-equation also changes non-trivially in the continuum limit [17,26]: where Q is the U (1) (topological) charge of the continuum model. The large volume limit of equations (2.34), (2.35) and (2.36) implies that the positions of the holes correspond to the rapidities of the fermions (solitons). This is why in the sequel we will refer to holes as fermions or solitons.
Finally, we note that the actual value 6 of the quantum number δ is important from the point of view of the continuum theory. Its value can make difference between fermions (δ = 1) of the MT model and the solitons (δ = 0) of the SG model in the odd U (1) charge sector of the theory [27]- [29]. In the even charge sector only the δ = 0 value is physical and there is no difference between MT fermions and SG solitons [27]- [29].
From the discussion above it follows that only the even charge sector of the MT and SG models can be regularized by the twistless 6-vertex model. The description of the odd charge sector requires a twisted vertex-model with a twist angle ω = π 2 [59]. However, in this paper we restrict ourselves to the twistless case.

Expectation values of the trace of the stress-energy tensor
In this subsection using the NLIE description of the finite volume spectrum given by (2.33) and (2.34), we compute the fermionic (solitonic) expectation values of the trace of the stressenergy tensor, which will be denoted by Θ T . The finite temperature 1-point functions, which correspond to the finite volume vacuum expectation value, has been previously computed and discussed in [66] and [60].
It has been shown in [18] that the diagonal matrix elements of Θ T can be computed from the volume dependence of the energy of the sandwiching state by the following formula: In our case E(ℓ) is given by (2.34) and thus can be computed from the NLIE equations (2.33). As a starting point, it is worth to compute the infinite volume or in other words the bulk expectation value: Θ ∞ T . Using Zamolodchikov's argument [18], it can be expressed in terms of the eigenstate independent bulk energy of the model by the formula: In the MT (SG) model the bulk energy term is of the form [24]: Inserting (2.39) into (2.38) one obtains: As a next step we express the non-bulk part of Θ T in (2.37) in terms of the solution of the NLIE (2.33). To do so, it is worth to introduce some useful notations. Let F ± (θ) denote the nonlinear combinations as follows: (2.41) 6 On the lattice the actual value of δ can be influenced by the parity of N 2 .
Then the derivative of L ± (θ) with respect to any parameter P is given by the formula: In practice P can denote one of the parameters of the NLIE equations (2.33). Namely, it can be the dimensionless volume ℓ, the spectral parameter θ or one of the positions of the holes H j . The second term in the right hand side of (2.37) consists of two terms. The first term ∼ E(ℓ) ℓ can be expressed in terms of dZ(θ) dθ and of F ± (θ), while the second term ∼ dE(ℓ) dℓ turns out to be the functional of dZ(θ) dℓ and of F ± (θ). Integrating the right hand side of (2.34) by parts, E(ℓ) ℓ can be rephrased as follows: They satisfy the set of linear integral equations as follows: Taking the derivative of (2.34) and (2.33) with respect to ℓ, leads the following expression for . They are solutions of the set of linear integral equations as follows:   .47) gives the required expectation value. Though, this representation for Θ T might seem strange for the first sight, but in the later sections it will turn out, that it fits very well for the structure of the lattice results.
In the rest of the paper our main goal is to reproduce the formula (2.47) from the lightcone lattice computation of the multi-fermion (soliton) expectation values ofΨΨ.

The lattice counterpart ofΨΨ
We close this section with a short discussion about the lattice counterpart of the operator ΨΨ in the MT model. A simple Jordan-Wigner transformation (2.11) shows, that certain bilinears of the lattice Fermi operators are simple expressions of the lattice spin operators: where σ ± n are the usual spin creation and annihilation operators corresponding to the nth site of the lattice.
Using the representation (2.2) for γ 0 , the following identification can be made for the unrenormalized bare operators on the lattice: where the term 1 a is introduced to account for the correct bare dimension of the continuum Fermi field.
A similar computation shows that the pseudo-scalar combination of the Fermi operators correspond to the antisymmetric combination of the spin operators: Thus (2.49) and (2.50) implies, that the determination of the expectation values of the bare scalar-and pseudo-scalar fermion bilinears is equivalent to computing the two-point functions of neighboring spin operators. This task is completed in the rest of the paper via the QISM [33]- [57].
We note that beyond the computation of 2-point functions σ ± n σ ∓ n+1 the 2-point function e n e n+1 with e n = 1 2 (1 n − σ z n ) can also be computed with the techniques presented in this paper. This later 2-point function contains a combination a 4-fermion term, as well: e n e n+1 = (ψ + n ψ n − 1 2 )(ψ + n+1 ψ n+1 − 1 2 ). A usual argument based on the bare dimensions of the operators implies that this operator has a nontrivial mixing under renormalization. This means that the correct implementation of the renormalization process requires the computation of the expectation values of further operators. This investigation is left for future work.

Computation of lattice correlators
The strategy of computing the fermionic (solitonic) expectation values ofΨΨ consists of two main steps. First, one has to compute the expectation values of the lattice operators σ ± n σ ∓ n+1 and then one has to take the continuum limit of the lattice results by sending the number of lattice sites N to infinity. In [17] the efficiency of this method has been demonstrated via the computation of the solitonic expectation values the U (1) current of the model. In this section we describe in detail the computation of our special lattice spin-spin 2-point functions.
Consider a vector of the Hilbert-space obtained by successive action of creation operators on the bare vacuum: Such a state is called Bethe-state if the numbers λ j are arbitrary and it is called Betheeigenstate if the set {λ j } j=1,..,m is equal to the set of roots of the Bethe equations (2.21). The corresponding "bra" vector can be defined by acting from the right with the annihilation operators on the "bra" bare vacuum: The determination of the diagonal form factors ofΨΨ andΨγ 5 Ψ requires the computation of the following two special 2-point functions: where here | λ denotes a Bethe-eigenstate. The determination of these 2-point functions can be achieved in a purely algebraic way [35,36] within the framework of the QISM [33], such that only the Yang-Baxter algebra relations and the expression of local spin operators in terms of the elements of the monodromy matrix [35] of the model are used.
The core of the algebraic computations is the relation between the local spin operators and the elements of the Yang-Baxter algebra [35]: where the operator E n is given in terms of spin operators as follows: The formulas (3.4) and (3.5) together with the Yang-Baxter algebra relations imply that any correlators of the spin operators can be represented as complicated linear combinations of scalar products of Bethe-states and Bethe-eigenstates.
The coefficients in such a linear combination follow from the Yang-Baxter algebra relations. For our actual computations we need only to know, how an operator B(ξ n ) with ξ n being an inhomogeneity of the vertex model, act on a "bra"-vector. This is given by the formula [36]: Here the functions f 0 , f 1 and f 2 are of the form: As a consequence of (3.8) and (2.21), r(λ) satisfies some important identities: where here λ k s are solutions of (2.21). We note that in general, in (3.6) and (3.7), λ k s are not necessarily the solutions of (2.21) and so M can be any positive integer number different from m, which denotes the number of Bethe-roots.
To carry out the computation of our special 2-point functions we need to know the scalar product of a Bethe-state and a Bethe-eigenstate. This is given by Slavnov's determinant formula [55]. Let | µ an arbitrary Bethe-state in the sense of (3.1) and | λ be a Betheeigenstate. Then their scalar product can be determined with the help of the formula [55]: where H( µ| λ) is an m × m matrix with entries: (3.11) An important special case of (3.10), when the scalar product of two identical Bethe-eigenstates are considered. This is called the Gaudin formula: where Φ( λ) is the Gaudin-matrix, which is related to the counting-function (2.23) by the formula: During the computation of the special 2-point functions considered in this work, such scalar products arise, in which the components of the vector µ take values either from the set of Bethe-roots {λ j } j=1,..m or from the set of inhomogeneities {ξ k } k=1,..N of the model. In these cases the matrix elements of H( µ| λ) remarkably simplify: (3.14) (3.15) These simplifications allow one to compute a typical scalar product arising in the computation of diagonal form factors: where µ (a 1 ,...,a K ) (ξ α 1 , .., ξ α K )| denotes a state, which is obtained from λ| = λ 1 , ...λ m | by replacing K pieces of λ j to certain inhomogeneities of the lattice model: such that both sets {a k } and {α k } contain distinct numbers. In (3.16) Y denotes a K × K matrix with entries as follows: where the m-component vector X(λ j |ξ) is the solution of a set of linear equations: In [17] it has been shown, that the discrete linear problem (3.19) can be transformed into a set of linear integral equations as follows. Let: with Z λ (λ) being the lattice counting-function in the λ-convention. 7 Then for pure hole states above the anti-ferromagnetic vacuum, G(λ|ξ) satisfies the linear integral equations as follows [17]: where h j s denote the positions of the holes, η is a small positive contour-integral parameter, and G λ (λ) is related to the kernel of NLIE equations (2.28) by: The discrete set of variables X j (ξ) in (3.22) are related to the continuous function G(λ|ξ) by the formula: The correlators of our interest are as follows: where here | λ denotes a Bethe-eigenstate. Using (2.22), (3.6) and (3.16) the following expressions can be obtained for these 2-point functions: In (3.29) and (3.30) we preserved the determinant structure implied by (3.16). Nevertheless for future computations it is better to shift the anti-symmetrization to the coefficient of the quadratic expression of X. Then the quadratic in X parts of (3.29) and (3.30) can be written as follows: (3.32) where f(λ, λ ′ |ξ) is an antisymmetric function given by the formula: In (3.31) and (3.32) the coefficient function f has better large λ and λ ′ asymptotics, than the coefficients of the quadratic terms of (3.29) and (3.30). This property proves to be very useful, when the discrete sums in (3.29) and (3.30) are transformed into integral expressions.
As it was done in [17], this transformation can be achieved by the lemma as follows [24,26]: .,m solutions of the Bethe-equations (2.21) and let f (λ) a meromorphic function, which is continuous and bounded on the real axis. Denote p (f ) its pole located the closest to the real axis. Then for |Im µ| < |Im p (f ) | the following equation holds:

34)
where h j and c j denote the positions of holes and complex Bethe-roots, respectively and F (λ) ± (λ) is given by (3.24). η is a small positive contour-integral parameter which should satisfy the inequalities: where p ± λ denotes those complex poles of F (λ) ± (λ), which are located the closest to the real axis. The validity of formula (3.34) in µ can be extended to the whole complex plane by an appropriate analytical continuation method [26].
The typical sums arising in (3.29), (3.30), (3.31) and (3.32) are of the form: where the functions f 1 , f 2 , f 3 and f ± are of the form: and we exploited the concrete inhomogeneity structure of the model, namely that ξ 2n = ρ 0 −i γ 2 and ξ 2n+1 = −ρ 0 − i γ 2 , with ρ 0 given by (2.17). Formula (3.39) together with (3.40)-(3.44) and the integral representations given in appendix A constitutes our final result for the lattice expectation values of the scalar and pseudo-scalar fermion bilinears. To get the expectation values of these operators in the continuum theory, the lattice formulas (3.39) should be evaluated in the continuum limit. This will be discussed in the next section.

Continuum limit
In this section the expectation values (3.39) are evaluated at the continuum limit. This task reduces to the evaluation of the sums (3.36) and (3.37) entering (3.39) in the large ρ 0 limit. The ρ 0 dependence of these terms is given by the ρ 0 dependence of the functions f 1 , f 2 , f 3 , f ± given in (3.40)-(3.44), and the ρ 0 dependence of X(h j |ξ ± ) and G(λ|ξ ± ). Latter is governed by the linear integral equation (3.22). First, it is worth to discuss the continuum limit of the variables X(h j |ξ ± ) and G(λ|ξ ± ). They are solutions of the set of equations (3.22)-(3.26). The continuum limit means, that one has to take the number of lattice sites N to infinity, such that the inhomogeneity parameter ρ 0 is tuned with N according to the formula (2.17). This means, that in the continuum limit procedure ρ 0 also tends to infinity, but logarithmically in N or a. The large N limit of G(λ|ξ ± ) is determined by the large ρ 0 expansion of S 0 (λ|ξ ± ) in (3.22): According to (2.17) e − π γ ρ 0 ∼ 1 N ∼ a. Since apart from S 0 , all terms in the equations (3.22)-(3.26) are proportional to G(λ|ξ) and X(h j |ξ), (4.1) implies that: As a consequence of (4.2), the terms being quadratic or multi-linear in G and X in (3.39) will be suppressed in the continuum limit by a factor of 1 N ∼ a with respect to the leading order. This implies that the leading order expression for O ± 2n λ can be represented as a linear functional of G(λ|ξ ± ) and X(h j |ξ ± ), and the quadratic terms become negligible in the continuum limit 8 .
To get the leading order term of O ± 2n λ in the continuum limit, first it is worth to eliminate the finite part of the 1 N term by the definitions: The finite parts G ± and X ± (h j ) satisfy the equations: , j = 1, ..., m H .

(4.4)
With help of the large ρ 0 expansions summarized in appendix B, one obtains the following leading order result for O ± 2n λ in the attractive regime: where It turns out, that the variables G d (θ), X (d) j of (2.44) and G ℓ (θ), X (ℓ) j of (2.46) can be expressed in terms of G (±) (λ) and X (±) j of (4.4) as follows: (4.8) Using (4.8) and changing variables from λ to θ in (4.5), one ends up with the final result: where dots mean next to leading order terms in the N tends to infinity limit of the attractive regime. Comparing (2.47) and (4.9) one can easily recognize the proportionality of O + 2n λ and Θ T : According to (2.49) the expectation value for the bare fermion bilinear is given by: where we exploited the relation (2.17) between the lattice constant a and the inhomogeneity parameter ρ 0 . Using the relation between p and β in (2.19), one can see that Ψ Ψ is proportional to the expectation value of the stress energy tensor, and it scales as a β 2 /4π−1 as it is expected from (2.8) obtained via purely field theoretical considerations.

Large volume expansion
In this section we rephrase the leading order in N terms of (4.9) and (4.10) in the form of a systematic large volume series. To get rid of the unnecessary constants, we consider the following quantities: where G d (θ), X (d) and G ℓ (θ), X (ℓ) are defined by the equations (2.44) and (2.46) respectively. From (2.47) and (4.9) it can be seen that O + is simply related to the fermionic expectation value of the trace of the stress-energy tensor: On the other hand O − is proportional to the fermionic expectation value of the renormalized pseudo-scalar fermion bilinear Ψ γ 5 Ψ . Here we do not care about the actual value of the proportionality factor, since it will turn out, that this expectation value is zero between multi-fermion states. The process of the evaluation of (5.1) and (5.2) in the large volume limit is very similar to the method used for computing the diagonal matrix elements of the trace of the stress-energy tensor in purely elastic scattering theories [11]. The reason for this is that formally the NLIE equations (2.33) are very similar to the Thermodynamic Bethe Ansatz (TBA) equations of a purely diagonally scattering theory of two types of particles. This analogy, the actual form of the large volume series of the U (1) current of the theory [17] together with the all order large volume series conjectures for the diagonal form factors of purely elastic scattering theories [10,11,[14][15][16], led to the following large volume series conjecture for the diagonal multifermion (soliton) expectation values of local operators in the MT (SG) models [17]:

4) is called the dressed form-factor [11] and it is given by an infinite sum in terms of the connected diagonal form-factors of the theory:
where F O c denotes the connected diagonal multi-fermion (soliton) form factors of O(x) in pure fermion (soliton) states and F ± (θ) are the nonlinear expressions of the counting function given by (2.41). We note that the structure (5.4) is the same for the purely elastic scattering theories and for the MT (SG) model. The difference arises in the concrete form of the exact Gaudinmatrix 9 and in the actual form of the dressed form factors. However up to exponentially small in volume corrections the formulas of purely elastic scattering theories are also appropriate to describe the multi-fermion (soliton) expectation values of local operators [9], [17].
So far, conjecture (5.4)-(5.7) has been checked against the diagonal fermionic (solitonic) form factors of the U (1) current of the theory [17] and now by rephrasing O + as a large volume series, we will show that this conjecture remains valid in the case of the trace of the stress energy tensor, too. Thus, our purpose is to bring O ± into the form of (5.4) and check whether the coefficients of ρ({H − }|{H + }) agrees with D O ({H + }) given by (5.7).
The first step of the computation is to rewrite G d (θ), X (d) and G ℓ (θ) X (ℓ) in terms of the solutions of some "elementary" linear problems. For any function f, let f [±] (θ) = f (θ ± iη), then by definition an "elementary" solution indexed by A satisfy the linear equations as follows: A (θ), α = ±, (5.8) where the symmetric kernel ψ αβ (θ) is given by: and f A (θ) is the source term specifying G

[α]
A (θ). An "elementary" solution with unshifted argument satisfies the equations as follows: From the defining linear equations (5.8) it can be shown, that the "elementary" solutions satisfy the following identities: In the linear problems (2.44), (2.46) and (3.22) belonging to the real physical problem, there are X-vectors associated to the continuous functions denoted by G ... (θ). This is why, it is also convenient to introduce an X-vector to each elementary solution by the analogous definitions as follows: Then using the linear equations (2.44), (2.46), (3.22) and (5.8) the pairs of quantities G d (θ), X (d) and G ℓ (θ), X (ℓ) can be expressed in terms of the elementary solutions as follows: whereΦ kj ( H) is the exact Gaudin-matrix defined by the formula: Using the formulas aboveΦ kj ( H) and O ± can be expressed in terms of the elementary solutions of (5.8) as follows: (5.20) To bring (5.19) and (5.20) into the form of (5.4) one needs to use two theorems. Theorem 1. The inverse of the Gaudin-matrix can be expressed in terms of its principal minors and sequences of its matrix elements [58], [11] by the formula as follows: with C ij being the co-factor matrix with entries:  , α 1 , ...,α n }).
(5.23) Theorems 1. and 2. allow one to rewrite the double and single sums respectively in (5.19) and (5.20) into a more convenient form. Using (5.23), the single sums of (5.19) and (5.20) can be represented as: with A, B ∈ {s, c}. It is convenient to determine first the large volume series expansion of the first terms in the right hand sides of (5.19) and (5.20). They are called the vacuum contributions [10,11] since they correspond to the {H + } = ∅ case. These terms can be rephrased as an infinite series similar to that of LeClair and Mussardo [14][15][16]. To get this series representation, first one has to construct the all order large volume solution of (5.8) for A = s, c. This can be obtained by a simple iterative solution of the equations. Then one has to insert these large volume series into (5.19), (5.20). At the end of this process one gets a bulky sum of terms, such that a lot of terms are identical under certain permutations of the integrating variables. Taking into account these permutational symmetries by appropriate symmetry factors, one ends up with the formula for the vacuum contributions as follows: n + +n − ,c (θ 1 +i η, ..., θ n + +i η, θ n + +1 −i η, ..., θ n + +n − −i η), where the "connected" form factors of O ± are given by the definitions:
one 10 can conclude that the fermionic expectation values ofΨγ 5 Ψ are zero. Next we discuss the result obtained for O + . According to (5.3), O + is proportional to the fermionic (solitonic) expectation values of Θ T . This operator belongs to a conserved current, this is why its connected form-factors between pure fermion (soliton) states can be determined by using the arguments of references [14] and [15]. The actual computations lead to the following simple result: Then, (5.31) together with (5.3), (5.29) and (5.30) imply that the conjecture (5.4) of ref. [17] is valid also for the diagonal fermionic (solitonic) matrix elements of the trace of the stress-energy tensor.

Summary, outlook and discussion
In this paper, using the light-cone lattice regularization, we computed the finite volume expectation values of the composite operatorsΨΨ andΨγ 5 Ψ between multi-fermion (soliton) states in the Massive Thirring (sine-Gordon) model. In the light-cone regularized picture, these expectation values are related to such 2-point functions of the lattice spin operators, in which the operators are located at neighboring sites of the lattice. The operatorΨΨ is particularly interesting operator, since it is proportional to the trace of the stress-energy tensor. Thus, its continuum finite volume expectation values can be computed [18] also from the set of non-linear integral equations (NLIE) governing the finite volume spectrum of the theory. The final result, which was obtained after a lengthy computation of the spin-spin 2-point functions of neighboring operators, reproduced the pure NLIE result.
In general the finite volume matrix elements of local operators are computed via their large volume series. Thinking in this framework, previously in [17], an all order large volume series representation similar to that of [10,11,14] was conjectured for the finite volume diagonal matrix elements of local operators between multi-fermion (soliton) states. To check the conjecture of [17] we rephrased the diagonal multi-fermion (soliton) matrix elements of the trace of the stress-energy tensor as a large volume series. The form of the series was conform to the conjecture of [17].
Nevertheless, one has to note that, so far the large volume series conjecture of [17] for the finite volume diagonal matrix elements of local operators between pure fermion (soliton) states, has been checked in cases when the local operator belongs to a conserved quantity of the theory. This is why it would be interesting to test the conjecture for such operators, which do not belong to the conserved quantities of the model. The results of [7,8] indicate that the truncated conformal space approach could be an appropriate method for such investigations.
Beyond the results of [17] and the present paper, several questions are still open. The light-cone lattice approach gives access to all eigenstates of the Hamiltonian. Thus, diagonal matrix elements between non-pure fermion or soliton states can also be computed in this framework. It would be interesting to see, how the large volume series representation of [17] should be modified in that case. Finally, a much more difficult but still open problem is the determination of non-diagonal finite volume form factors.
Beyond our approach, there is another approach to the form-factors of the SG model in cylindrical geometry [61][62][63][64][65][66][67]. We close the paper with some discussion concerning the comparison of our method with that of the series of papers [61][62][63][64][65][66][67]. In [61][62][63][64][65][66][67] the hidden Grassmannian structure of the XXZ model was exploited to determine the finite temperature 1-point functions [66] and infinite volume form-factors [67] of the local operators of the SG theory. In this approach the compactified direction is time and the compactification length corresponds to the inverse temperature. The 1-point functions and the form-factors are computed as the continuum limit of appropriate partition functions of the 6-vertex model.
Our approach is more conventional. We consider the inhomogeneous 6-vertex model as a lattice regularization of the MT (SG) model and the operators we consider is the set of composite operators of Fermi fields and their derivatives 11 . In our case the compactified direction is space which makes easy to consider matrix elements of operators between excited states of the model. The Fermi fields are expressed in terms of local spin operators and the form-factors are given by such correlation functions of the spins, in which the spins are located at neighboring positions on the lattice. The correlators are evaluated by usual Algebraic Bethe Ansatz methods [35]- [39].

Acknowledgments
The author would like to thank Zoltán Bajnok and János Balog for useful discussions. This work was supported by OTKA grant under K109312. The author also would like to thank the support of an MTA-Lendület Grant.

A Integral representation of some typical sums
In this appendix we summarize the integral representation of the typical sums (3.36), (3.37) arising in the computation of the expectation values (3.27) and (3.28). For the sake of completeness we recall them in this appendix, too: In order to be able to transform these sums into integral expressions, we require that in each of its arguments f should not have worse than constant asymptotics at infinity. In (A.2) we also require for f to be an antisymmetric function of its arguments. Using the lemma (3.34) and the linear integral equations (3.22) the following formulas can be derived for Σ (1) where J 0 [f ](ξ) and J G [f ](λ) are functionals of f and are of the form: such that ⋆ denotes the convolution with conventions given by (B.3). The formula for Σ (2) λ [f ](ξ, ξ ′ ) is composed of six terms: The lower index of each term on the right hand side of (A.5) refers to the internal structure of the expression as it becomes clear from their explicit form: where the two functions F X and F XX , which are functionals of f, take the form: with the "elementary" functionals: We note, that at the derivation of (A.3) and (A.5), it was important to bring the sums into a sum of such integral expressions which contain G(λ) only in the combination G(λ) F α (λ). The reason for the preference of such a form is, that this combination of variables has a well defined continuum limit. This convenient form could be derived by eliminating the single G(λ) terms with the help of (3.22).

B Large ρ 0 expansions
In this appendix we summarize, how one can obtain the kernels K ± (λ|ρ 0 ) and the bulk term of (4.5) via the computation of the large ρ 0 (2.17) limit of the functionals entering (3.39).
The key point of the computations is that one should work in Fourier space. This is why as a first step we fix our conventions for the Fourier-transformations. The Fourier-transform of a function f is given by:f The inverse transformation reads as: The Fourier-transform of the convolution of two functions f and g is given by the product of individual Fourier-transforms 12 : 12 Provided they exist.
Using (3.39), (A.3) and (A.5) the bulk term of (4.5) can be represented as follows: where here the upper index ± refers to the upper index of the operators O ± 2n . The other important terms are the coefficients of X(h j |ξ ± ) in (3.39). They are denoted by C (±) (h j |ξ ± ), such that the upper index of C (±) (h j |ξ ± ) refers to the upper index of the operators O ± 2n , and the lower index agrees with the index of ξ ± corresponding to the actual value of the inhomogeneity parameter in the 2nd argument of X(h j |ξ). With the help of (A.3) and (A.5) one can obtain their form: Then the following kernels can be defined: The bulky computations detailed below show that in the large ρ 0 limit the upper index of K where dots stand for terms tending to zero, when ρ 0 → ∞ and p < 1, and K ± (λ|ρ 0 ) is given by (4.6). Now, we compute K (±) ± (λ|ρ 0 ) of (B.7) and C ± bulk of (B.4) in the large ρ 0 limit and derive the formulas (4.6) and (4.7). To complete this task, the following functionals should be computed in the large ρ 0 limit: • J G [f ](λ) of (A.4) taken at the functions f 1 , f 2 , and f 3 given by (3.40), (3.41) and (3.42).
We note that in principle one should compute F XX [f ± ](λ, λ ′ ), too. Nevertheless, at the end of this appendix we will argue, that we do not need to compute these functionals explicitly. The reason is that, as implied by (A.10)-(A.11), F XX [f ± ](λ, λ ′ ) arise together with multi-linear or quadratic combinations of G(λ|ξ ± ) and of X(λ|ξ ± ) and these combinations always become next to leading order in the continuum limit. The strategy of the large ρ 0 evaluation of the above listed functionals is as follows. One can write them in a very special form, namely as a linear combination of convolutions, such that ρ 0 appears in the argument of the convolutions. This means, that the large ρ 0 expansion of the functionals of our interest is equivalent to determine the large argument series expansion of the convolutions appearing in them. It is worth to represent these convolutions in Fourierspace, since by using the property (B.3) they can be written as a single Fourier-integral. The large argument series expansion of these Fourier-integral expressions can be computed by using the residue theorem. Thus the positions of the poles of the integrand will determine the large argument decay of the convolutions of our interest.
Let us start with the computation of the bulk (B.4) or in other words the global constant in rapidity term. The representation of the building blocks of (B.4) as linear combinations of convolutions read as: where T 0 (ρ 0 ) and I(ρ 0 ) are given by: The large argument expansions of appendix C lead to the following leading order large ρ 0 expressions for these building blocks in the attractive regime: . Now, we can continue with the computation of the kernels (B.7). These kernels can also be represented as appropriate linear combinations of some convolutions. This type of representations of the elementary building blocks of (B.7) are given by the formulas: The large argument expansions presented in appendix C lead to the following leading order large ρ 0 forms for the functionals (B.30)-(B.36): of (3.39) contains the multi-linear and quadratic in G and X terms of present interest. Using (A.5), (A.9), (A.10) and (A.11), it is straightforward to show that the multi-linear and quadratic terms in G and X are given by the terms as follows 13 : From (A.9), (A.10) and (A.11) and from the magnitude of f ± (λ, λ ′ ) it follows that the numerator of (B.43) is of order e −2 p ρ 0 . Then putting everything together it turns out that multi-linear and quadratic in G and X terms are of order: e −(1+p)ρ 0 ∼ a, which tend to zero in the continuum limit.