The LeClair-Mussardo series and nested Bethe Ansatz

We consider correlation functions in one dimensional quantum integrable models related to the algebra symmetries $\mathfrak{gl}(2|1)$ and $\mathfrak{gl}(3)$. Using the algebraic Bethe Ansatz approach we develop an expansion theorem, which leads to an infinite integral series in the thermodynamic limit. The series is the generalization of the LeClair-Mussardo series to nested Bethe Ansatz systems, and it is applicable both to one-point and two-point functions. As an example we consider the ground state density-density correlator in the Gaudin-Yang model of spin-1/2 Fermi particles. Explicit formulas are presented in a special large coupling and large imbalance limit.

Different methods have been developed for integrable relativistic Quantum Field Theories. In these models the form factors are obtained as solutions to the so-called form factor bootstrap program [31][32][33][34]. Originally the form factor approach was developed to compute two-point functions at zero temperature, but later it was found that even finite temperature correlations can be obtained with it. A key development was an integral series derived for finite temperature one-point functions in [35], which is today known as the LeClair-Mussardo (LM) series. The validity of the LM series for two-point functions was questioned first in [36,37], and it was pointed out in [38] that it is not even well defined. Later a well defined prescription for the two-point functions was found in [39]. Interestingly, in the case of the 1D Bose gas the LM series was discovered already in the very early works [4,5], which treated static two-point functions. These early results were partially forgotten and not discussed in most of the papers dealing with the LM series.
Most of the theoretical work (both on the spin chains and in field theory) focused on ground state and finite temperature correlations. However, the last 10 years saw tremendous progress in the field of non-equilibrium dynamics and quantum quenches in particular, and this led to the study of Generalized Gibbs Ensembles (GGE's), which arise as steady states in these systems [40,41]. It is now generally accepted that in Bethe Ansatz solvable systems the information stored in a GGE is equivalent to specifying the rapidity distributions of the quasi-particles [42], and this motivated the study of correlation functions in Bethe states with arbitrary root distributions [43,44]. Such studies were further spurred by the development of Generalized Hydrodynamics (GHD) [45,46], which describes the large scale transport properties of integrable models. A central assumption in GHD is the existence of fluid cells such that each cell is described by a space and time dependent GGE, thus the rapidity distributions can become complicated dynamical functions. An especially interesting result for dynamical two point functions was derived in [47], using the GHD and the LM series.
Despite all of the progress described above, most of these works are concerned with models of one-component gases or spin chains related to gl(2) symmetry. Much less is known about the correlation functions of spin chains with gl(N) spin symmetry and multi-component gases. There are partial results for short-range correlation functions in higher rank spin chains [48][49][50]. Furthermore, form factors of local operators were computed in models related to the algebra symmetry gl(3) and gl(2|1) in [51][52][53][54] using the nested ABA approach. This allowed to treat the asymptotic behavior of correlation function using form factor expansion [55]. On the other hand, up to now there have been no results for generic correlation functions at small or intermediate distances, and it was not known whether some generalization of the LM series exists in the nested models.
Our goal is to contribute to the solution of these problems. We develop a form factor expansion in nested Bethe Ansatz systems, which can be applied to one-point and (static) two-point functions in arbitrary equilibrium ensembles. Thus it can be used both for the ground state, in Gibbs and Generalized Gibbs Ensembles, but also in non-equilibrium situations where the GHD can be applied. Our result is the generalization of the LeClair-Mussardo formula to nested non-relativistic models. We use standard tools and also recent developments of the algebraic Bethe Ansatz to derive the LM series. We also explain how to compute the so-called symmetric form factors of the two-point functions, which are the key ingredients of the LM series.
The method is applicable for an arbitrary model related to algebra symmetry gl(2|1) or gl (3).
We focus on one specific example, namely the 1D spin-1/2 Fermi gas (a.k.a. Gaudin-Yang model, which is related to the algebra symmetry gl(2|1)). A similar example of the two-component 1D Bose gas will be considered separately. The paper is organized as follows. In Section 2 the main results of the paper are described, omitting the technical details in the various definitions. Notation and short description of Bethe Ansatz are given. In section 3 we introduce the symmetric form factors and prove the finite volume expansion theorem for algebra symmetry gl(2|1) related models. The case of algebra symmetry gl(3) related models can be obtain in a similar way and the difference is minor between these cases. In Section 4 the thermodynamic limit of the form factor expansion is considered. In Section 5 we present an efficient way for the computation of the symmetric form factors, focusing on the Gaudin-Yang model. Finally, this Section includes our explicit results for the correlators in the combined large coupling and large imbalance limits. Appendix A includes technical details of the computations involving the symmetric form factors. In Appendix B we consider the mean values and form factors of the conserved charges, and we perform a check of the expansion theorem. Finally, in Appendix C the comparison of the the integral series with result of [3] for the correlation functions is given, considering the example of 1D Bose gas (Lieb-Liniger).

Main results
In this Section we summaries our main results and apply them to the one dimensional two component quantum gas models. Most of our results for correlation functions are valid in a much more general setting: they can be applied in Yang-Baxter integrable models related to the symmetry algebras gl(3) and gl(2|1). This will be explained in later Sections. However, in this summary we focus on the concrete applications.
Similarly, we can also consider the two component Bose gas defined through a Hamiltonian which has a form identical to (2.1) but with bosonic fields ψ ↑ (x), ψ ↓ (x) satisfying commutation relations ψ † α (x), ψ β (y) = δ αβ δ(x − y). The bosonic model is related to the algebra gl(3). In both cases periodic boundary conditions are assumed.
Both models can be solved by the nested Bethe Ansatz [10,11,57,58,61]. Eigenstates can be characterized by two sets of rapiditiesv = {v 1 , . . . , v b } andū = {u 1 , . . . , u a }. The first set of the rapidities describes the momentum of the particles, and the second set describes the orientation of the wave function in the internal spin space.
Let us however note that in further Sections we call the spin variablesū the first level and rapidity variablesv the second level. The only consequence of our choice will be the unusual convention for the numeration of levels. To summarize our conventions: the first setū describes the spin degrees of freedom, and the second setv describes the physical momenta of the particles.
We will focus mostly on the fermionic case. In this model the rapidities satisfy the Bethe equations [61]: where S, S 12 are phases of (quasi)particles scattering and p(v) is the single particle momentum, which in this model is simply p(v) = v. For the Bethe equations in the bosonic case see [10,57].
In this work we are interested in static two point correlation functions in these models. We focus on the density-density correlation functions and similarly in the bosonic case. Furthermore, in the actual computations we will consider the correlators where q(x) = q ↑ (x) + q ↓ (x) is the total density operator. In the following we present expansion theorems for the mean values in both finite and infinite volumes; in order to simplify the notations we will denote the operator product as O(x, y). Specific details of the local operators will enter the form factors involved.
Our first main result is the following finite volume expansion theorem: Here O(x, y) denotes the normalized expectation value and the sum on the r.h.s. is taken over all bipartite partitionsū → {ū I ,ū II },v → {v I ,v II } 1 , #ū = a, #v=b. The matrix G is the so-called Gaudin matrix, which is the Jacobian of the Bethe equations of the system with some finite size L. The cardinalities of the sets are #ū I = m, #v I = n and #ū II = a − m, #v II = b − n. The object F m,n is the so-called symmetric form factor of the composite operator, which is a special limit of an off-diagonal form factor defined as Precise details of this definition, including the normalization N are presented in Section 3.3. We should note that in a particular model the symmetric form factors F a,b can be non-zero even if there is no Bethe state corresponding to the particle numbers a, b. This happens because the symmetric form factor describes the amplitude of a multi-particle process in the presence of the other particles, and thus the selection rules for them can be different from those for the states. This type of finite volume form factors expansion was first used in [62] for generic local operators in integrable QFT. It was later shown in [43] that in the thermodynamic limit this expansion reproduces the so-called LeClair-Mussardo series [35], which is an integral series for the mean values. It was argued in [39] that in QFT the expansion theorem should hold also for composite objects such as the two-point functions. We should also note that in spin chains related to algebra symmetry gl(2) the theorem was proven rather generally in [63,64].
The expansion theorem treats the operator product as a single object, and the form factors depend on the operators and the distance x − y. In this respect our method is different from the direct approaches, where a complete set of states is inserted between the two operators. The advantage of this method is that in the thermodynamic limit it leads to integral series involving form factors with low particle numbers, which can be determined from relatively simple computations.
The resulting integral series can be considered the LeClair-Mussardo series for nested Bethe Ansatz models, which we describe in the following. Let us consider for simplicity a situation where in ground state the solutions to the Bethe equations are purely real. Extensions to the general case with so-called string solutions will be considered elsewhere. We take the thermodynamic limit of the formula (2.6) such that the particle densities are kept finite. Then we obtain the following result: where ω (1,2) are weight functions defined in (4.10). This series can be considered as the generalization of the LeClair-Mussardo formula, and it is proven in Section 4.3. The series involves the "symmetric" diagonal form factors. In contrast, the original LeClair-Mussardo formula used a different prescription for the diagonal limit, which are called the "connected" diagonal form factors. Furthermore, the original LM series does not involve the additional ω weight functions. The equivalence of the two integral series follows from relations between the two different evaluations of the diagonal limit, as discussed in length in [62]. We should note that in the one-component case the LM series for two-point functions was already derived in [4,5], and these works also used the connected form factors.
We also present the LM series in the nested case using the connected form factors, which is written as Here F c m,n (t;s) are the connected form factors, defined as Precise details of this definition are presented in Section 3.4. Using (2.6) with (2.7) the static zero-temperature two-point correlation functions can be calculated. We apply them to a specific large imbalance limit of the two-component fermionic model, which we describe in the following.
In the two-component Fermi gas with c > 0 there are two natural small parameters, which can be chosen as Q/c and B/c, where Q and B are the Fermi rapidities in the ground state, which involves only real roots. These two parameters can be related to the particle densities, the exact formulas are presented later in Section 4. Smallness of Q/c means large coupling limit, whereas smallness of B/c is equivalent to large imbalance, or smallness off the ratio n ↓ /n, where n = q(x) , n ↑,↓ = q ↑,↓ (x) . If the parameter B/c is small, then the densities can be expressed as (2.11) Note that these formulas are exact in Q/c. However in the case of two-point correlation functions the series (2.8) converges only when the ratio Q/c is small too. Thus we will develop a Taylor series in the two small parameters. Our result for the leading terms of the correlation functions is These leading terms describe free fermionic behavior in both cases. We expect that this behavior only holds at the particular level of approximation. Note that the correction terms to the formula (2.12) are quadratic; the absence of linear terms (carrying an additional factor of B/c or Q/c) is a non-trivial result of our computations.
In the x → 0 limit both formulas above approach zero. This agrees with result of Section 5.6 where leading terms of the local correlators were found from the Hellmann-Feynman theorem. It is shown there that Clearly, this is sub-leading compared to the result (2.13).

Notations
In order to simplify the formulas we use the following notational conventions in the paper. We denote sets of arbitrary variables asū = {u 1 , u 2 , . . . },v = {v 1 , v 2 , . . . }. Subsets are labeled by roman or Arabic numbers, e.g.v I ,v II ,ū 1 ,ū 2 , etc. We use the notationū j =ū \ u j , u j,k =ū \ {u j , u k } for sets with the certain elements are omitting. For an arbitrary function F (x), F (x, y) and for an arbitrary setsx,ȳ we use the notation In a similar way Skew-symmetric products of an arbitrary function g(x, y) over the setx are defined as

Bethe Ansatz
We consider models with algebra symmetry gl(3) and gl(2|1) simultaneously. We present detailed computations mostly in the case of gl(2|1); the case of algebra symmetry gl (3) is very similar and we will just give the final results. Throughout the paper we use conventional notations for the following functions: . (2.18) We will also use the kernel K(x, y) defined as We define a Z 2 -graded vector space C 2|1 with a grading Our computations will be performed using the Algebraic Bethe Ansatz (ABA) technique. In ABA the system is given by certain monodromy matrix T that satisfies the so-called RTT relation [9,11,65,66] with an appropriate R-matrix that holds in a tensor product of spaces V 1 ⊗ V 2 ⊗ H (subscripts of R 12 and T 0k denote spaces in which correspondent matrix acts non-trivially). Here the so-called auxiliary spaces V 1 and V 2 are C 3 in case of algebra symmetry gl(3) and C 2|1 in case of algebra symmetry gl(2|1) related models. So-called quantum space H coincides with the Hilbert space of the Hamiltonian of considering system. The cases of rational gl(3) and gl(2|1)-invariant R-matrices are considered in this paper with a (graded) permutation operator P and a unity operator I [65]. The Gaudin-Yang model, discussed in introduction, is related to algebra symmetry gl(2|1) while the two-component Bose gas to algebra symmetry gl(3). In this paper we do not give the explicit form of monodromy matrix entries for the models under consideration, instead we refer to [10,11].
We denote Bethe Ansatz vacuum as |0 (and dual vacuum as 0|). Vacuum eigenvalues of the diagonal entries of the monodromy matrix T ii 3 are denoted by λ i : (2.23) Following [1] we use the concept of the two-site (also called partial) model for the description of subsystems. Thus we consider the system of length L and split it into two subsystems [0, x] and [x, L]; this will be useful for the computation of the two-point correlation functions. We denote the partial quantities belonging to the first (second) subsystem by superscript (1) (correspondingly (2)). For example Q (1) i denotes the number of particles of type i, i = 1, 2, in the first subsystem. Furthermore, r (k) i (t), i = 1, 3, k = 1, 2 denote the ratios of eigenvalues of the diagonal elements of partial monodromy matrix belonging to subsystem k. Obviously r i = r (1) i r (2) i , i = 1, 3. We apply the following notations for logarithmic derivatives of r i (w): We extend notations F (ū) = F (u 1 )F (u 2 ) . . . introduced for functions to operators, that commute among themselves. For instance in case of the algebra symmetry gl(2|1) related models Operators T 13 , T 31 , T 32 and T 23 do not commute among themselves. For example In these cases we introduce the operator products that are symmetric with respect to the parametersv. The transfer matrix, defined as is a generating function of (mutually commuting) integrals of motion In ABA the Bethe vectors are special polynomials in the monodromy matrix entries acting on vacuum. For the simplest case of the gl(2) related models they are given by 4 For models related to higher rank symmetries we need to apply the so-called nested Bethe Ansatz, where we have different types of Bethe roots corresponding to the different levels of the nesting procedure. In concrete physical models the different types of spectral variables correspond to the different physical degrees of freedom. In the gl(2|1) related models the states are described by two types of variables [10,67]. The explicit form of the Bethe vectors is given by 30) and the dual Bethe vectors are given by In the following we will denote #ū = a and #v = b. Above the sum is taken over partitions In the general case {ū,v} are sets of arbitrary complex numbers, and we call Bethe vectors off-shell. In the case when the sets {ū,v} satisfy the system of Bethe Ansatz equations (BAE) Bethe vectors become eigenvectors of the transfer matrix, we call them on-shell. In our notations the (twisted) Bethe equations can be written as Following [68] we use the concept of the generalized model, in which sets {r 1 (u k )}, k = 1, . . . , a, {r 3 (v j )}, j = 1, . . . , b are treated as sets of free parameters, without any reference to the particular model. For the off-shell Bethe vector instead of two sets of parameters {ū,v} we have now two additional sets of parameters The generalized model is a convenient instrument for computations in the ABA approach and will be used throughout the paper. It allows to establish expansions (2.6)-(2.8) for arbitrary integrable model related to algebra symmetry gl(3) and gl(2|1).
In sections 5.4 we specify functions {r 1 , r 3 } as r 1 = 1, r 3 (v j ) = e iv j L and produce results of section 2 for the Gaudin-Yang model from the general formulas.

