Low-density limit of dynamical correlations in the Lieb-Liniger model

We derive explicit expressions for dynamical correlations of the field and density operators in the Lieb-Liniger model, within an arbitrary eigenstate with a small particle density ${\cal D}$. They are valid for all space and time and any interaction strength $c>0$, and are the leading order of an expansion in ${\cal D}$. This expansion is obtained by writing the correlation functions as sums over form factors when formally decomposed into partial fractions.


Introduction
The Lieb-Liniger model is a key paradigm of many-particle systems [1][2][3]. In the repulsive regime, it is considered as one of the simplest interacting quantum integrable model, for having the simplifying feature of involving only real rapidities [4,5]. The main objects of interest are the correlation functions of local observables in the thermodynamic limit, that are the macroscopic output of the model resulting from the short-range interactions between the bosons. Moreover their Fourier transform is directly measurable in cold-atoms experiments [6][7][8]. Although exactly solvable, the computation of correlation functions in quantum integrable models is a notoriously difficult problem. In particular the computation of dynamical correlations for an arbitrary interaction strength c and within arbitrary eigenstate is an open problem.
Special cases of this problem were studied and solved in the past. The first of these special cases is the impenetrable bosons limit c → ∞, that can be reformulated in terms of free fermions [9]. Here, both because of the free fermionic nature of the model and the simple structure of the form factor 1 , the Lehmann representation of dynamical correlations can be fully resummed into a Fredholm determinant [11][12][13][14] and its asymptotics extracted with differential equations [15,16]. Another important special case is the ground state static correlation at finite c, that was treated within the Algebraic Bethe ansatz framework [17] and led to the first ab initio calculation of the critical exponents previously predicted by Conformal Field Theory and Luttinger Liquid theory [33][34][35][36]. This approach was then generalized to static correlations within arbitrary eigenstates [31,32]. The particular case of static correlators within thermal eigenstates was also studied with quantum transfer matrix methods [28][29][30]. The full asymptotics of ground state dynamical correlations were derived in [18][19][20] from form factor expansions, confirming predictions of Non-Linear Luttinger Liquid theory [37][38][39][40][41][42].
Progress on the general case of dynamical correlations within arbitrary states at finite c has been much more limited. This calculation is very different from the special cases of zero-temperature or static cases and poses important problems. The successful methods for static correlations such as the algebraic Bethe ansatz approach or the quantum transfer matrix methods do not apply to the dynamical case, and the form factor expansion used to compute ground state dynamical correlations relied on combinatorial identities that are usable for zero-entropy states only. This general case has however been studied through several approaches. Firstly, Generalized HydroDynamics (GHD) provide predictions for the leading asymptotics of the dynamical correlations of conserved charges such as the density within arbitrary macrostates [45][46][47]. However the approach cannot a priori be applied to compute the next corrections, restricting the dynamical structure factor to small frequency and mo-mentum, and does not apply neither to semi-local operators such as the field correlations. Secondly, numerical summations of dominant form factors proved very efficient and permitted to obtain numerical estimations of the dynamical structure factor on the full plane [43,44]. On the Bethe-ansatz calculations side, one-point functions within arbitrary eigenstates were studied in [48][49][50]. An approach based on thermodynamic form factors was started [51][52][53][54][55][56], but still involves non-integrable singularities and so requires a particular understanding of this feature. A regularized form factor expansion was derived for the XXZ spin chain in [57]. In [58] was computed the full spectral sum at order c −2 for all root densities and all momentum and frequency, involving one-and two-particle-hole excitations. It showed the necessity for a fine-tuned treatment of the non-integrable singularities that is crucial for e.g. detailed balance to be satisfied in thermal states.
The objective of this paper is to derive the full dynamical correlations for all space and time and all interaction strength c, in the limit where the particle density D of the averaging state becomes small. This low-density limit is defined in terms of a partial fraction decomposition of the form factor, initially introduced in [10] in a model that can be reformulated into free fermions. This decomposition naturally organizes the spectral sum as an expansion in the particle density D of the averaging state, and the low-density limit is defined as the leading term of this expansion.
This result provides another limiting case where the initial problem becomes solvable, namely the computation of dynamical correlations for all space and time and arbitrary c within finite-entropy macrostates. Moreover the framework also enables to compute the subleading corrections in the particle density, as was explicitly shown in [10] for the Transverse Field Ising Model. The computation of these subleading corrections in the interacting case however comes with higher technical difficulties and should be the object of further work. Finally, this low-density limit calculation sheds light on the structure of the spectral sum and the nature of the states contributing in the thermodynamic limit.
In Section 2 we introduce the Lieb-Liniger model and recall known results on its form factors. In Section 3 we define what is meant by low-density limit of the dynamical correlations, taking the field correlations as an example. In Section 4 we compute the low-density limit of the field two-point function (70), and in Section 5 the low-density limit of the density two-point function (96).