Expansion theorem for the mean values
Here we give a proof LeClair-Mussardo series using standard tools of nested algebraic Bethe Ansatz.

Analytic properties of a scalar product
We define the scalar product S a,b (ū C ;v C |ū B ;v B ) of two generic (off-shell) Bethe vectors as Our expansion theorem introduced below is built on certain singularity properties of the scalar products and form factors. Thus we start with investigating the poles of the scalar product. It was proven in [69] that the scalar product can be expressed via highest coefficients Z n,m

(3.2)
Here and further #ū B = #ū C = a and #v B = #v C = b. The sum is taken over partitions The pole structure of the highest coefficient is known. We are interested in residues that appears in the limitū Let us consider now the analytic properties of the scalar product of two Bethe vectors. . .
The second case is (3.6) The simplification of the last line gives .
Taking the sum of both contributions (3.4) and (3.6) we come to where in the r.h.s. the vacuum expectation values are modified as .
We can take now the limit of first factor The first contribution to the residue comes from the partitions where v C a ∈v C I , v B a ∈v B I and the second from the partition v C a ∈v C II , v B a ∈v B II . The first contribution is (3.11) The simplification of the last line gives . (3.12) The second contribution is (3.13) The simplification of the last line gives . (3.14) Taking the sum of contributions (3.11) and (3.13) we arrive at (3.15) Here the vacuum expectation values at the r.h.s. are modified as and we can take the limit of the first factor there

Analytic properties of form factors
It is our main goal to study form factors of products of local operators. It is known that in the fundamental models local operators can be expressed directly via the entries of the monodromy matrix [70,71]. Therefore, we first study the analytic properties of them.
In models of quantum gases the solution of the quantum inverse problem is not known. However, taking into account that such models can be derived from the fundamental one using a special type of scaling limit [72,73], we assume that the properties of form factors in Gaudin-Yang model are similar to the properties of form factors in the fundamental model. This was explicitly demonstrated in models related to algebra symmetry gl(2) in [74].
We start with the analytic properties of form factors of the matrix elements of a monodromy matrix. They follow from the formulas of multiple action of operators T i,j on the Bethe vectors and from the analytic properties of the scalar product, established earlier. Formulas for the multiple action in the case of algebra symmetry gl(2|1) were derived in [75].
Consider for instance the multiple action of the operator T 31 (z) Here the sum is taken over partition {ū,w} =η → {η I ,η II , . . . }, {v,w} =ξ → {ξ I ,ξ II , . . . } with the restrictions #ξ I = #η I = #z and K n (x|ȳ) is the Izergin-Korepin determinant, it is defined for any two setsx,ȳ with #x = #ȳ = n: The first pole in (3.18) is present if u B a ∈η III . The singular part is proportional to the singular part of products ū C ;v C |η III ;ξ III . We write separately only the part under the sum over partitions that depends on u C a and u B a . The following factor appears in front of these singular parts . (3.20) After an elementary transformation, we get The second pole is present if v B b ∈ξ III . We write only the part under the sum over partition that .