Definition
The Lieb-Liniger model [1] is a non-relativistic quantum field theory model with Hamiltonian where the canonical Bose field ψ(x) satisfies equal-time commutation relations We will impose periodic boundary conditions. For later convenience we define the time-t evolved version of the field ψ(x, t) = e iHt ψ(x)e −iHt . We also define the density operator at position x and its time-t evolved version σ(x, t) = e iHt σ(x)e −iHt .

The spectrum
The Lieb-Liniger model is solvable by the Bethe ansatz: an eigenstate |λ λ λ with N bosons can be written as with the B(λ)'s some creation operators, |0 the pseudo-vacuum and the λ i 's some rapidities that satisfy the following set of 'Bethe equations' The energy E and the momentum P of this state read It is convenient to express the Bethe equations in logarithmic form with I k an integer if N is odd, a half-integer if N is even. For c > 0 all the solutions to this equation are real [3]. We will denote the particle density of the eigenstate |λ λ λ .
Our strategy is to use a Lehman representation in terms of energy eigenstates |µ µ µ = |µ 1 , ..., µ N ′ , where {µ 1 , . . . , µ N ′ } are solutions to the Bethe equations (7) to rewrite the correlation functions as sums of form factors over the full spectrum. For the two-point function of a generic operator O this representation reads The (normalized) form factors of the field and density operators between two Bethe states |λ λ λ , |µ µ µ with respective numbers of Bethe roots N, N ′ have been computed previously [21][22][23][24][25][26].
In the case of the field operator, it reads for any p, s = 1, ..., N . The various terms entering this expression are and the N × N matrix , (14) and finally N λ λ λ = det G(λ λ λ) , with the Gaudin matrix [27] The form factor of the density operator reads for any p = 1, ..., N , with and with now

Root densities
In the thermodynamic limit, any sum of a non-singular function over the Bethe roots can be expressed in terms of a root density that characterizes a macrostate as far as such quantities are concerned However if the function is singular the result will in general depend on the representative state of the macrostate, see [58]. It is customary to introduce the hole density ρ h (λ) defined by 3 Definition of the low-density limit The purpose of this section is to define what is meant by the low density limit of correlation functions. It is defined as the leading order of an expansion in D, obtained by decomposing the form factor in partial fractions. As a consequence it is an expression valid for all x, t and c, that becomes closer to the dynamical correlations as the particle density D of the averaging state becomes smaller.
This definition requires some technicalities, but is rigorous and allows for a computation of the next orders, as shown in [10] for a model that can be reformulated into free fermions. However, it a priori lacks some intuitive picture. For that reason we provide an interpretation of this low-density limit so defined as a Lehmann representation in terms of the low density limit of the form factor. The reasoning is here rather different, and consists in first approximating the form factor by the thermodynamic limit value they take when one of the two states is a dilute state, i.e. such that for any pair i, j we have L(λ i − λ j ) → ∞. In this limit, the spectral sum of the dynamical correlations indeed matches the leading order of the expansion in D, providing an interesting and intuitive consistency check. But it must be stressed that the right definition of the low-density limit is more general than this intuitive calculation, since it only requires D to be small, not the root density ρ(λ) to be small everywhere.

Recall
We recall that the partial fraction decomposition (PFD) of a ratio of two polynomials with P 0 (X) a polynomial, and B i,ν coefficients given by The polynomial P 0 (X) can be determined by studying e.g. the large X behaviour of the ratio of the two polynomials on the left-hand side of (22).