(3.22)
After an elementary transformation and the replacement r 1 (η II ) by Since we have the symmetry w.r.t to permutations of pairs {u k , r 1 (u k )} ←→ {u n , r 1 (u n )} for arbitrary k, n = 1, . . . , a and {v k , r 3 (v k )} ←→ {v n , r 3 (v n )} for arbitrary k, n = 1, . . . , b the same property is true for each rapidity in the setū andv.
Using the action rules for all T ij (w), i, j = 1, . . . , 3 it is easy to obtain that the same analytical properties will hold for arbitrary T ij (w). It is also easy to check in the same way that the same analytic properties holds for superposition of operators T ij (w I )T kl (w II ) . . . . Thus, the form factor of an arbitrary operator O will satisfy the same properties: The only difference in the case of algebra symmetry gl(3) will be replacement of the factor

Expansion formula
The proof of the expansion formula is loosely based on the proof for the algebra symmetry gl(2) case given in [43,64]. Consider the diagonal form factor Here it is understood that the diagonal limit is performed on the off-shell form factor. It follows from the above that the diagonal form factor is a linear function of the variables Z j , Y k with the coefficients cf (ū C a , u C a )f (ū B a , u C a ) and (3.27) Now we introduce the symmetric form factor. As given in the previous Section, its definition in the general case is Here we regard {ū C ,v C } and {ū B ,v B } as solutions of the Bethe equations, i.e. limitū C →ū B , v C →v B is taken after the r 1 (u k ), r 3 (v k ) are expressed via r.h.s. of (2.32) (without twists).
In the case of models related to symmetry gl(2) it was proven already in [64] that the symmetric form factor is the mean value in the case when the parametersZ are chosen to be zero. In the nested case a similar property holds: (3.29) 6 Factors f (w, u C a ) and g(v C b ,w) are absorbed here into the normalization of the form factor (see (3.26)).

31)
and G is the Gaudin determinant. The cardinalities of the sets are #ū I = m, #v I = n and #ū II = a − m, #v II = b − n.
Proof. Formula (3.30) can be proven by induction. Let us start from the first nontrivial case #ū II ≡ a 2 = a, #v II ≡ b 2 = b for which there is one non-zero term in the r.h.s. of (3.30) and this is F a,b . This term isZ-andȲ -independent. The l.h.s. is, of course,Z-,Ȳ -independent too, it is clear from (3.27).
The matrix elements of the Gaudin matrix are given later, see (4.3)-(4.4). It is clear from that explicit form that the derivatives of S a,b w.
here Z k , Y j are modified as in first case and as in the second case. Supposing that (3.30) holds for some arbitrary a − 1 we can prove that it holds for a also. Take the derivative of (3.30) w.r.t. Z j . Using (3.32) in the r.h.s. of (3.30) we immediately reproduce the same expansion but with setū replaced byū a that is, according to the induction assumption, F O (ū a ;v|Z a ;Ȳ ). At the l.h.s. we use (3.27) and obtain exactly the same. Therefore theZ-dependent part coincides. The same procedure can be done forȲ -dependent part. In order to compareZ-andȲ -independent part it is enough to send allZ andȲ to zero. In this case only one contribution F a,b (ū;v) survives and it is equal to F O a,b (ū;v|{0}; {0}) in the l.h.s. by definition of F a,b (ū;v).
After dividing (3.30) by norm of eigenvectors ū;v|ū;v (see [76]) we reproduce formula (2.6) Here O in the l.h.s. denotes the normalized expectation value and F m,n denotes normalized F m,n .

(3.36)
It is easy to prove, that formula (3.35) is still valid in the case of algebra symmetry gl(3) up to replacement of the Gaudin matrix G(ū;v) and F a,b by the corresponding analogs.

Correlation functions via the connected form factors
The original LeClair-Mussardo series [35] was built on the so-called connected form factors, which were defined as the pole free part of the diagonal limit of the on shell form factors. The relations between the symmetric and connected form factors were later clarified in [62]. In this work we present the LM series using the connected form factors also in the nested Bethe Ansatz. The connected form factors can be defined similarly to [35,62] through the limit where the left Bethe vectors is defined as and finite part Fin{. . . } is defined in such a way that all singular terms including ratios of form ε i /ε j , ε ′ i /ε ′ j for i = j and ε ′ i /ε j , ε i /ε ′ j are discarded. Combining the arguments of [43,62] and our computations with the scalar products and form factors in the nested case it is possible to prove the following expansion. . (3.39) HereḠ (ū I ,v I |ū;v) is defined as the matrix built from the columns and rows of the Gaudin matrix G(ū;v) that belong to the subsetsū I ,v I , and we stress explicitly that in contrary to G(ū I ;v I ) this matrix depends on the full sets of parametersū,v. The normalization of F c here coincides with the normalization of (3.36).

Thermodynamic limit
In this Section we take the thermodynamic limit of the expansion theorem, leading to the LM series in nested Bethe Ansatz systems.

Distribution of Bethe roots
The distribution of Bethe parameters in the thermodynamic limit is given by integral equations for densities, which follow from the Bethe equations [61]. We restrict ourselves to cases where, at least in the ground state, solutions of the BAE are real; the Gaudin-Yang model with c > 0 belongs to this class of models. In such a case the integral equations for the ground state are where ρ 1 , ρ 2 are densities of particles, p ′ 1 (x) = −i∂ log r 1 (x), p ′ 2 (x) = −i∂ log r 3 (x) and B (resp. Q are Fermi boundaries for the first (resp. second) level of Bethe rapidities.
In the Gaudin-Yang model total densities n ↑/↓ depend on Q and B, and in principle this relation can be inverted to find Q and B in terms of the physical quantities n ↑/↓ . In practice this involves the solution of the linear integral equations above. The relations simplify in the large imbalance limit, as discussed in detail in Section 5.6.

Ratio of Gaudin determinants: Fermi case
Here we calculate the ratios of the Gaudin determinants det G(ū I ,v I ) and det G(ū,v) in the thermodynamic limit. We use the approach developed earlier [1,3,43]. The Gaudin determinant for the algebra symmetry gl(2|1) related models was calculated in [76].
G is an (a + b) × (a + b) block-matrix. The diagonal blocks are The anti-diagonal blocks are After extracting these diagonal parts from the Gaudin determinant, it can be presented as a product det G = det θ det G, (4.6) and det θ is a determinant of a diagonal matrix with components {θ a+1 , . . . , θ a+b }. Consider first the caseū I =ū j ,v I =v (orū I =ū,v I =v j ). Denote as det θ ′ j the determinant extracted from det G(ū j ,v) (or from det G(ū,v j )). The ratios of det θ ′ j /θ for the first and the second case can be presented as (4.8) In the thermodynamic limit and (4.8) correspondingly leads to where we defined the weight functions (4.11) Note that these weights contain an explicit factor of c −1 ; this factor comes from our choice of normalization of the Bethe states and thus the form factors. This factor is not included in the ω-function found in the works [1][2][3]7]; instead it is carried by their form factors. In the thermodynamic limit the remaining part of (4.6), detĜ(ū,v) leads to a Fredholm determinant with the kernel K, similarly to the algebra symmetry gl(2) related models [3]: (4.13) Since detĜ(ū j ,v) and detĜ(ū,v j ) become the same Fredholm determinant in the thermodynamic limit, detĜ(ū j ,v)/ detĜ(ū,v) = detĜ(ū,v j )/Ĝ(ū,v) = 1 + O(L −1 ). (4.14) In the caseū I =ū j 1 ,...,j k (orv I =v j 1 ,...,j k ) the same recipe can be applied with minor changes. Each u j 1 , . . . , u j k (and each v j 1 , . . . , v j k ) will produce a corresponding density function in the denominators of (4.10) and the corresponding weight in the exponent. The approach remains valid as long as we restrain ourselves to #v I ∼ #v and #ū I ∼ #ū. In practice we always truncate the series (3.35) to small values of #ū II , #v II , so this restriction is always respected.

Thermodynamic limit of the expansion formula
Now we are able to take thermodynamic limit of (3.35). Repeating the steps from [43] and taking into account (4.6), (4.10), (4.14) we arrive at the following expression where #t = m and #s = n. This is one of our main results; it can be considered as a generalization of the LeClair-Mussardo formula to nested Bethe Ansatz systems. The computation of the symmetric form factors entering the series is presented in Section 5. It can be shown that formally the same series also holds in the case of algebra symmetry gl(3). The only difference is the replacement of the functions ω (2) (v) by under the integral and in the explicit form of the coefficients F m,n .