The poles of the normalized form factor
We consider λ λ λ and µ µ µ two sets of respectively N and N − 1 rapidities, and would like to apply a partial fraction decomposition to the square of the normalized field form factor, with respect to each of the µ i 's successively. The first task is to identify the poles of a µ i at fixed other µ j 's. There are a priori three types of poles for µ i 1. Double poles in λ j for all j 2. Simple poles in µ j ± ic for all j 3. Poles corresponding to zeros of the determinant of the Gaudin matrix N µ µ µ We remark that the last two types of poles come from the fact that we consider normalized form factors. We also remark that since some of the entries of the Gaudin matrix diverge when µ i − µ j → ±ic, the second type of pole could happen to be absent from the full normalized form factor, but this possibility will not be relevant to our discussion. Lastly we notice that the last two types of poles are in fact never attained when all the roots µ i are real, which is always the case if c > 0. Indeed, when the roots are real, the Gaudin matrix is strictly However, when performing the PFD of the form factor it has to be considered as a mere fraction of polynomials in µ i 's that do not necessarily satisfy the Bethe equations, and these poles have to be taken into account indeed. Among these three types of poles, the zeros of N µ µ µ are the most problematic since their location is a complicated function of the other µ j 's. For this reason we are going to consider the PFD of the form factor without these factors. Namely we define F ψ (λ λ λ, µ µ µ) by and consider the PFD of F ψ (λ λ λ, µ µ µ) only, with respect to the µ i 's.

The shape of the PFD of the reduced form factor
The reduced form factor F ψ (λ λ λ, µ µ µ), seen as a function of µ 1 , is a ratio of two polynomials with double poles in each of the λ i 's and simple poles in µ j ±ic, so one can apply the decomposition written in Section 3.1.1. Since the reduced form factor goes to zero when µ 1 → ∞, we have P 0 (X) = 0 and so one can write with where B i,ν (µ 2 , ..., µ N −1 ), C ± i (µ 2 , ..., µ N −1 ) are 'coefficients' independent of µ 1 , but that still possess a dependence in the remaining µ k 's. They have the same pole structure, except that each µ i for i = 1 does not anymore has a pole in µ 1 ± ic. The function h µ 1 (µ 2 , ..., µ N −1 ) is a function of µ 1 that has no poles in µ 1 when all the µ j 's are real. We now apply the same procedure to B i,ν (µ 2 , ..., µ N −1 ) and h µ 1 (µ 2 , ..., µ N −1 ) with respect to µ 2 to obtain the writing where a function of µ 1 , µ 2 without singularities in real µ 1 , µ 2 . We note that the poles in µ 2 arising from the function h µ 1 (µ 2 , ..., µ N −1 ) that has no poles in real µ 1 are counted through the case ν i = 0. Proceeding recursively, one obtains the writing where each ν i takes the value 0, 1 or 2, and where f are functions defined on the points The coefficients A(λ λ λ, µ µ µ, {ν}, f ) crucially do not depend on any µ i whenever ν i > 0, and are bounded regular functions of real µ i when ν i = 0.

Computing the coefficients
In the special case where all ν i > 0, the coefficients A(λ λ λ, µ µ µ, {ν}, f ) = A(λ λ λ, {ν}, f ) do not depend on any µ i 's and can be computed according to the formula If now there is a subset K ⊂ {1, ..., N − 1} such that ν i = 0 for i ∈ K, one first defines the following function of the µ i 's for i ∈ K The functionĀ({µ i } i∈K |λ λ λ, {ν}, f ) still has poles in µ i for i ∈ K since it also contains all the cases when ν i > 0. To compute the coefficient A(λ λ λ, µ µ µ, {ν}, f ) one has to remove from this function any pole in real µ i for i ∈ K. Although these formulas are not very explicit, they are still of practical use to compute the simplest terms in the PFD, as we will see below.
The functions f over which we sum in (28) are actually rather constrained. First, since the form factor vanishes whenever two µ i 's or two λ i 's coincide, we see from (30) and (31) that if ν i = 2 or ν j = 2 and f (i) = f (j), then the corresponding coefficient vanishes (namely, If we have ν i = ν j = 2 it directly follows from the absence of derivative in (30); if ν i or ν j is equal to 1, then it follows from the fact that there is a zero of order 2 in the numerator with only one derivative. Similarly, if k indices have a ν i = 1 and take the same value through f , then there is a zero of order k(k − 1) in the numerator, with only k derivatives. Hence this imposes k = 2. Thus one can impose in (28) the two following constrains: (i) that f (i) = f (j) whenever ν i = 2 or ν j = 2, and (ii) that f can take at most twice the same value.
In the following we will denote the number of elements of a set E by according to the most readable choice in the context.

A density expansion
Let us now rewrite the Lehmann representation (11) in the following way. Instead of summing over the Bethe roots µ i , we sum over their Bethe numbers J i , and trade the ordering of the Bethe numbers for a non-ordered sum with a 1 (N −1)! factor. Whenever two Bethe numbers coincide, the form factor is zero so that the two representations are indeed equivalent. Using (28), we obtain The sum over J 1 , ..., J N −1 is invariant under any changef = f •(i j) with (i j) the permutation of indices i, j, whenever ν i = ν j > 0 and f (i) and f (j) are attained the same number of times by f . Hence this sum only depends on the set of points attained a given number of times by f , not on the particular realization of the function f . To rewrite the sum without these functions f , let us define I k for k = 0, 1, 2 the set of points i in {1, ..., N } that are attained k times by f from points where ν j = 1, namely As a consequence the points in {1, ..., N } attained by f from points where . These subsets I 0 , I 1 , I 2 ⊂ {1, ..., N } have to be disjoint and to satisfy |I 0 | = |I 2 | + 1 + p with p = |{i|ν i = 0}| the number of points with ν i = 0, because f can take at most twice the same value. We will denote and parametrize When rewriting (33) in terms of these subsets, one picks a combinatorial factor corresponding to the number of functions f with such an output. Choosing the set of points where ν i = 0 yields a factor N −1 attained only once by f and where ν i = 1 a factor m! 2n+m m . Finally those attained twice by f yield a factor (2n)!!n!. Writing (2n)!! = (2n)! n!2 n yields a total combinatorial factor (N − 1)! 2 n p! .
We conclude that we can write with S n,m,p = 1 2 n p!L 2N −1 (39) The specific ordering of the µ i 's in this expression is irrelevant.
Here we have A(I 0 , I 1 , Since each choice of an index in {1, ..., N } comes with a factor D, the term S n,m,p is of order O(D 1+2n+m+p ). Hence expression (38) is an expansion in the particle density D of the averaging state.

Definition of the low-density limit of the correlation function
The low density limit of the dynamical correlations is defined as retaining only the leading term S 0,0,0 in (38). It is obtained with p = 0 and I 1 = I 2 = ∅, and so I 0 = {a} for a = 1, ..., N . Namely, reparametrising µ µ µ = {µ 1 , ..., µ a−1 , µ a+1 , ..., µ N } for convenience Since all the other terms in (38) have at least a multiplying factor D, we have In the rest of paper, we will use the sign ∼ l.d.
to indicate a low-density limit. Namely Using formula (30), one finds Hence the low-density limit (44) A few comments are in order. Up to now, nothing has been said of the Bethe equations, and in principle the µ i 's in this expression should satisfy exactly the Bethe equations (7). If one wishes to determine the dynamical correlations at leading order in D, there remains only the term (44) in the full expansion (38); but one can also satisfy the Bethe equations (7) only at leading order in D, since their exact solution will involve higher orders in D that are of the same order as the terms discarded in (38). Stated differently, the leading order in D of the dynamical correlations is obtained by both retaining only (44) and satisfying (7) at leading order in D, while the higher orders in D will require both taking into account higher terms in (38) and satisfying (7) at higher orders in D in (44).

Interpretation: low-density limit of the field form factor
The low-density limit is defined above as the leading term in an expansion of the correlation functions obtained by decomposing the form factors in partial fractions, that turns out to be an expansion in the density of the averaging state D. This definition allows for a systematic calculation of the next corrections in the density by taking into account more terms in (38).
The low-density limit of the correlation functions can however be recovered more intuitively but less rigorously with the following reasoning. If the root density ρ(λ) of the averaging state is small, then the distance between two consecutive roots L(λ i − λ j ) is 'typically' 2 large in front of 1, and in the limit of vanishingly low density becomes infinite. We will say that a sequence of states satisfying this property is dilute, namely For notational convenience we will drop the L dependence of the sequence of states. Let us consider a dilute state λ λ λ and investigate the consequences it has on the form factors between λ λ λ and another state µ µ µ.
Let us first note that because of the assumption of diluteness, a Bethe number J i of µ µ µ can be at a distance O(1) of at most one Bethe number I j of λ λ λ. Hence given a root µ i , the quantity L(µ i − λ j ) can be of order 1 for at most one λ j . Consequently, it is seen from the expression (12) that in the low density limit the form factor (12) is non-zero only if each Bethe number J i of µ µ µ is at a distance of order 1 from exactly one Bethe number I j of λ λ λ, and reciprocally if these N − 1 Bethe numbers I j are each at a distance of order 1 from exactly one Bethe number J k of µ µ µ. Since there are N roots in λ λ λ, there is one remaining root λ a that is not close to any µ i 's. We can re-label the roots µ µ µ = {µ 1 , ..., µ a−1 , µ a+1 , ..., µ N } so that L(µ i − λ i ) is of order 1 for all i = a. Let us then investigate the value taken by the normalized form factor (12) in this regime.
Let us first study the determinant N λ λ λ in the low-density limit. We have so that N λ λ λ = exp tr log G(λ λ λ) Here, we did not use the diluteness of λ λ λ, but only evaluated the leading order in D of the determinant.
Second, again because λ λ λ is dilute, all the L(µ i − λ i ) are negligible in front of any λ a − λ j . It follows that We also have and As for the matrix U , setting λ p = λ s = λ a , the dominant entries are while the other entries are of order O(L −1 ). Hence We obtain the following low-density limit of the form factor µ µ µ|ψ (0) |λ λ λ λ λ λ |λ λ λ µ µ µ| µ µ µ ∼ l.d.
Such an expression implies an ordering of the roots µ 1 , ..., µ N −1 according to the ordering of the λ i 's. Hence in the spectral sum (11), there is no factor 1/(N − 1)! once expressed in terms of the Bethe numbers of the µ j 's. Using this expression for the form factor, one indeed recovers the low-density correlation function (44) properly defined from the partial fraction decomposition, with N λ λ λ , N µ µ µ = 1 + O(D) already imposed at leading order in D.

Field two-point function
In this section we compute S 0,0,0 in (40) in the low-density limit.

States contributing to the thermodynamic limit
In order to carry out the sum over the Bethe numbers in (40), let us first investigate which values of Bethe numbers give a non-zero contribution in the thermodynamic limit.
Let us consider Bethe number configurations such that µ i − λ i = O(L −b i ) for i = a with b i ≥ 0. Since the parity of the number of roots λ i and µ i is different, the Bethe numbers of λ i are integers (resp. half-odd integers) if those of µ i are half-odd integers (resp. integers). Hence from the Bethe equations it follows that one has b i ≤ 1.  (40) is O(L −2N +1+2 i =a b i ). Hence the contribution of these configurations is O(L −N +1+ i =a b i ). Given that 0 ≤ b i ≤ 1, the only possibility to have a non-vanishing result in the thermodynamic limit is to have ∀i, b i = 1. Hence the only non-vanishing configurations in (40) are those for which the Bethe numbers of µ i differ from that of λ i by O(L 0 ). We will denote the difference between the Bethe numbers of µ i and λ i , with n i an integer of order O(L 0 ).

Decoupling of the spectral sum in the low-density limit
The Bethe roots µ i involved in (40) depend on the difference of Bethe numbers n i and are all coupled through the Bethe equations. Taking the difference of the Bethe equations (7) for µ k and λ k , one obtains with G = G(µ µ µ) introduced in (16), and the vector In the low-density limit, one has Gx ∼ l.d.
x for all vectors x, so that the roots decouple and are expressed as where we introduced Finally, at leading order in D, the determinants N λ λ λ , N µ µ µ are equal to 1 according to (47). The spectral sum (40) can thus be expressed as a product of N one-dimensional sums (61)

Thermodynamic limit of the spectral sum
In order to proceed we need to determine the thermodynamic limit of each of the onedimensional sums, that are of the type for α, w, τ reals, α non integer. Let us first consider the case τ = 0. The quantity n∈Z e iW n (n+α) 2 is exactly the Fourier series of a certain 2π-periodic function of W . Noticing that for n integer we conclude that To treat the case τ = 0, we write 14 The first term was computed previously, and the second one is 1/L times the Riemann sum of a regular function, hence in the thermodynamic limit converges to its integral. This yields with ± = sgn (τ ) and Hence we obtain with ± = sgn (t), using the expression (60) and trigonometric relations Plugging this expression in (61), we use to obtain in the thermodynamic limit This expression is valid for all space x and time t, and all interaction strength c > 0. It is the leading term in the expansion (38), and constitutes the low-density limit of the dynamical correlations of the field.

Density two-point function
In this Section we apply the same reasoning as in Sections 3 and 4 but to the density two-point function.

Partial fraction decomposition
We consider |λ λ λ and |µ µ µ eigenstates of the Lieb-Liniger Hamiltonian with N rapidities, and define F σ (λ λ λ, µ µ µ) by | µ µ µ|σ (0) |λ λ λ | 2 λ λ λ |λ λ λ µ µ µ| µ µ µ = F σ (λ λ λ, µ µ µ) Similarly to the field case F ψ , the reduced form factor F σ (λ λ λ, µ µ µ) is a ratio of two polynomials in the Bethe roots so that it can be decomposed into partial fractions. Repeating the arguments that apply to the field case, we similarly obtain the writing with S n,m,p = 1 2 n p!L 2N (75) The leading term is a priori obtained with ν i = 2, and the next terms by replacing some of the ν i 's by 1 or 0, each time picking a factor D. However, in contrast to the field case, the coefficient of the a priori leading term is actually found to vanish, using the analogue of (30) for the density operator case A(∅, ∅, ∅|∅) = 0 .
If one ν i is set to 1, the coefficient still vanishes However, if it is set to 0, the coefficient is non-zero and reads Hence the leading order in the density expansion, i.e. the low-density limit, is obtained with It yields the following low-density limit of the density two-point function (80)

Interpretation: low-density limit of the density form factor
Similarly to the field operator, the low-density limit of the density two-point function is defined in terms of a partial fraction decomposition of the form factor. In the field case, one was also able to obtain a low-density limit of the field form factor by considering dilute states, see Section 3.4. In the density case, we are going to follow a different route in order to show the usefulness of (53), by recovering the low-density limit of the density correlation (80) from the low-density approximation of the field form factor (53). The rationale is that the density form factor can be obtained as a limit of a field two-point function between different eigenstates λ λ λ and µ µ µ with N roots µ µ µ|σ(0)|λ λ λ µ µ µ|µ µ µ λ λ λ|λ λ λ = lim x→0 µ µ µ|ψ † (x, 0)ψ(0, 0)|λ λ λ µ µ µ|µ µ µ λ λ λ|λ λ λ .
The Bethe equations allow us to write with n j , p j integers. The integer p j is related to µ j and λ j through which is a parameter of the problem that is not to be summed over. We obtain the following factorization In order to determine the thermodynamic limit of the field four-point function we need to carry out the sums over the integers, that reduces to sums of the type We notice that n∈Z e iW n n+α is the Fourier series of a 2π-periodic function of W , and π −π π sin πα e iπα sgn (W )−iW α e −iW n dW = 2π n + α .
To use this formula we write .
Plugging this expression into the Lehmann representation for the density two-point function (11), one obtains exactly the low-density limit (80) previously obtained from the partial fraction decomposition. This provides an interesting consistency check of the interpretations of the low-density limit as obtained from low-density approximations of form factors.

The thermodynamic limit of the spectral sum
We now wish to compute the spectral sum (80) in the low-density limit. Repeating the arguments performed for the field case, in the thermodynamic limit only the states with µ j − λ j = O(L −1 ) contribute in (80). Besides, like for the field correlations, the Bethe equations decouple in the low-density limit. By taking the difference of the Bethe equations, we have for j = a with p j an integer. Using N λ λ λ , N µ µ µ ∼ l.d.
1, we obtain the following factorization of the spectral sum We now use (66) to compute the sum over the integers p. We find (95) Plugging this relation in (94), the product over j becomes an exponential in the thermodynamic limit, according to (69). The sum over the roots of the representative state λ a becomes an integral over the root density ρ(λ), while the sum over ν a that has to differ from any λ j becomes an integral over the hole density ρ h (λ) defined in (21). Hence one obtains (96) This expression is valid for all space x and time t, and all interaction strength c > 0. It is the leading term in the expansion (74), and is the low-density limit of the dynamical correlations of the density.

Bare particle-hole excitations
In the spectral sum (11) the intermediate states µ µ µ can be seen as excited states above the averaging state λ λ λ. In this picture it is natural to expand the spectral sum (11) in terms of the number of particle-hole excitations that µ µ µ has over λ λ λ. Such an expansion consists in writing with S n = µ µ µ |{Ia}∩{Ja}|=N −n | µ µ µ|σ(0)|λ λ λ | 2 λ λ λ |λ λ λ µ µ µ| µ µ µ e it(E(λ λ λ)−E(µ µ µ))+ix(P (µ µ µ)−P (λ λ λ)) , which is the spectral sum restricted to intermediate states µ µ µ that share N − n Bethe numbers J a with the Bethe numbers I a of λ λ λ. Ideally, each individual S n would have a well-defined and finite thermodynamic limit L → ∞ that could be represented as a multiple integral over root densities and hole densities The function F n (λ 1 , µ 1 , ..., λ n , µ n ) appearing in this expression would thus be identified as a thermodynamic form factor for n particle-hole excitations [51][52][53][54]. This idea was backed by calculations of lim that is finite and well-defined in the thermodynamic limit, by various means [52,56]. However, this picture seems to be a priori contradicted by results of [58], where the thermodynamic limit of the spectral sum (11) was computed exactly at order c −2 . At this order, the spectral sum exactly truncates to one-and two-particle-hole excitations, hence to S 1 and S 2 . But it was found that these two separate sums in the thermodynamic limit individually diverge and depend on the representative state λ λ λ of the root density, at order c −2 . Namely we have with A 1,2 some reals of order c −2 and f 1,2 (ρ, γ −2 ) functions of the root density and the pair distribution function (see [58] for a precise definition). However, their sum S 1 + S 2 is not divergent and depends only on ρ, as we expect. To fit in the general picture (99), the only solution would be that these divergences in L are an artefact of the 1/c expansion, i.e. that the finite-size corrections to S 1,2 in the thermodynamic limit would be e.g. of the form ϕ(L/c) with a function ϕ(x) → 0 when x → ±∞. In the framework of the low-density expansion, one can investigate this issue since the calculations are performed non-perturbatively in 1/c, and compute the thermodynamic limit of S 1 for a finite fixed c, in the low-density limit.
and assume that µ µ µ only involves one particle-hole excitation above λ λ λ. Hence for all j = a, the Bethe numbers J j of µ µ µ are equal to the Bethe numbers I j of λ λ λ. In the low-density limit we have then ∀j = a, It follows that in the thermodynamic limit, the low-density form factor squared (102) becomes | µ µ µ|σ(0)|λ λ λ | 2 µ µ µ|µ µ µ λ λ λ|λ λ λ ∼ l.d. with The sign of the logarithm can be determined by applying the following inequality that is valid for all real u u 2 ≥ sin 2 u , to the value u = arctan x − arctan y , for x, y reals. Indeed, using trigonometric relations, Eq (106) is exactly with equality only for x = y. This permits to deduce that the logarithm in (105) is always negative whenever λ = µ. One concludes that we always have if λ = µ φ(λ, µ) < 0 , and if λ = µ we have It follows that in the framework of a multiple integral representation in terms of bare particlehole excitations in the thermodynamic limit (99), the function F 1 (λ, µ) would be in fact zero everywhere except if µ = λ where it takes the value 1 as deduced from (102). Namely we would have This analysis shows that the integral representation (99) in terms of 'bare' particle-hole excitations (98) is singular. If the leading behaviour at large space and time of (99) is indeed obtained when λ 1 = µ 1 in (99) for n = 1, where the function F 1 (λ, µ) takes a finite value, the singular behaviour of this representation should be an obstacle to the computation of the subleading orders.

Dressed particle-hole excitations
Let us now investigate the nature of the states summed over to obtain the expression (96) in the low-density limit. The double integral over the root density and hole density are oneparticle-hole excitations with a macroscopic amplitude (i.e. the difference between the Bethe numbers is O(L)). The exponentials however arise from the product of a macroscopic number of one-dimensional sums over all the other remaining roots of the averaging states. These take into account an arbitrary number of particle-hole excitations, but with only a microscopic or mesoscopic amplitude (i.e. the difference between the Bethe numbers can be O(L ν ) for any ν < 1). These configurations are represented in Figure 1. They correspond to what is called a dressed one-particle-hole excitation. They entail expanding the spectral sum (11)  | µ µ µ|σ(0)|λ λ λ | 2 λ λ λ |λ λ λ µ µ µ| µ µ µ e it(E(λ λ λ)−E(µ µ µ))+ix(P (µ µ µ)−P (λ λ λ)) .
. . . . . . . . . . . . . Figure 1: Sketch of a dressed one-particle-hole excitation: positions of the momenta of the averaging state λ λ λ (empty circles) and the intermediate state µ µ µ (filled circles) respectively, and position of the holes (dots). In red is indicated the 'soft modes' corresponding to microscopic excitations, and in blue the only macroscopic excitation.
Loosely speaking, S dr n includes n macroscopic particle-hole excitations and any number of microscopic particle-hole excitations. Since these configurations are not disjoint, one requires the expression (113) for a precise definition. In the low-density limit the full spectral sum truncates to S dr 1 , and this expression is found to be finite and well-defined in the thermodynamic.
We note that these dressed particle-hole excitations are actually in principle what is computed in [51][52][53]. But therein the 'soft modes' contribute as a numerical factor that multiplies the form factor. Here, within the low density expansion, it is seen in (96) that these soft modes actually carry an x and t dependence as well. Namely we have a representation for S dr 2n F x,t n (λ 1 , µ 1 , ..., λ n , µ n )e it(E(λ λ λ)−E(µ µ µ))+ix(P (µ µ µ)−P (λ λ λ)) × ρ(λ 1 )ρ h (µ 1 )...ρ(λ n )ρ h (µ n )dλ 1 dµ 1 ...dλ n dµ n , where the 'dressed thermodynamic form factor' F x,t n (λ 1 , µ 1 , ..., λ n , µ n ) carries a x, t dependence coming from the soft modes summation, although the soft modes do not modify the energy and momentum of the state in the thermodynamic limit. This dependence emerges from the double poles that lift to order L 0 some finite-size corrections. The function F x,t 1 (λ, µ) can be directly deduced from (96).
From Sections 5.4 and 5.5, we conclude that the right expansion scheme of the spectral sum is an expansion in terms of dressed particle-hole excitations (112), in the sense that the truncated spectral sums S dr n are finite and well-defined in the thermodynamic limit, and are smooth functions of the macroscopic excited rapidities. Moreover, as shown in [58], they admit a well-defined and uniform 1/c expansion without any spurious divergences in L coming from c-dependent finite-size effects. These two properties are not satisfied by the 'bare' particle-hole excitations expansion.

Summary and conclusion
We computed the dynamical two-point function of the field and density operator averaged within a state with a small particle density D, given by Equations (70) and (96). They are valid for an arbitrary interaction strength c > 0 and for all space and time -hence one can deduce the spectral function and dynamical structure factor in the full momentum-frequency plane, in the same regime of small D. This low-density limit is defined as the leading term in the correlation functions obtained by decomposing the form factors into partial fractions.
Besides the explicit expressions obtained in the low-density regime, this work also provides interesting insights on the nature of states that contribute to the spectral sum, and on its possible expansions as detailed in Sections 5.4 and 5.5. The low-density regime is indeed naturally interpreted as a single dressed particle-hole excitation, i.e. a macroscopic particlehole excitation with an arbitrary number of microscopic particle-hole excitations. In contrast, an integral representation in terms of 'bare' particle-hole excitations in the thermodynamic limit is found to be singular, in the sense that the 'thermodynamic form factor' for oneparticle-hole excitations is non-zero only in the limit where the amplitude of the macroscopic excitation vanishes. Hence this work indicates that the right expansion of the spectral sum is in terms of dressed particle-hole excitations, as already suggested by the 1/c expansion developed in [58]. Importantly, this dressing by microscopic particle-hole excitations also comes with an x and t dependence.
This PFD framework allows in principle for a computation of the next orders, which constitutes the most natural direction of improvement of this work. The fact that the next orders can indeed be computed with the PFD has been shown in [10] for a model that can be reformulated in terms of free fermions. In the Lieb-Liniger case, the interaction introduces technical but not fundamental difficulties, and we hope to be able to pursue this program in following works.