Thermodynamic limit of the connected form factor series
It can be shown absolutely similarly to a case of algebra symmetry gl(2) related models [43] that the thermodynamic limit of the expansion theorem using the connected form factors is given by

Correlation functions
The techniques developed above allow to expand correlation functions in a series of form factors (symmetric or connected). The representation is given by an infinite integral series, and it is a natural question whether it is possible to obtain meaningful expressions for correlators via this series.
In the case of the Lieb-Linger gas it was observed that the series similar to (4.15) has very good convergence in the strong coupling limit [1,3,43]. This was also justified from the field theory limit and checked numerically in [77,78] for the case of correlator (ψ † ) k (ψ) k , k = 2, 3. Furthermore, the expansion of correlation functions using relatively similar to our approach studied by [2,3] also possess good convergence.
We expect that the similar properties will hold in the two-component models. We show below, that in the spin-1/2 fermion model the convergence is guaranteed if the two dimensionless parameters Q/c and B/c are small. Smallness of the two combinations corresponds to the strong coupling and strong imbalance between the components with different spin projections.

Generating functions
As it was pointed in Section 3.2 we expect that the LM-series can be applied also in the nonfundamental models, despite that it is often not easy to express physical operators via the monodromy matrix entries in the general case. Often it can be done by proper scaling limit from the fundamental model; here, however, we use an approach that allows to compute the two-point correlation functions in a more direct way.
We denote by L the system size, the segment [0, x] will be referred as the first subsystem while [x, L] will be the second. Consider the generating function for the form factors of diagonal operators exp αQ (1) , here αQ (1) k , k = 1, 2 are operators that measure the numbers of particles of type k in the subsystem 1 and α 1 , α 2 we call twist parameters. The generating function exp (αQ 1 ) contains enough information so that it yields both the onepoint and two-point correlators. It is proven in Appendix A, that using the concept of the partial model (see description after (2.22)), similarly to the case of the algebra symmetry gl(2) related models, the (un-normalized) form factor of O = exp αQ (1) can be presented as where it is understood that {ū B ,v B } and {ū C ,v C } satisfy BAE (2.32) with a trivial twist κ 1 = κ 2 = κ 3 = 1. In the formula above Θ α m,n is a quantity that will be called the highest coefficient following [1,3] 7 . It is proven in Appendix A that it coincides with a scalar product of the on-shell and the twisted-on-shell Bethe vectors with a twist α = {α 1 , 0, α 2 }: Here ℓ i = r (1) i , i = 1, 3 (partial r i , that describe the first subsystem, see the description above (2.24)).
Correlation functions of densities can be expressed via the generating function using the generation functions. Thus in the case of Gaudin-Yang model 2 is a particle number in the second subsystem [x, L], and in the last formula we replace α 2 → −α 2 that just changes the common sign. Notation α = 0 here and further means that α 1 = α 2 = 0.
In a similar way correlators in other models can be expressed via derivatives of exp(αQ (1) ) (for example correlators of electrons densities in super-symmetric t-J model, densities in the Fermi-Bose mixtures). For the lattice models the derivatives w.r.t. x would be naturally replaced by the finite differences.
Taking the second derivative of (5.1) w.r.t. α s , s = 1, 2 at α = 0 and taking into account, that Θ α=0 m,n (ū B ;ū C |v B ;v C ) = δ m,0 δ n,0 due to the orthogonality of the on-shell Bethe vectors, we arrive at here the derivatives are taken at α = 0, we omit superscript α in this case, and the operator in question is O = exp(αQ (1) ). The terms with #ū C II = #ū B II = a 2 = a, #v C I = #v B I = b 1 = b and #ū C I = #ū B I = a 1 = a, #v C II = #v B II = b 2 = b are written separately and 0 ≤ a 1 ≤ a, 0 ≤ b 1 ≤ b with the restriction that the summation over {a 1 , b 1 } excludes the cases {0, 0} and {a, b}.
As an obvious consequence of the proposition 5.1 we can rewrite (5.6) as (5.9)

Computation of the symmetric form factors
In order to calculate correlators using (2.8) the explicit expressions for the symmetric form factors F m,n is required. We do not have an efficient algorithm for the calculation of F m,n andÎ a,b; k,l , however particular cases of small particle numbers can be considered. We can apply the formula (5.9), and the remaining task is to substitute the representation for highest coefficient of correlation function there. These explicit representations of Θ α m,n were derived in a compact determinant form in [76]: Here the block matrix N is defined as , j = 1, . . . , m, k = 1, . . . , m, , j = 1, . . . , m, k = 1, . . . , n, , j = 1, . . . , n, k = 1, . . . , m, , j = 1, . . . , n, k = 1, . . . , n.

(5.11)
Here we denote κ i = e α i , where α is defined in (5.9). Thus, it remains to calculate the sums over partitions in (5.9). In general this is a complicated problem even for the algebra symmetry gl(2) based models. However, for small a, b it is easy to obtain the results, since the sums contain only few terms in these cases. Further in this section we restrict ourselves to the Gaudin-Yang model. We consider here the case of repulsion, i.e. coupling constant c > 0. It was shown by [61] that there are no string solutions of the Bethe equations in this case. The case of attractive interaction will be considered in a separate work since it requires a consideration of the string solutions. Notice, that in comparison with (2.2) our Bethe equations (2.32), (4.1) differ by shift v → v + ic/2. Thus, further in this Section we will assume the shift of v by −ic/2 every time we discuss the Gaudin-Yang model 8 . We remind that in the Gaudin-Yang model r 3 (v) = e iLv , r 1 (u) = 1 and ℓ 3 (v) = e ixv , ℓ 1 (u) = 1.

Symmetric form factor
Consider the symmetric form factors for the correlators of densities. We take the diagonal limit v C → v B , u C → u B in (5.9) as defined in (2.7).
In models where ℓ 1 = 1, such as the Gaudin-Yang model, we are not able to distinguish different "modes" in Fourier series. Thus we make the expansion w.r.t. ℓ 3 . Then A n a,b has only one superscript n = #v, and the expansion is defined as 5.4 Symmetric form factor F n,m for q(x)q(0) We denote Fourier coefficient in front of oscillating parts of correlation function as I a,b . One and two particles oscillating parts are given by I 0,0 = 0, I 1,0 = 0, I 0,1 = 0, (5.14) where and Here A i k,l are the Fourier coefficients introduced in (5.13). Three particle oscillating parts are given by I 3,0 (ū) = 0, I 0,3 (v) = 0, I 2,1 (ū;v) = 0, (5.18) where (5.20) I 3,1 , I 2,2 , etc. can be calculated too, but the explicit formulas for them become quite lengthy. One particle F are F 1,0 = 0, F 0,1 = 0. (5.21) Two and three particle form factors are where the second (non-oscillating) term is (5.25) and where the second (non-oscillating) term is We can now argue that the series (4.15) converges. We follow again the arguments of [1,3]. We observe, that at large c F a,b ω (1) ω (2) are suppressed by the power c 2−a−b . Using explicitly expressions for F m,n it is easy to estimate now that n, m particle contributions behave as where Taking into account terms proportional to J ≡ ℓ 1 (ū C )/ℓ 1 (ū B ) − 1 we obtain the symmetric form factors.
where the second (non-oscillating) term in (5.36) is We observe, that F a,b ω (1) ω (2) at large c is suppressed by the power c 2−a−b in the case of the mixed densities correlator too.

The limit of large imbalance
Here we show that if we fix Q/c, then the limit of small B/c corresponds to large imbalance.
In the repulsive Gaudin-Yang model the ground state consists of only real roots and the densities of Bethe parameters satisfy the system of integral equations (4.1) with p ′ 1 (x) = 0 and p ′ 2 (x) = 1. The densities of fermions are given by The integral equations (4.1) can be solved iteratively by assuming a small B/c parameter. In the first approximation the integrals over ρ 1 (λ) can be approximated simply by ρ 1 (0). This gives Combining the two equations we get: Regarding the total density we have The last two equations determine n ↑/↓ including up to first order in B.
Note that these equations are exact in Q/c. It is straightforward to perform a further expansion in Q/s, if it is needed for some purpose.
Let us also compute the local 2-body correlations in this expansion. It follows from the form of the Hamiltonian and the Hellmann-Feynman theorem that where E 0 is the energy of the ground state at fixed particle densities. For the energy density the leading terms in the expansion in B/c are known, see for example formula (2.123) of [79] which includes even the O(B 2 /c 2 ) corrections. The linear term can be computed using the expansion of the densities above. In our conventions it reads Taking a derivative with respect to c, while keeping n and n ↓ fixed we obtain Combining this with the expansion (5.43)-(5.44) we see that

Correlation functions
Now using (4.15), (5.4), (5.5) and F, calculated in the section 5.2, we are able to calculate the correlation functions of the highly polarized gas in a strong coupling regime explicitly. For the first two terms we get Nontrivial symmetric form factors are F 0,2 (∅;v)ω (2) (v) = −2 e iv 12 x + e iv 21 x + 4 e −8B/(πc) + O(1/c), After substitution to (5.49) we get immediately In the same way substituting non trivial symmetric form factors of q ↑ (x)q ↓ (0)

Conclusion
In this work we computed the LeClair-Mussardo series for the correlation functions in models related to algebra symmetries gl(2|1) and gl (3). We presented two versions involving the symmetric and connected form factors. It is important that even though we focused on the particular example of the spin-1/2 Fermions, the expansion theorems are model independent and could be applied in other situations as well.
An important finding is that the symmetric and connected form factors that enter the LM series do not satisfy the same selection rules as the on-shell Bethe states i.e. form factors can be non-zero even when the corresponding on-shell Bethe states do not exist. In our computations this is just a result of the algebraic procedure with which we obtained them. We believe that the physical reason for this is that such form factors describe multi-particle processes in the presence of the other particles, therefore it is not required that the same selection rules should hold as in a situation where the same particles are placed in the vacuum [36,38,43,62].
In our concrete example we only treated the ground state of the Gaudin-Yang model, which consists only of real roots. In a future work we also plan to treat cases with string solutions. This allows to proceed to computations of finite temperature and dynamical correlation functions. Another focus is to study the two-component Bose gas, related to the algebra symmetry gl (3).
Furthermore, we plan to extend this method to other models. At present there are very few results for correlation functions in spin chains with higher rank symmetries, especially in models with o(N) symmetry (see for example [50]). Our method could be helpful in finding new results. Also, it would be interesting to consider the sine-Gordon model, whose finite volume mean values were investigated in [80,81]. As far as we know, it is not yet clear whether a finite volume expansion theorem exists in that model. Our results could be helpful in settling this issue.

Acknowledgments
A.H. is grateful to Frank Göhmann and Andrii Liashyk for numerous discussions. The authors are also thankful to Gábor Takács andÁrpád Hegedűs for useful discussions. A.H. was partially supported by the BME Nanotechnology and Materials Science TKP2020 IE grant of NKFIH Hungary (BME IE-NAT TKP2020). The authors were also supported by the grant K-16 no. 119204 of NKFIH. B.P. acknowledges further support from the the János Bolyai Research Scholarship of the Hungarian Academy of Sciences, and from theÚNKP-19-4 New National Excellence Program of the Ministry for Innovation and Technology. Most of this research work was performed while B.P. was employed by the "MTA-BME Quantum Dynamics and Correlations Research Group", at the Budapest University of Technology and Economics.

A Appendix
Here we give a proof of (5.1). We follow the method used in [1,3,82].
In models based on the algebra symmetry gl(2|1) the off-shell Bethe vectors satisfy the following co-product property here the superscripts (i), i=1,2 label explicitly that these quantities belong to the corresponding subsystems. The similar property for a dual Bethe vector is The operator Q i , i = 1, 2 counts the number of particles of type i in subsystem 1. Operator exp αQ (1) acting on Bethe vectors of first subsystem produces the factor κ a 1 1 κ b 1 2 = exp(α 1 a 1 + α 2 b 1 ), where a 1 , b 1 are numbers of particles in the first subsystem, i.e. the eigenvalue of Q (1) i . Let us consider simplest example i = 1, the general case exp α 1 Q can be proven in the same way. From the co-product property of Bethe vectors in the algebra symmetry gl(2|1) based models we can write the following representation for the (un-normalized) matrix elements of the exponent exp αQ where S ℓ , ℓ = 1, 2 are scalar products of the partial Bethe vectors (A.1). Substituting now the scalar products S 1 , S 2 in form (3.2) into (A.3) we get 9 and we express all r 3 (v B i ) and total r 3 (v B i ), and substitute Bethe equations for r 3 . (A.7) Further we collect the rest of the factors . (A.8) In a similar way simplifications for sets {ū B ,v B } can be performed. Finally we make some simplification of pre-factor with r 1 (ū C 23 )r where h α (v) are the one-particle charge eigenvalues. Note that here only the second level rapidities v contribute; they are also called the momentum-carrying particles. The formula (B.2) also holds for the total particle density operator q = q ↑ + q ↓ , in which case h(v) = 1. In the case of the density q ↓ we have formally Comparing (B.2) to (3.35) we can compute the symmetric form factors recursively. These results can then be compared to direct computations of the form factors, eventually confirming (3.35). We consider the simplest cases of this procedure.

B.1 Symmetric series
Consider the algebra symmetry gl(2|1) related model with ∂ v log r 3 (v) = ip ′ L, and r 1 = 1. Let us compute explicitly F of conserved charges for the small values a, b.
For case b = 1, a = 1 For case b = 2, a = 0 we have similarly And for case b = 2, a = 1 Substituting in the last expression ρ(u 1 ,v) = p ′ (v 1 )p ′ (v 2 )L 2 (t(v 1 , u 1 ) + t(v 2 , u 1 )) + t(v, u 1 )L (p ′ (v 1 ) + p ′ (v 2 )) , (B.12) we obtain Continuing this line of argument every symmetric form factor can be expressed using the charge eigenvalues, the kernel t(u, v) and some graph theoretical tools. This leads to a generalization of the graph theoretical results of [64] to nested Bethe Ansatz systems. However, further details of this computation are not important for our present purposes.
Nevertheless we remark that it is also possible to find the symmetric form factors of the current operators, which describe the flow of the charges. This too can be done using the methods of [64]. The resulting form factors could be used to give an alternative proof of the current mean values in nested systems, first found in [84]. However, this is beyond the scope of the present paper, and we do not pursue it.
Instead, let us construct a check of the above formulas, in the the case of the density operators q ↓ (x). Taking the first derivative of (5.1) w.r.t. α 1 , and taking into account that Θ α=0 a,b = δ 0,a δ 0,b we obtain the off-diagonal form factors (B.14) where notation exp αQ (1) is the same as in 5.1, and Θ α a,b is given explicitly by (5.10). Applying the symmetric limit (see 5.3) we can reproduce (B.6), (B.8), (B.13). This confirms the expansion theorem (2.6) for the specific charge operators.

B.2 Connected series
Applying (B.2) to (3.39) we can compute connected form factor using the recurrence procedure similar to used in (B.1). Consider again model related to algebra symmetry gl(2|1) with ∂ v log r 3 (v) = ip ′ L, and r 1 = 1.
In a case b = 2, a = 0 we have taking into account that ρ(∅;v) = L 2 h α (v 1 )h α (v 2 ) and using results of (B.15)-(B.17) it is easy to obtain that F c (∅;v) = 0. In case b = 2, a = 1 we have and using (B.12) andρ(u 1 ; v 1 |v 2 ) = p ′ (v 1 )L[t(v 1 , u 1 ) + t(v 2 , u 1 )] + t(v 1 , u 1 )t(v 2 , u 2 ) we arrive at Continuing a similar procedure we can calculate further F c and even establish a general graph rule, similarly to case of symmetric form factors [64]. Again using (B.14) and applying to this formula corresponding symmetric limit (3.37) we can easily reproduce results that follow from (B.15)-(B.20). This confirms the expansion (3.39) for the charge operator.

C Appendix
Correlation function of densities in the Lieb-Liniger model was computed in [1,3].
q(x)q(0) = Q 2 π 2 + Here we reproduce this result using LM-theorem. The Lieb-Liniger model is related to the algebra symmetry gl(2). Using (4.15) and (5.4) we arrive at q(x)q(0) = 1 2 where I k are irreducible parts for algebra symmetry gl(2) related models and ω is a weight function.