A systematic $1/c$-expansion of form factor sums for dynamical correlations in the Lieb-Liniger model

We introduce a framework for calculating dynamical correlations in the Lieb-Liniger model in arbitrary energy eigenstates and for all space and time, that combines a Lehmann representation with a $1/c$ expansion. The $n^{\rm th}$ term of the expansion is of order $1/c^n$ and takes into account all $\lfloor \tfrac{n}{2}\rfloor+1$ particle-hole excitations over the averaging eigenstate. Importantly, in contrast to a 'bare' $1/c$ expansion it is uniform in space and time. The framework is based on a method for taking the thermodynamic limit of sums of form factors that exhibit non integrable singularities. We expect our framework to be applicable to any local operator. We determine the first three terms of this expansion and obtain an explicit expression for the density-density dynamical correlations and the dynamical structure factor at order $1/c^2$. We apply these to finite-temperature equilibrium states and non-equilibrium steady states after quantum quenches. We recover predictions of (nonlinear) Luttinger liquid theory and generalized hydrodynamics in the appropriate limits, and are able to compute sub-leading corrections to these.


Introduction
The Lieb-Liniger model [1] is a key paradigm of integrable many-particle systems [2]. Moreover, it is directly relevant to a range of cold atom experiments both in and out of equilibrium, see e.g. [3][4][5][6][7][8]. While the excitation spectrum at zero temperature [9] and thermodynamic properties [10] have been known for a long time, the exact solution does not provide easy access to correlations functions as these encode more detailed information about the exact energy eigenstates. An exception is the case of impenetrable bosons [11][12][13][14][15][16][17][18][19][20][21][22][23][24], which can be mapped onto non-interacting fermions. In absence of full analytic solutions valuable insights on the large space and time asymptotic behaviours of correlation functions at zero and low temperatures were gained by combining exact results on spectral properties obtained from the Bethe Ansatz with with conformal field theory (CFT) [25,26] and Luttinger liquid theory [27,28] and its recent extensions [29][30][31][32][33][34]. The last two decades then witnessed remarkable progress in the computation of zero temperature dynamical correlation functions by expressing them in terms of spectral representations over the energy eigenstates of the model. On the one hand it became possible to numerically evaluate the spectral sums to very high precision for large, finite systems [35,36]. On the other hand remarkable analytic progress led to a fairly complete understanding of the asymptotic behaviour at late times and large distances [37][38][39].
To make analytical progress it is essential to identify the classes of states that give the dominant contributions in a given range of frequencies and momenta or space and time [37]. Known results suggest that in interacting theories this generally requires the summation over an infinite number of states. Firstly, the large space and time asymptotics of zero temperature dynamical correlators in interacting models has been shown to be determined by an arbitrary number of (soft) particle-hole excitations over the ground state around the Fermi points and the saddle points of the dispersions of elementary excitations [37,75]. Secondly, it has been shown that the asymptotic behaviours of dynamical correlations of semi-local operators in thermal and other finite entropy states involves an arbitrary number of (soft) particle-hole excitations [76] over the macro state of interest. Truncating this sum to a finite number of particle-hole excitations leads to a result that diverges in time. In the zero temperature case it has been shown that it is possible to take the thermodynamic limit of (partial) spectral sums and obtain a representation in terms of (dressed) excitations in the thermodynamic limit [37,75]. An analogous result for the finite temperature/entropy case would be highly desirable, but is not known at present. In Refs [70,72,73] such an expansion in terms of thermodynamic particle-hole excitations was conjectured. It is based an phenomenological assumptions on how partial sums over states in the finite volume combine into thermodynamic form factors. It also exhibits singularities, whose regularization is not presently known.
Given this state of affairs it is highly desirable to obtain explicit results through ab initio calculations that do not require any assumptions, i.e. carrying out the spectral sum in a finite volume and then taking the thermodynamic limit exactly. In order to make progress in this direction we consider the spectral sum in the framework of an expansion in the inverse interaction strength c −1 around the impenetrable limit. Strong coupling expansions have have previously been used at zero temperature and for static correlators at finite temperatures [77][78][79][80]. More recently the 1/c contribution for the finite temperature dynamical densitydensity correlation function was determined in [69]. This contribution has a particularly simple structure similar to that of the impenetrable limit, that does not carry over to the next orders, and as a consequence until now it has been unclear how to determine higher orders in this expansion. In the following we develop a method for calculating the higher orders of this expansion and apply it to obtain the contribution to the dynamical densitydensity correlator at order c −2 . The general idea of the 1/c expansion, and more generally of strong coupling expansions in integrable models, is as follows. A consequence of integrability is that N -particle energy eigenstates in a finite volume can be labelled by N rapidity variables |λ = |λ 1 , . . . , λ N . (1) These rapidities are in a one-to-one correspondence with sets of (half-odd) integers {I j } through the quantization conditions in the finite volume The energy and momentum of these states are given by where (λ) and p(λ) parametrize the energy and momentum of a single-particle excitation over the vacuum (reference) state. For the Bose gas we have (λ) = λ 2 and p(λ) = λ. Twopoint correlation functions of a local operator O(x) in a given energy eigenstate |λ thus have spectral representations of the form where the first sum runs over the particle number and the second over all M -particle energy eigenstates. The matrix elements F O (λ, µ) = λ|O(0, 0)|µ λ|λ µ|µ (5) are also known as form factors and, as we will see, admit a 1/c-expansion Denoting the truncation of the sums to order O(c −j ) by F We stress that the expansion sums certain 1/c contributions to all orders by virtue of the fact that although the (exactly known) energies and momenta are expanded inside the exponentials, the exponentials are not expanded in 1/c. In this sense the expansion is nonperturbative, and in fact rather different from more standard (diagrammatic) approaches pursued in [81]. As discussed in detail below (8) is in fact both a 1/c expansion and an expansion in terms of number of particle-hole excitations. At order n in the expansion (i) only excitations that involve at most n 2 + 1 particle-hole pairs contribute, and (ii) all terms up to O(c −n ) contribute. Importantly, this "mixed" expansion has a well-defined thermodynamic limit and is uniform in space and time. This is in contrast to both the bare 1/c expansion that is non-uniform, and the bare expansion in the number of particle-hole excitations that is divergent in the thermodynamic limit.
Expectation values of the form (4) are relevant in two contexts.
1. By working in a micro-canonical ensemble dynamical response functions at finite temperature can be cast in this form. In the following we will use this to determine the finite temperature dynamical structure factor in the Lieb-Liniger model.
2. At late times after quantum quenches local observables relax to non-thermal stationary values [91][92][93][94]. It follows from the quench action approach [89,90] to quantum quenches that expectation values in the stationary state in fact involve non-thermal energy eigenstates at finite energy densities. This has been used to study the stationary behaviour of certain one-point functions after (particular) quantum quenches [95][96][97]99]. A natural extension is then to consider linear response functions in such steady states [100,101]. These can be expressed in the form (4), where |λ corresponds to the non-equilibrium steady state relevant to the quench of interest.
In the following we will consider both these cases and evaluate (4) for the density operator and general |λ . A brief summary of some of our key technical results is as follows. We show that the 1/c-expansion corresponds to an expansion in the number of particle-hole excitations. This leads to a dramatic reduction in the complexity of the spectral sum that needs to be carried out. Interestingly, the contributions of one particle-hole and two particle-hole excitations are individually divergent in the infinite volume limit L → ∞. Moreover they individually depend on details of the "averaging state" |λ beyond the root distribution function in the thermodynamic limit. Crucially, their sum is not divergent and is independent of the choice of representative state |λ , and is well-defined.
The manuscript is organized as follows. In Section 2 we introduce the Lieb-Liniger model and recall the key elements of its Bethe Ansatz solution. In Section 3 we report some important intermediate results on the thermodynamic limit of expressions computed within the Bethe Ansatz. In Section 4 we discuss the 1/c-expansion up to and including O(c −2 ) of the Bethe Ansatz equations, energy eigenvalues, form factors and the spectral representation of the density-density correlation function. These results are then used in Section 5 to obtain a fully explicit expression for the dynamical density-density correlator (and the related dynamical structure factor) in the thermodynamic limit, cf. equations (171), (166), (178) and (187). This constitutes the main result of our work. In Section 7 we obtain the asymptotic behaviour of the correlator and structure factor in various regimes. In particular we perform non-trivial consistency checks of our formulas, and recover known results from (nonlinear) Luttinger liquid theory and generalized hydrodynamics (GHD) [83,84,88].

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 In the following we set = 2m = 1 and impose periodic boundary conditions. In first quantization (9) corresponds to a quantum mechanical system of N particles with positions 0 ≤ x 1 , ..., x N ≤ L and Hamiltonian For later convenience we 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: the energy E and the momentum P of an eigenstate |λ λ λ with N bosons read where the rapidities λ λ λ = {λ 1 , .., λ N } satisfy the following set of "Bethe equations" It is convenient to express them in logarithmic form with I k an integer if N is odd, a half-integer if N is even. For c > 0, which we will assume in this paper, all the solutions to this equation are real [2].
The (normalized) form factors of local operators between two Bethe states have been derived in Refs [102][103][104]. In the case of the density operator σ, the (square of the normalized) form factor between two eigenstates |λ λ λ , |µ µ µ with respective numbers of Bethe roots N, N reads Here p ∈ {1, ..., N } can be freely chosen, and N λ λ λ is given by 3 Thermodynamic description of eigenstates In a finite system of size L all eigenstates of the Hamiltonian are fully characterized by a set of N Bethe numbers I k , or equivalently a set of N Bethe roots λ k . The purpose of this section is to explain how to turn this description into one based on (continuous) distribution functions of these roots in the thermodynamic limit L → ∞ when N scales like L. In particular, contrary to a common misconception, we emphasize that the usual "root densities" defined below do not fully characterize an eigenstate in the thermodynamic limit; this observation turns out to be of crucial importance in our calculation.

Root density
In the thermodynamic limit, any sum of a non-singular (piece-wise continuous) function f over the Bethe roots or Bethe numbers is independent of the precise values taken individually by each I k or λ k , and depends only on the number of Bethe roots or Bethe numbers in any given interval. This information is encoded in the so-called root density ρ(λ) ≥ 0 and filling function 0 ≤ χ(ι) ≤ 1. They are defined by the requirement that in the large L limit Lχ(ι)dι = number of Bethe numbers I k /L per length in [ι, ι + dι].
In the thermodynamic limit the sums (21) can be turned into integrals over these functions The same holds for multidimensional sums of a multivariate non-singular function f , with converging to As far as expressions of the form (21) and (24) are concerned, an eigenstate in the thermodynamic limit is entirely characterized by the root density ρ(λ), or equivalently the filling function χ(ι). To relate these two equivalent quantities, we introduce the function ϑ(λ) as the L → ∞ limit of a function ϑ(λ k ) ≡ χ( I k L ) of the Bethe roots, where I k is the integer associated with λ k . Using the Bethe equations (15) ϑ(λ) can be expressed in terms of χ and ρ as The filling function χ(ι) and the root density are then related through It is customary to introduce the so-called hole density ρ h (λ) defined by which again contains equivalent information to ρ(λ) or χ(ι). When expressed in terms of the particle and hole densities (27) is known as the thermodynamic limit of the Bethe Ansatz equations [105]. Finally, the particle density is given by We introduce the Fermi momentum q F defined by Although there is a simple relation between D and q F , we will in the following sometimes use D and sometimes q F , depending on the physical context at hand. We also denote (in the units where = 2m = 1)

Definition
Root densities entirely characterize the value of sums of the type (21) and (24) in the thermodynamic limit. However, some functions of the Bethe roots cannot be expressed solely in terms of root densities in the thermodynamic limit, and as a consequence can take different values in the thermodynamic limit for states that have the same root density. An example is provided by that we will encounter below 1 . The sum in (32) by definition depends on the joint distribution function of pairs of roots separated by O(L −1 ), and the latter clearly contains information beyond that contained in the root density (which does not distinguish between roots separated by O(L −1 )). We first note that if we impose the constraint |λ i − λ j | > for a > 0 then Σ L [g] vanishes in the thermodynamic limit. Hence, it only depends on g(λ, λ) and its derivatives at λ. Taylor expanding g(λ i , λ j ) for λ i close to λ j reduces the order of the pole and makes the next terms vanish in the thermodynamic limit, so it depends only on g(λ, λ). Being a linear functional of g it can be written in the thermodynamic limit in the form where the function γ −2 (λ) depends on the state. We call γ −2 (λ) a pair distribution function as it encodes information about the joint distribution of pairs of Bethe roots. The index −2 relates to the fact that we are summing over the inverse square of the difference between two Bethe roots. The pair distribution function γ −2 (λ) characterizes certain properties of the thermodynamic limit of an eigenstate and is unrelated to the root density ρ(λ). Two states can have the same ρ(λ) but different γ −2 (λ). The simplest example is that of (translationally invariant) free fermions, where the Bethe roots reduce to the single-particle momenta. Here we may construct two sequences of eigenstates labelled by an integer n, with momenta {λ i = ni L |i = 1, . . . , N } and {µ 2i = 2ni L , µ 2i+1 = 2ni+1 L |i = 1, . . . , N/2} respectively. In the thermodynamic limit both states are described by a root density ρ(λ) = 1/n, but the pair distribution functions are different: γ −2 (λ) = π 2 3n 2 for the first state and γ −2 (λ) = 1 + π 2 12n 2 + m =0 1 (2nm+1) 2 for the second one.

(Generalized) micro-canonical ensemble and representative states
The (generalized) micro-canonical ensemble average of a local operator O(x) is a priori defined as where the sum is over an appropriate "shell" of simultaneous eigenstates of the Hamiltonian and the local conservation laws of the theory. C L is the number of terms in the sum. In a large but finite volume this means that for thermal averages we fix the energy within a window that contains an exponential (in system size) number of eigenstates. In the case of generalized micro-canonical ensembles we fix the eigenvalues of (some or all) of the local conservation laws in an analogous fashion [89,106]. It is believed that almost all states in the sum in (34) have identical local properties, and hence the sum over states can be replaced by an expectation value with respect to a single typical state |λ in the thermodynamic limit The state |λ is sometimes called a representative state and we follow this terminology here. We note that in practice there is a great deal of freedom in choosing a representative state in a large, finite volume.

Average over representative states
As we have seen above, the thermodynamic limit of the sum (32) cannot generally be expressed as an integral over the root density, but depends on the choice of representative state in the finite volume. The thermodynamic limit of these sums involves the separate function γ −2 (λ) defined in (33). As we will see in the following, in our calculations of the density-density correlation function the dependence of certain intermediate quantities on γ −2 (λ) eventually compensate and the end result depends only on the root density. However, it is a priori possible that in other calculations involving sums of form factors no such cancellations will occur and the end result will indeed depend on the choice of representative state through γ −2 (λ) or an analogous quantity.
We now make the following observation. As we have discussed above, averages with respect to a Bethe state |λ λ λ often emerge upon simplifying averages over exponentially (in system size) many representative states corresponding to a given root density ρ(λ). By construction such averages will depend only on the density. This then poses the question what value (32) takes after averaging over all representative states with same root density in the thermodynamic limit. We now address this issue.
First, we need to define properly which states in a large finite volume L are acceptable representative states for a given root density. We define a sequence of sets of states to be complete for the root density ρ if the corresponding sequence of sets of solutions to the Bethe equations (S L ) L∈N all give rise to the density ρ in the thermodynamic limit, and if the number of elements of the set S L satisfies where S YY [ρ] is the Yang-Yang entropy [10] In order to build such a set in a large finite volume let us consider a root density ρ(λ) at a given particle density D = ρ(λ)dλ. Given the root density we may introduce a particle counting function by Next we choose a "coarsening function" L with the property that L → 0 and L L → ∞ when L → ∞ -for example one can take L = 1 √ L . We now split the real axis into n L "bins" [x L,j , x L,j+1 ] containing L L Bethe roots by defining x L,1 , ..., x L,n L +1 such that z(x L,i ) = i L for 1 ≤ i ≤ n L + 1 = D/ L .
Finally we define S L as the set containing all the states in a finite volume L that contain exactly L L Bethe roots in each of the n L bins [x L,i , x L,i+1 ]. All states in S L have N L = ( D/ L −1) L L Bethe roots, which for L → ∞ by construction are distributed with density ρ(λ). The number of elements of S L will depend on the number of "vacancies" in each of the bins, which in turn depend on the values of all the Bethe roots since they interact via the Bethe equations. However, asymptotically in L, we have K L,i = L(x L,i+1 −x L,i )(ρ(x L,i )+ρ h (x L,i )) vacancies in each of the bins, so that Using Stirling's formula in the large-L limit one has which shows that S L is indeed a complete set of representative states for a given root density ρ(λ).
We can now state our result for the average of (32) over all representative states with root density ρ(λ): A proof of (41) is given in Appendix B.
We note that if we instead sum over rapidities distributed regularly according to the inverse of the counting function z −1 (λ) without imposing that the rapidities are solutions of the Bethe equations, the sum takes a different value: If we sum over rapidities distributed regularly according to the inverse of the counting function z −1 (λ) and impose the Bethe equations, the sum (32) is not easily expressed in terms of ρ, but takes a value different from either (42) or (41). Hence formula (41) is both non-trivial and non-intuitive.

Single principal value
The sums (24) can be expressed in terms of root densities in the thermodynamic limit, provided f is non-singular. We have seen in the previous section that for functions f with a quadratic singularity the thermodynamic limit value of the sum cannot be expressed in terms of the root density. We now turn to functions that are singular but integrable in a principal value sense. This is the case of the sum We will assume that g and ρ are continuous. Symmetrizing the sum, we havẽ The function F (x, y) = g(x,y)−g(y,x) x−y is regular, so that it has the form of (24) and its thermodynamic limit can be expressed in terms of ρ according tõ Since the integrand is finite, one can remove a small shell |λ − µ| < with an error of O( ), and then un-symmetrize the sum. This yields with the following usual definition of the principal value integral Hence, sums of type (43) can indeed be expressed in terms of root densities. In contrast partial sums like 1 L at fixed j cannot be expressed in terms of the root density in the thermodynamic limit.

Double principal values
Higher-dimensional sums of the form can be treated likewise, but with subtleties hiding in the fact that i can be equal to k. Separating out the term with i = k and symmetrizing the remaining sum gives (50) The first term is regular so that (25) can be used, while the second term is of the type (32) and can be expressed in terms of γ −2 (λ). In the first term we can remove the region where |λ − µ| < or |λ − ν| < or |ν − µ| < with an error that is O( ), and then un-symmetrize the integral. One obtains where the simultaneous principal value in the triple-integral is defined as As shown in Appendix A.1 this can be expressed in terms of the successive principal value triple-integral according to a Poincaré-Bertrand-like formula where we defined It can also be expressed as and as shown in Appendices A.2. Using these principal value integral identities we can rewrite (51) in the form

Examples of root densities
The calculations presented in this paper hold for a generic piece-wise continuous root density ρ(λ). Two applications we have in mind is to thermal states and non-equilibrium steady states after quantum quenches, and we now discuss specific root densities that arise in these contexts.

Thermal states
Thermal states are characterized by root densities that maximise the Yang-Yang entropy at inverse temperature β [10]. Defining the so-called dressed energy ε dr (λ) by the filling function n(λ) of a thermal state is such that Here h is a chemical potential that is used to fix the desired particle density D. In practice one first solves the nonlinear integral equation (59) and then uses (58) to determine ρ(λ) from the linear integral equation (27). A particular case of thermal states is the zero temperature ground state, obtained in the limit β → ∞. Its root density satisfies with Q defined such that

Non-equilibrium steady states
Refs [95,97] considered a particular interaction quench in the Lieb-Liniger model, where the system is initially in the ground state of (9) for c = 0, and is subsequently time-evolved with the Lieb-Liniger Hamiltonian at a finite value of c. The root density characterizing the steady state reached at late times was determined in [97] and remarkably allows for a closed form solution

1/c expansion of the Lieb-Liniger model
In this section we perform an expansion around the limit c → ∞ at order 1/c 2 of the energy levels and form factors in the Lieb-Liniger model, at fixed L and fixed Bethe numbers. We then expose the consequences it has on the spectral sum (17) in Section 4.3.3.

The Bethe equations
The Bethe equations (15) admit a regular 1/c expansion at large c. In the following, in order to expand the form factor at order 1/c 2 we will need the value of the Bethe roots at order 1/c 3 . The Bethe equations (15) at order 1/c 3 read This gives the following expression for the Bethe roots in terms of the Bethe numbers at order 1/c 3 The alert reader will have noticed that some of the terms contain higher powers of 1/c than the order at which we are working, that is 1/c 3 . We find it useful throughout the manuscript to retain certain "resummed" expressions of 1/c as they appear in calculations, both for clarity and convenience since they often happen to compensate each other. In any case, keeping these resummed expressions in 1/c does not affect the validity of the equations at the order considered.

Leading order in 1/c of the form factor between two generic states
The behaviour of the 1/c expansion of a form factor (18) between states |λ λ λ and |µ µ µ depends on the "relative positions" of the Bethe numbers of one state to the other. To see this, let us determine the leading order in 1/c of the form factor (18) without making any assumptions on the eigenstates |λ λ λ and |µ µ µ . It is then straightforward to see that when c → ∞ while the non-diagonal term in the determinant appearing in the form factor is of order We see that the order in 1/c of this expression entirely depends on the roots λ k and µ k . To be specific, let us now denote by I k and J k the Bethe numbers of λ λ λ and µ µ µ respectively, and define the number of Bethe numbers present in λ λ λ and absent from µ µ µ. If λ i and µ j have different Bethe numbers, then from (65) Hence expanding in 1/c naturally orders the Lehman representation (17) into an expansion in terms of number of particle-hole excitations of µ µ µ above λ λ λ, i.e. of the number of changes in the Bethe numbers of µ µ µ compared to those of λ λ λ. This means that if one considers (17) at order c −m , then only intermediate states µ µ µ with ν ≤ m 2 + 1 contribute to the sum. We note however that the converse is not true: restricting (17) to e.g. one-particle-hole excitations would still involve arbitrarily high orders in 1/c. Since our goal is to compute correlations at order 1/c 2 , we only need to investigate the restriction of (67) to one-and two-particle-hole excitations.

Order 1/c 2 of form factors involving a single particle-hole excitation
In this section we consider one-particle-hole excitations of the state |µ µ µ above |λ λ λ . Up to reordering the roots, we can assume that the Bethe numbers I k of λ λ λ differ from those J k of µ µ µ only at a single position a: Since the excited particle cannot coincide with an already existing particle, we also have the constraint ∀i = a I a + n = I i .
This has the following consequences at order 1/c 3 on the value of the Bethe roots. Using (65) we have while for i = a, we obtain Using (72) and (73) we can determine the various terms entering the expression of the form factor at order c −2 i =j Putting everything together we have at order c −2

Order 1/c 2 of form factors involving two particle-hole excitations
We now consider two particle-hole excitations. Up to re-ordering the roots of µ µ µ, we can assume its Bethe numbers differ from those of λ λ λ only at positions a and b = a, and thus assume Since the excited particles cannot coincide with an already existing particle, we also have the constraints Moreover we must also exclude the case where one of the excited particles fill the hole left by the other, since this reduces to a single particle-hole excitation and is therefore already covered. The corresponding constraint is Finally we have to exclude the case where the two excited particles coincide From (65) we obtain We can now investigate the form taken by (67) for these values of roots. At leading order in 1/c we have which, when substituted in (67) yields the following leading order expression of the form factor for two-particle-hole excitations

The Lehmann representation
We can now write the Lehmann representation (17) for the density-density correlation functions at order 1/c 2 . As explained in the previous section, only one and two particle-hole excitations contribute to (17) at order 1/c 2 , and the corresponding form factors were computed at this order in the previous subsections. This leaves us with working out the phases in the corresponding terms in (17) at order 1/c 2 .

The phase for a single particle-hole excitation
For excitations with one particle and one hole, it follows from (72) and (73) that It will be convenient to perform the following change of variable x defined as where Then the phase becomes For later convenience we define 4.3.2 The phase for two particle-hole excitations We can express this in terms of x as well

The sum over intermediate states
So far we have expanded all the terms arising in (17) at order 1/c 2 , at a fixed L for arbitrary eigenstates |λ λ λ and |µ µ µ with fixed Bethe numbers. We have shown that the sum truncates to one-and two-particle-hole excitations, and that the resulting terms are well-defined functions of the excitation parameters n and m. However, as the Lieb Liniger model is a field theory and not a lattice model it features an infinite number of particle-hole states even if L is finite, so that (17) is still an infinite sum even if it involves only one-and two-particle-hole excitations. This creates two notable problems. The first one is that we encounter infinite sums of the type ∞ k=−∞ k n e ik 2 t+ikx for n = 0, 1, 2 which are ill-defined as functions of x, t (except if n = 0 and t = 0). The explanation for this behaviour is that σ (x, t) σ (0, 0) , similarly to the propagator of a quantum particle, should be understood as a probability amplitude that is meant to be integrated against a smooth and localized function of x and t, or, stated differently, that it must be understood as a distribution in x, t. The second problem is that the 1/c expansion of a form factor λ λ λ|σ(0)|µ µ µ has been performed for fixed Bethe numbers, whereas in the spectral sum at fixed c there are always excited states with Bethe roots larger than c. This poses a potential problem of commuting two limits.
In order to address these problems we are going to impose that all the rapidities involved in the spectral sum (17) are smaller than a certain cut-off Λ, that can be taken as large as desired. Firstly, this imposes a restriction of the state |λ λ λ , in which we are calculating our expectation value. We require that for all roots |λ j | < Λ, i.e. that the density ρ(λ) vanishes for |λ| > Λ; this is a mild restriction in the following sense. In practice we are interested in the dynamical response in macro states characterized by root distributions ρ(λ) that decay faster than |λ| −2 for λ → ∞, which is a necessary condition for the energy density of the state to be finite (for example in the thermal state the decay is Gaussian). We therefore can always approximate ρ(λ) to any given accuracy by a root density ρ Λ (λ), which vanishes outside the interval [−Λ, Λ]. Moreover this truncation can be done in an infinitely differentiable way, so it does not affect the regularity of the root density ρ(λ). Secondly, this cut-off also restrains the sum (17) to excited states |µ µ µ such that the |µ i | < Λ, which removes the problem of possible excited rapidities becoming larger than c. Hence we define a Λ-regularized correlation function The correlator σ (x, t) σ (0, 0) Λ defined in this way and expanded in 1/c has a regular thermodynamic limit L → ∞, as we will see below. Now, in order to recover the true correlation functions (17), one would like to then take the limit Λ → ∞. It turns out that such a limit of σ (x, t) σ (0, 0) Λ seen as a function of x, t does not exist. To be more specific one encounters problematic terms of the form for which the limit Λ → ∞ does not exist (except for n = 0 if t = 0). However, the limit exists in a distribution sense, i.e. the integral of I n (Λ|t, x) over any smooth localized function of x, t has a well-defined limit when Λ → ∞. This is all we require, since the correlation function is in any case meant to be integrated with a smooth localized function of x, t.
To take the limit we perform an integration by part and obtain In particular we have where Terms like e −itΛ 2 ∓ixΛ Λ n and e ixΛ do not have limits when Λ → ∞ as a function of x, t. In a distribution sense however, they vanish when Λ → ∞ in the sense that their integral with any smooth localized function of x, t vanishes when Λ → ∞. Hence we obtain that when Λ → ∞ I n (Λ|t, x) tends to I n (t, x) with I n (0, x) = 0 and One notices that these limits are exactly those obtained by introducing a small imaginary part in time and taking Λ → ∞ However, such a small imaginary part cannot be incorporated from the beginning in (91), since E(λ λ λ) − E(µ µ µ) can take both signs when |λ λ λ is not the ground state. These limits will be useful in the following sections in order to take the limit Λ → ∞ of the Λ-regularized correlation functions.
At order 1/c 2 we therefore have the following decomposition where C Λ 1,2 (x, t) are defined in the following. Introducing the convenient notations and we have the following contribution at order c −2 of the one-particle-hole excitations We already neglected a global factor (1 + 2 cL ) 2 that is 1 in the thermodynamic limit, as well as a O(L −1 ) contribution in the exponential. We also used 1 In Figure 1 we show the distribution of Bethe numbers for the particle-hole excitations that are summed over in (101). Compared to the representative state we have changed a single integer. For the two particle-hole excitations the sum in (17) is over the set {a, b} and over n, m with the constraint µ a < µ b . Since the form factor is symmetric upon swapping a, b and n, m simultaneously, this constraint can be taken into account with a factor 1/2 and with imposing µ a = µ b . The sum over {a, b} as a set can be transformed into a sum over a = b as a couple with a factor 1/2 as well. Hence we have the leading contribution of the two particle-hole In Figure 2 we show the distribution of Bethe numbers for the two particle-hole excitations that are summed over in (101). Compared to the representative state we have changed two integers.

Examples of root densities
In this subsection we complete the 1/c expansion of the model by determining the expansions of the root densities introduced in Section 3.4.

Hole density
We introduced earlier the hole density ρ h (λ) in (28) with ϑ(λ) given in terms of ρ(λ) in (27). The hole density is a function of the root density ρ(λ), and for a generic ρ(λ) it reads at order where we recall that D is defined in (29).

Thermal states
Thermal states at finite inverse temperature β < ∞ are defined in terms of the nonlinear integral equation for the dressed energy (59) and the thermodynamic limit of the Bethe Ansatz equations (27). These can be expanded in 1/c without difficulty, and we obtain the following result for the particle density at order 1/c 2 where We recall that h is the chemical potential used to fix the particle number D. In order to derive (104) we used the following relations

Zero temperature ground state
Equation (60) for the ground state root density can be expanded in 1/c to yield Here 1 1 1 P is the indicator function, equal to 1 if the affirmation P is true and 0 if it is false.

The thermodynamic limit of correlation functions
In this section we perform explicitly the sum over intermediate states in (98) at order 1/c 2 in the infinite volume limit L → ∞.

One particle-hole excitations
Our starting point is C Λ 1 (x, t) as defined in (101). In the following we consider the different orders in the 1/c-expansion and derive integral representations of the corresponding contributions to C Λ 1 (x, t) in the thermodynamic limit. As we have noted before, we retain certain resummed expressions in this expansion for convenience, an example being the factor (1+ 2D c ). When we refer to a given order of the 1/c-expansion this should be understood modulo such factors.

Order c 0
Let us focus first on the leading order O(c 0 ), namely We rewrite this as Using (25) the sums over a and k can be turned into integrals over the root density ρ(λ), and the sum over n into an integral with density . Altogether we find where we used the expression (102) for the hole density ρ h at order c −2 .

Order c −1
We next turn to the c −1 term We rewrite this as This term involves either a sum over regularly spaced integers n that becomes an integral with density 1+2D/c 2π in the thermodynamic limit, or sums of the type (43) that can be expressed as principal part integrals over the root density. We obtain Introducing the Hilbert transformρ of ρ bỹ permits us to rewrite this contribution in the form

Order c −2 : first contribution
We now consider contributions involving the factor i =a which are more delicate. The first term on the right hand side gives rise to a contribution We rewrite this as The two terms are of the form (49) and we apply (57) to express them in terms of the root density ρ(λ) and the pair distribution function γ −2 (λ) defined in (33), with a triple integral with successive principal values defined in (54). This yields The definition of the successive principal value integral allows us to rewrite it in terms ofρ, to give

Order c −2 : second contribution
The second term on the right hand side of (116) is particularly cumbersome to deal with. We first treat the case i = j separately, and for all terms with i = j we apply a partial fraction decomposition with respect to n, so that we have only one n appearing in the denominator. Finally we split the sum over n as the difference of sums over vacancies and particles. Specifically, we have for f (u) = u 2 e ix u+it(λ 2 a −(λa+u) 2 )) i,j,n ∀k, λa,n =λ k |λa,n|<Λ i,j =a In all these terms, the conditions i, j = a only give rise to subleading contributions in L, so that they can be discarded. The first term on the right hand side of (121) gives rise to a contribution to C Λ 1 (x, t) of the form As this is proportional to L −4 and only involves three sums the dominant contribution arises from the double pole. Using n =0 1 n 2 = π 2 3 for the sum over n, we obtain

Order c −2 : third contribution
The second term on the right hand side of (121) gives rise to a contribution The sum is of the form (32) and according to (33) in the thermodynamic limit gives rise to integrals over the pair distribution function

Order c −2 : fourth contribution
The third and fourth terms in (121) are "hybrid" terms mixing sums over λ i 's and sums over regularly distributed n's. They give rise to a contribution to C Λ 1 (x, t) of the form By symmetrizing over i, j, one obtains a sole pole in n, but since n is regularly distributed and avoids only the pole one can convert the sum into a principal value integral. This leads to integrals with two successive principal values This can be simplified further by expressing the rightmost double integral in terms ofρ(λ).
To that end, let us consider the integral of this term as a function of µ with an arbitrary continuous function ϕ(µ). Using (53) we have Under the simultaneous principal value triple integral it is legitimate to decompose and split the integral into two since |µ − v| > : We then use (53) to express the two simultaneous principal value triple integrals in terms of successive principal value integrals The first integral on the right hand side is ϕ(µ)ρ(µ) 2 dµ while the third equals minus the left hand side. Using that this identity holds for any continuous function ϕ(µ) we conclude that Putting everything together we obtain

Order c −2 : fifth contribution
The fifth and fourth terms on the right hand side of (121) are of the form (49) and give rise to a contribution Applying (57) with successive principal values and then using (131) we find

Order c −2 ': sixth contribution
Finally, the last term in (116) gives rise to a contribution By again decomposing the sum over n as a sum over vacancies minus a sum over particles we find 5.1.9 Result for the contribution of one particle-hole excitations We leave the remaining contributions to C Λ 1 (x, t) untouched, i.e. in sum form, since they will be cancelled by contributions from two particle-hole excitations to the correlator. Our final result for C Λ 1 (x, t) is thus given by where we have defined 5.2 Two-particle-hole excitations

A partial fraction decomposition
The computation of C Λ 2 (x, t) defined in (101) is slightly different. In order to proceed we decompose the form factor into partial fractions with respect to n, and then m: We now carry out the sums over m and b in order to bring this to a form similar to the contribution from one particle hole excitations. We will denote the resulting five terms by Σ i for i = 1, . . . , 5 and treat them one at a time.

First term Σ 1
In this subsection we take the thermodynamic limit of We begin by splitting the exponential factor Performing the sum over m for the second term in (142) gives wherẽ (144) The advantage of this representation is that the pole in m is now only of order 1. Writing the sum over m as a sum over m = 0 minus sums over particles one obtains The first term is a sum over regularly spaced integers m with only a simple pole. In the thermodynamic limit it can therefore be expressed in terms of a principal value integral with a constant density 1+ 2D c 2π . The second term is of type (43) and gives rise to an integral over the root density ρ(λ) in the thermodynamic limit. The last term is negligible in L. We find where we defined and

Second term Σ 2
The next contribution is Writing the sum over m again as sums over vacancies minus particles we have In the first two lines of (150) the sums over n are regular. The first line involves only sums of the form (43), while the second line is of the form (49) and the thermodynamic limit can be worked out using (57). The third term, after splitting the sum over n as as sums over vacancies minus particles is seen to be negligible in L. Hence we obtain where B x,t (λ) is defined in terms of principal value integrals by

Third term Σ 3
In this subsection we take the thermodynamic limit of Expressing the sum over m as the difference of sums over vacancies and particles Σ 3 reduces to terms of the form (43) that can be readily expressed as integrals over root densities. We obtain where we have defined

Fourth term Σ 4
The next contribution is given by and can be treated in complete analogy with Σ 3 . We obtain where we have defined

Fifth term Σ 5
The final contribution to C Λ 2 (x, t) is In order to take the thermodynamic limit we rewrite this as The first two lines are of type (43) while the third and fifth lines are of types (49) and (32) respectively. Finally, in the fourth line we use that n =0 1 n 2 = π 2 3 to arrive at

Result for the contribution from two particle-hole excitations
Combining the results of sections 5.2.2-5.2.6 we arrive at the following expression for the two particle-hole contribution to the density-density correlation function where Ω Λ 2 has been defined in (148).

Compensation of divergent parts
As explained above, the O(c −2 ) contributions due to one-and two particle-hole excitations are individually divergent in the thermodynamic limit. The divergent parts are given in (138) and (148) respectively. Their difference is Crucially it vanishes for the class of root densities we use in our Λ-regularization discussed in Section 4.3.3, i.e. ρ(λ) = 0 for |λ| > Λ. We conclude that within our regularization scheme all divergences cancel at order O(c −2 ), but they do so in a non-trivial fashion: divergent contributions from intermediate states with one particle-hole excitation precisely cancel those arising from intermediate states with two particle-hole excitations.

Compensation of contributions that depend on the choice of representative state
As we have seen above, the contributions from both one-and two-particle-hole excitations in the thermodynamic limit individually depend on the choice of the representative state through the pair distribution function γ −2 (λ). Importantly, these contributions exactly cancel one another and the full correlation function does not depend on the representative state.

Λ-regularized correlation function
Combining the results for the one and two particle-hole excitations we obtain the following result for the dynamical density-density correlator in the Λ-regularization where the integrand is given by Here the contributions due to one and two particle-hole excitations are respectively The functionρ(λ) is defined in (114) and the four functions A x,t , B x,t (λ), C x,t and D x,t are given in (147), (152), (155) and (158) respectively. Some comments on the term 4π 2 c 2 (µ − λ) 2 ρ(λ)ρ(µ)ρ h (µ) 2 are in order. This term arises from the sum of the contributions involving the pair distribution function γ −2 (λ) in both C Λ 1 (x, t) and C Λ 2 (x, t). Strictly speaking it therefore involves two particle-hole excitations as well one particle-hole excitations. Since it does not involve double integrals, as is the case for the other contributions from C Λ 2 (x, t), we have chosen to include it entirely in χ (1) x,t (λ, µ). It can be interpreted as a "dressing" of contributions arising from one particle-hole excitations by two particle-hole excitations.

Dynamical correlations
The result (164) gives the thermodynamic limit of the Λ-regularized correlation function. We now remove the cutoff dependence by taking the limit Λ → ∞. The resulting ill-defined integrals (92) are to be understood as distributions following (96). To express the limit Λ → ∞ in terms of well-defined integrals we consider the expansion of f it follows from (96) that we can express the limit Λ → ∞ of (164) as a function of x and t = 0 For the energy of a macro state to be well-defined we need ρ(µ) = o(µ −2 ) for µ → ∞. From this we have ϕ (2) x,t (λ) = 4 Alternatively, one can also write, using (97)

Static correlations
The result (169) is singular for t → 0 since it behaves as 1 √ t e −x 2 /t . However, in a distribution sense we have 1 we have the following representation of the static correlator as a function of x Alternatively, one can also write 5.4 Dynamical structure factor in arbitrary macro states for all ω, q at order c −2 Given the correlation function (169) for all x and t we can determine the dynamical structure factor (DSF) S(q, ω) by taking the Fourier transform It is convenient to decompose S(q, ω) in terms of the contributions of one and two particle-hole excitations, which we denote by S (1) (q, ω) and S (2) (q, ω) respectively: In practice we determine the dynamical structure factor by first computing the Fourier transform (175) of the Λ-regularized correlator (164) and then taking the limit Λ → ∞, which turns out to be straightforward.

One particle-hole contributions to the dynamical structure factor
The contribution of the one particle-hole excitations to the DSF S (1) (q, ω) is obtained from the relation The Λ → ∞ limit of (177) is straightforward and yields at order O(c −2 ) where we have defined q . (179)

Two particle-hole contributions to the dynamical structure factor
The two particle-hole contributions involve the functions A x,t , B x,t (λ), C x,t and D x,t given in (147), (152), (155) and (158) respectively. Their simple dependence on x and t allows for a straightforward computation of their contribution to the DSF. For example, by first integrating over x, then over v, then over t and finally over u we find Here we have set The limit Λ → ∞ of this expression is again routine. It is however not immediately obvious that the double integral over λ and µ in (180) is well-defined, since one of the factors in the integrand exhibits a non-integrable singularity. A closer inspection reveals that this singularity is cancelled by the product of root densities. We will show below by means of a change of variable that the double integral is indeed well-defined. All other terms involving the functions B x,t (λ), B x,t (µ) and D x,t can be computed analogously. The term involving A x,t however requires a slightly modified approach, since following through the same steps as before would split the 1/v 2 term into a sum of two quantities that are individually divergent. In order to circumvent this problem we replace the 1/v 2 by 1 v 2 + 2 and send → 0 in the final result. We obtain Putting everything together we obtain the following result for the contribution of two particlehole excitations to the DSF In order to make the convergence of this integral explicit we perform a change of variables from λ, µ to and define In terms of the new variables we have The integral over p only has singularities that are integrable in a principal value sense, and after the integral over p has been carried out the integral over z only has a singularity that is integrable in a principal value sense. We conclude that (187) is well defined and can be straightforwardly evaluated numerically.

Numerical evaluation of the dynamical structure factor
In this section we numerically evaluate the integral representations (178), (187) in order to determine S(q, ω) for the two examples of root densities introduced in Section 3.4, namely thermal states and the non-equilibrium steady state after a quantum quench from the ground state at c = 0.

Zero temperature
We first consider the zero temperature case for density D = 0.404 and c = 3 in (106). The value D/c ≈ 0.13 is well within the expected range of validity of the 1/c-expansion. The same holds true for all other cases considered below. In Figure 3 we present numerical results for the DSF at order c −2 as well as for the one particle-hole contribution S (1) (q, ω). It is well known that at zero temperature the one particle-hole contribution to the DSF is non-zero only in a certain region of the q, ω plane for kinematic reasons, and exhibits (not necessarily divergent) singularities at the edges of its support [29,31,37,75]. We note that although the DSF is expected to diverge near the upper threshold, the divergence near the lower thresholds is a consequence of the 1/c expansion, that produces logarithms instead of a finite behaviour with a fractional c-dependent exponent. Comparing the full result (left panel) to S (1) (q, ω) (right panel) we observe that the contributions due to two particle-hole excitations significantly modify the numerical values of the DSF within this region. S (2) (q, ω) is also non-zero outside the region, but this effect is barely visible in the plot.

Finite temperature
We next turn to the DSF at finite temperatures. Figure 4 presents numerical results for the full DSF S(q, ω) at order c −2 for thermal states with β = 5, c = 4 and D = 0.396. For comparison we also plot the one particle-hole contribution S (1) (q, ω). Like in zero temperature case, for these parameter values the one particle-hole contribution already gives a fairly good account of the full DSF. The two-particle-hole contribution modifies some details that become increasingly significant for q > 2q F . The main difference to the zero temperature case is the emergence of spectral weight at negative frequencies and the "washing out" of the threshold singularities.
In Figure 5 we consider the DSF for a different thermal state characterized by a higher temperature β = 1 and D = 0.38, c = 4. The differences between the S(q, ω) and S (1) (q, ω) are difficult to discern in these plots. In order to get a more precise notion of the relative  (1) (q, ω). The color scale is the same for both plots.
contributions of S (1) (q, ω) and S (2) (q, ω) to the DSF for these values of D, c and β, we show a number of "constant momentum cuts", i.e. plots of S(q, ω) as a function of ω for fixed q, in Figs 6 and 7. Fig 6 gives representative results at "small" momenta, defined as q q F . We see that the contribution from two particle-hole excitations is negligibly small. This is in perfect agreement with observations made in Ref. [69] based on comparisons with numerical computations for a finite number of particles. Our results makes this observation fully quantitative in the thermodynamic limit in the framework of a 1/c-expansion.    We see that it grows with q and for the values shown is no longer negligible.

DSF in a non-equilibrium steady state
In Figure 8 we show numerical results for S(q, ω) and S (1) (q, ω) for the root density given in Section 3.4.2. The latter describes the stationary state reached for the interaction quench of Refs [95,97,98], where the system is initialized in the ground state at c = 0 and density D = 1/π and time-evolved with the Lieb-Liniger Hamiltonian at c = 4. We observe that the two particle-hole contributions lead to a slight narrowing of the DSF for q > q F .

Analysis of the result in limiting cases
In this section we report a detailed analysis of our results for the density-density correlation function (171), (166) and the dynamical structure factor (178), (187). The details of the derivations of the results in this section are reported in Appendix C.

Asymptotics of equal-time correlations at zero temperature
At zero temperature with the root density (106), we obtain the following asymptotic behaviour at large x and at order c −2 where and γ E is Euler's constant. This expression is the large x behaviour of the 1/c expansion of the correlation functions, hence one has c large first and then x large. Combining CFT/Luttinger liquid theory with exact results provides the following prediction for the correlations at large x at fixed c [25][26][27][28]37] σ(x, 0)σ(0, 0) = with K given in (107), and with A 1 a known constant [33]. If one wishes to compare this expression with (188), one is a priori faced with two problems: (i) If one expands (190) in powers of c −1 one has to take first x large and then c large, which is the reverse of (188). Hence comparing (188) and (190) entails to commute two limits. This commutation is possible if our expansion in c −1 (188) is uniform in space.
(ii) There could be corrections to (190) that are subleading in x at fixed c, but become of the same order as the dominant term once expanded in c −1 (i.e. give rise to log(x) terms). An example would be a term ∝ x −4K+2 . These corrections would be visible in (188), but not in (190).
In the case of density correlations in the Lieb-Liniger model it follows from the exact large x expansion at fixed c [37] that there are no subleading corrections with the property described in (ii). We thus expand (190) in powers of c −1 . Since K → 1 when c → ∞ the power-law x −2K becomes x −2 corrected by logarithms and we find The coefficient A 1 depends on c but as its representation is rather complicated [33] we left calculating its c −1 expansion for future work. We see that it agrees with (188) if we identify In particular the critical exponents are reproduced at order c −2 . This both provides a check of our formula, and shows that our 1/c expansion is uniform in space.

Dynamical correlations asymptotics at zero temperature
At zero temperature (with the root density (106)) we can evaluate the asymptotic behaviour of the dynamical correlation function at large x, t at fixed It is convenient to define and set We obtain with Ref. [37] derived the full asymptotic expansion at large x, t for any value of c at zero temperature. The c −1 expansion of this result at order c −2 (without expanding the prefactors) is in agreement with (196). In particular the critical exponents ν ± are reproduced at order c −2 . This both provides a check of our calculation and shows that our 1/c expansion is uniform in time as well, since the large x, t and large c limits commute.

Asymptotics of dynamical correlations for a generic root density and Generalized Hydrodynamics
For a generic continuous root density ρ in the large x, t regime at fixed α (193) we obtain the following asymptotic behaviour where we recall the definition of α in (194). The first line arises from one particle-hole contributions, while the second and third lines are two particle-hole contributions. If the root density is not continuous the leading term in 1/t is still correct, but the higher order corrections may change. GHD [83,84] makes predictions for the coefficient of the 1/t term in the density-density correlator for any value of c [87,88]. For the sake of completeness we summarize the 1/cexpansion of the GHD results in Appendix C.1.4. The leading term proportional to 1/|t| of (198) is in perfect agreement with the order c −2 expansion of the GHD results. To the best of our knowledge this constitutes the most non-trivial check to date of GHD predictions in an interacting integrable model. Importantly, we can assess the accuracy of the GHD approximation outside the asymptotic large space and time regime by comparing it to the full correlations at order c −2 . In Figure 9 we show our results for the real part and the modulus of σ(x, t)σ(0, 0) for two thermal states at c = 4 together with the GHD approximation. We see that at high temperature β = 1 the GHD approximation is surprisingly good even at short times. At lower temperatures β = 3 the correlation is still very well approximated by GHD, but is seen to display damped oscillations in the absolute value that arise from the imaginary part of the correlations that decay as t −2 and is not accounted for by GHD. Figure 9: Correlation function C(x, t) ≡ σ(x, t)σ(0, 0) in a thermal state for x = 2αt as a function of t at c = 4, for β = 1 and D = 0.38 (left) and β = 3 and D = 0.386 (right). The three curves depict the real part (red), the modulus (green) and the GHD approximation (dashed blue).
In Figure 10 we present the analogous comparison for a non-equilibrium steady state with root density (62). This root density is "less regular" than thermal densities in the sense that it has a narrower peak at zero. As a consequence, we expect higher Fourier-like corrections to the oscillatory integral, whereas GHD describes saddle-point-like corrections only. We indeed observe a more pronounced discrepancy for short or intermediate times, but the agreement at later times is still excellent, and globally remains very good. Figure 10: Correlation function C(x, t) ≡ σ(x, t)σ(0, 0) in the non-equilibrium steady state (62) for x = 2αt as a function of t at c = 4, for D = 1/π (left) and D = 1/(2π) (right). The three curves depict the real part (red), the modulus (green) and the GHD approximation (dashed blue).

A simplified expression at zero temperature
Equations (178), (187) for the DSF can be simplified at zero temperature. The one particlehole contribution can be written as where we have defined The contribution from two particle-hole excitations can be simplified by carrying out the integral over p in (187) Here we have defined and
Here C 0,1,2,3,4 are c-dependent constants and the exponents µ ± have simple expressions that depend on c. The dots encompass less singular pieces |δω| µ with a c-dependent exponent µ > µ ± , and regular pieces C 5 δω + O((δω) 2 ). In the framework of the 1/c expansion these power-laws give rise to logarithms These expansions are valid if we take ω close to ω ± first and then consider the large-c limit, and in order to compare with our result we have to commute these two limits. Importantly, the less singular pieces that are subleading in δω can also produce logarithms if their exponent goes to 0 when c → ∞. However, it follows from our asymptotic analysis in real space that there are no such terms. Comparing (206) with (204) and (205) we find that our result is in agreement with the non-linear Luttinger liquid predictions if we identify In particular we obtain the correct exponents at order c −2 . This provides a check of our result for the DSF and shows that it is uniform in q and ω.

Sum rule at zero temperature
The f-sum rule for the dynamical structure factor in equilibrium states reads [107] ∞ −∞ S(q, ω)ωdω = 2πDq 2 .
In our calculation, this sum rule has to be perfectly satisfied at order c −2 . It is a stringent test of validity of our formula since it has to be satisfied for all q and encompasses every single piece of the DSF. At zero temperature we obtain from Equation (199) that This means that the two-particle-hole DSF (201)S (2) We computed this integral numerically from (201) for several values of q. We find that (211) is indeed satisfied within the numerical accuracy of our calculation. The relative deviations of our results from (211) are around 10 −4 , which is quite satisfactory.

Detailed balance for thermal states
The dynamical structure factor of a thermal state at inverse temperature β should satisfy the detailed balance relation for all values of q, ω S(q, −ω) = e −βω S(q, ω) .
In our calculation the detailed balance relation for S(q, ω) should be perfectly satisfied at order c −2 . We note it is a very stringent test of validity of our formulas for S(q, ω), given that a thermal state at finite temperature corresponds to a generic root density with a complicated c dependence, while for an arbitrary root density there is no particular relation between S(q, −ω) and S(q, ω).
In order to check that our formulas for the DSF satisfy detailed balance at order c −2 , we need to evaluate (212) with ρ(λ) given at order c −2 in (103). We found convenient to definẽ i.e. the one-particle-hole DSF without the dressed piece coming from two particle-hole excitations. We recall that ω and q were previously defined in (179). It is straightforward to check numerically thatS (1) (q, ω) satisfies detailed balance at order c −2 Hence the following quantitỹ evaluated at c = ∞ should also independently satisfy detailed balancẽ We find that this indeed holds within the accuracy of our numerical computation, i.e. within a relative error of 10 −5 . This is quite satisfactory.

Behaviour at small q, ω
At small q, ω with fixed we find the following behaviour of the DSF We have set Here the term proportional to 1/|q| arises only from the one particle-hole contribution, while the constant term is due to two particle-hole excitations. This result can again be compared to GHD predictions, which at order c −2 give [83,84,88] S GHD (q, ω) = 2π 2 1 + 2D This is indeed in agreement with the leading term in (219) at small q. It would be interesting to see whether the subleading terms in (219) can be obtained by considering corrections to GHD following Ref. [108].

High frequency tail
Finally we consider the large-ω behaviour of the DSF S(q, ω) (at fixed q) in an arbitrary eigenstate |λ λ λ with a root density ρ(λ) that decays faster than any power law |λ| −n at infinity. For such states we find where ε = u 2 ρ(u)du is the energy of the state and δ its momentum defined in (89). The result (222) arises entirely from the two particle-hole contribution since the one particle-hole contribution decays faster than any power in ω for ω → ∞. The corrections to this leading behaviour can all be computed and expressed as a series in ω −1/2 . For example the next term is For eigenstates |λ λ λ corresponding to root densities that instead decay like a power law at infinity it is straightforward to see that the one particle-hole contribution to the DSF decays at large ω with the same power-law. For such root densities the large-ω expansion of the two particle-hole contribution breaks down at some order because the coefficients would diverge.
It has been shown some time ago that the large ω behaviour of the DSF S(q, ω) in equilibrium states is universal, with a decay ∝ q 4 ω −7/2 for quantum fluids with short-range interactions [109][110][111]. This behaviour was also observed to be in good agreement with scattering experiments [110]. Our result (222) is in perfect agreement with these findings, which again confirms that our 1/c expansion is indeed uniform in q, ω.

Conclusions
In this work we have introduced and developed an ab initio expansion of dynamical densitydensity correlation functions in the Lieb-Liniger model that can be performed within any energy eigenstate. It is a combined expansion in 1/c and in the number of particle-hole excitations taken into account in the spectral representation of the dynamical correlation function. The expansion has a well-defined thermodynamic limit and is uniform in all x and t, or equivalently all ω and q. We have obtained fully explicit and readily usable expressions for both the correlator and the dynamical structure factor at order O(c −2 ) which take into account all one-and two particle-hole excitations, Equations (171), (166), (178) and (187).
The main obstacle we faced in deriving these results occurs at order O(c −2 ). Indeed, the leading O(c 0 ) term of the expansion is simply the result for impenetrable bosons, which can be straightforwardly obtained using the mapping to free fermions [2,112]. In terms of the form factor expansion the only non-zero form factors are those involving a single particle-hole excitation, and they are all equal. The O(c −1 ) term is almost as simple since its form factor expansion is identical to the impenetrable limit case albeit with a root density dependent numerical modification of the form factors. In contrast the O(c −2 ) contribution comes with a number of complications.
As is well-known the form factor expansion generally exhibits non-integrable singularities whenever two rapidities coincide. In the framework of the 1/c expansion these first arise at order O(c −2 ) for contributions involving both one-and two particle-hole excitations. The presence of such singularities precludes directly taking the thermodynamic limit and expressing the spectral sum as integrals over root densities in a simple way. Indeed, we find that the contributions from both one-and two particle-hole excitations are individually divergent in the thermodynamic limit, but their sum is not. Even after compensating the divergent parts they individually depend on the particular choice of representative state and cannot be expressed in terms of the root densities. But remarkably, and reassuringly, their sum -and thus the correlation function -is representative-state-independent, i.e. depends only on the root density. These cancellations eventually leave a piece that can be interpreted as a dressing of the contribution due to one particle-hole excitations by two particle-hole excitations. Although this vanishes for the zero-temperature ground state as well as for any zero-entropy states it is non-zero in general and is crucial for detailed balance to be satisfied in thermal states. Such a fine-tuned "regularisation" of the divergences could only be achieved with a careful treatment of the thermodynamic limit of the exact spectral sum in a finite volume. Anticipating that for other quantities the representative-state-dependent parts may not al-ways compensate one another we derived a formula for their average over all representative states for a given root density.
We have verified that our results are in full accord with known results including CFT and (non-linear) Luttinger liquid theory predictions for zero-temperature critical exponents, thresholds singularities, sum rules, detailed balance relations and high frequency behaviour. We have also recovered the order O(c −2 ) GHD predictions for Euler scale density correlations in finite entropy states. This constitutes the most non-trivial verification of GHD in an interacting integrable model. We have also determined corrections to GHD and compared the GHD result to the full correlator at order O(c −2 ) outside the asymptotic regime. We found that GHD provides a rather good description of the correlator even at short times and distances.
The framework developed in this work is not restricted to density correlations in the Lieb-Liniger model but is expected to apply to any local operator in any integrable model that has a well-behaved expansion around a strong coupling limit. One example is the large anisotropy regime of the spin-1/2 Heisenberg XXZ chain [113,114]. A significant complication that occurs in that case is the presence of string solutions to the Bethe Ansatz equations. The restriction to local operators is crucial as the spectral representation of two-point functions of semilocal operators such as the field ψ(x) are dominated by a completely different set of excited states [76] and does not allow for an expansion in the number of particle-hole-excitations.
Our work opens up several interesting lines of further enquiry. First, our analysis should be extended to higher orders in the expansion. The O(c −3 ) term still involves at most two particle-hole excitations, but the expansions of the Bethe equations and the determinant in the expression for the form factors become more involved. Second, the repulsive Lieb-Liniger model is particularly simple in that the Bethe equations have only real roots. It would be very interesting to extend our analysis to a model with complex roots, e.g. the spin-1/2 Heisenberg XXZ chain. Third, our framework is readily generalized to quench dynamics [115] by combining it with the quench action approach [89,90]. Here the novel feature is that the spectral sum involves "overlaps" that multiply the form factors. Finally, it would be interesting to recover results obtained from the 1/c-expansion considering corrections to GHD as well as the thermodynamic bootstrap program [116].
Hence successive principal value triple integrals can be written as If G(λ, µ, ν) is a function without singularities, then we have where and i, j, k range e.g. between −L 2 and L 2 . In (226) one has the freedom to exclude some values, e.g. consider i = j, since this only amounts to subleading corrections in L that vanish when taking the limit. The integrand of (225) is of this type. Hence one can write Separating the four sums and changing variables so that the argument of f is always Finally we turn this into a simultaneous principal value integral by adding the condition i = k and obtain Equation (53).

A.2 Proof of Equation (56)
We note that formulae (55) are direct consequences of Equation (56), as can be seen by interchanging the dummy variables. We start with representation (225). As the integrand is regular one can impose that |λ − µ| > and |µ − ν| > with an error O( ) + O( ). This allows one to separate the integral into four pieces and make appropriate changes of variables so that the argument of F is always λ , µ , ν . One sees that in the four cases one has |λ − µ | > and |µ − ν | > . Hence which is precisely (56).

B.1 Reduction to a combinatorial problem
For a given solution to the Bethe equations {λ i } i ∈ S L we define the set of pairs of rapidities that belong to the same bin We have Let us show that when the pairs of rapidities are not in B, the sum is negligible. We observe that 1 L 3 provided that |k − p| > 1, i.e. if the bins to which λ i and λ j belong are not adjacent. Here C 0 = min y [ρ(y)] −1 is a constant independent of the representative state and of the bins. Indeed, in this case we have |λ i − λ j | > (|k − p| − 1)C 0 L and there are (L L ) 2 pairs of roots.
Since there are D/ L bins, by summing over p and k these contributions are O( 1 L L ), and since L L → ∞, they are negligible in the thermodynamic limit.
If the bins are adjacent we have with C 1 another constant independent of the representative state and of the bins. Since there are D/ L bins and L L → ∞, these contributions are also negligible in the thermodynamic limit. Hence we have with the o(L 0 ) being independent of the representative state. Hence we also have we have (240) To go further and decouple the average over the representative states one needs to ensure that the modification of rapidities in one bin does not notably affect the distance between rapidities in another bin. From the Bethe equations, a modification of order L of DL rapidities modifies the distance between two rapidities i, j in the same bin by an order 1 L DL L (λ i − λ j ), which is indeed subleading compared to λ i − λ j . Hence one can assume at leading order in L that the rapidities decouple, and given a bin k, sum over the rapidities of the other bins without modifying the values of the rapidities inside bin k. Thus one can write with S k L ⊂ S L the subset of S L containing states whose rapidities outside the bin [x L,k , x L,k+1 ] are fixed to those of an arbitrary representative state, and is the number of vacancies in [x L,i , x L,i+1 ]. Since around λ two consecutive vacancies are separated by 1 L(ρ(x L,k )+ρ h (x L,k )) , for λ i and λ j in the same bin [x L,k , x L,k+1 ]. This yields This reduces the problem to evaluating the large K, M limit at fixed K/M of the following combinatorial quantity (244)

B.2 The generating functions
To simplify the expression of C M,K , we would like to recast the sum over pairs of integers into a sum over (next-nearest-)neighbouring integers. We exactly rewrite C M,K in the form Introducing with a 0 = 0, we have Let us now determine the asymptotic behaviour of C

[m]
M,K for large K, M at fixed K/M . Summing separately over a 1 , ..., a m , one obtains the following recurrence relation where we use conventions such that C this recurrence relation implies that Expressing S [m] (x, y) as we obtain the following generating function

B.3 Asymptotics of the coefficients
We now use Ref. [82] which shows how to determine the asymptotic behaviour of combinatorial coefficients from the analytic behaviour of their generating function 2 . One obtains This implies that and substituting this into (247) we obtain at leading order Using the asymptotics (255) one then finds in the limit M, K → ∞ at fixed K/M 2 Specifically, in order to have only a simple pole in the generating function as in [82], we definē We then integrate by parts

B.4 Conclusion
Coming back to (243), we have when L → ∞ which yields In the limit L → ∞ we then arrive at the result (41) lim C Derivations of the results presented in Section 7 C.1 Correlation functions C.1.1 Asymptotics of static correlators at zero-temperature The study of the asymptotic behaviour of (173) at large x at zero temperature reduces to the asymptotics of the Fourier transform of a given function f (u). These asymptotics depend on the regularity of the integrand, hence at leading order on points of non-analyticity of f on the real axis. We have the following behaviours.
• If f (u) has a discontinuity ∆ = f (u + 0 ) − f (u − 0 ) at u 0 and is otherwise regular, then for x → ∞f This is straightforwardly obtained with an integration by part.
• If f (u) ∼ ∆ log n |u − u 0 | for u > u 0 and is regular and bounded for u < u 0 , then Here the constants p 1,2 are given by where γ E is Euler's gamma constant. If f (u) ∼ ∆ log n |u − u 0 | for u < u 0 and is regular and bounded for u > u 0 , then the result is multiplied by −1 and p 1 , p 2 are changed to their complex conjugates p * 1 , p * 2 . These relations are obtained from the relation expanded around α = 0.
• If f (u) ∼ ∆ log n |u − u 0 | for both u < u 0 and u > u 0 , then we havê These equations directly follow from the previous results.
In order to determine the large x behaviour of the correlation functions, we also need zero-temperature result forρ(x) defined in (114) and the large x behaviour of the functions A x,0 , C x,0 and D x,0 defined in (147), (155) and (158) respectively The asymptotics of C x,0 and D x,0 follow from (298) and (299) for a generic root density. As for A x,0 , integrating (147) by parts we obtain for a generic root density Specializing to zero temperature at leading order in c −1 , it yields that is As for the B x,0 (λ) and B x,0 (µ) terms, they require a special treatment since they cannot be decoupled from the λ, µ integrals. The B x,0 (λ) term involves the following functions for n = 0, 1, 2, whose we wish to determine the asymptotic behaviour at large x, by computing its Fourier transformf n (q) = ∞ −∞ e −iqx f (x)dx. We have (at leading order in c −1 ) Specializing this relation to the ground state root density we obtain We note that the non-integrable divergence near λ = q is compensated by the argument of the log going to 1 in this limit. In the vicinity of q = Q we have for η > 0 where the last integral is Similarly we find andf n (q) does not have discontinuities elsewhere. This implies that Then one builds the full B x,0 (λ) term from the functions (273), in particular by taking into account the remaining oscillatory µ integral. One obtains that the B x,0 (λ) term gives contributions that decay as cos(2q F x)/x 2 and of order c −2 , which are encapsulated in the result in the O(c −2 ) term of the A in (189). A similar analysis shows that the B x,0 (µ) term also gives contributions that decay as cos(2q F x)/x 2 . From these various relations, it is straightforward albeit tedious to determine the asymptotics of the static correlation functions. Putting everything together we find that χ (1,2) x,0 (λ, µ) given in (166) contributes to the large-x behaviour of the density-density correlator as follows: • O(c 0 ) contribution of χ (1) x,0 (λ, µ) • O(c −1 ) contribution of χ (1) x,0 (λ, µ) • O(c −2 ) contribution of χ (1) x,0 (λ, µ) • Contribution of χ (2) x,0 (λ, µ) This establishes (188).

C.1.2 Asymptotics of dynamical correlations zero temperature
The study of the asymptotic behaviour of (164) at large x, t at fixed α = x 2t at zero temperature reduces to the study of an oscillatory integral of the type In this regime, the integral is dominated by the point where the phase has an extremum as a function of u, which is α defined in (194). If f is regular and α in the support of f , then we have I(x, t) = If f has singular points one has to combine (285) with the results of Section C.1.1. The correlation function (164) is expressed as a double integral, one over λ with a factor ρ(λ) and one over µ with a factor ρ h (µ). Because of the very particular structure of ρ(λ) at zero temperature (106), the saddle point α necessarily lies within the support of either ρ(λ) or ρ h (µ), but not both. Hence if |α| > q F , the λ integral is dominated by boundary effects as in the static case, while the µ integral is dominated by the saddle point. If |α| < q F , the converse holds true.
Let us detail the case |α| > q F (the case |α| < q F is very similar). We perform a change of variables λ → λ + α and µ → µ + α in (164) in order to move the saddle point to 0, which results in shifting the arguments of the root densities by α . The µ integral is then simply evaluated at µ = 0, while the λ integral is dominated by the vicinities of the points Q − α and −Q − α . Using the results for (262) with x = −2t(Q − α ) we obtain the leading contribution from Q − α to the integral over χ (1) x,t (λ, µ), with ± = sgn (t)

C.1.3 Euler scale asymptotic behaviour
In this section we will assume ρ to be continuous. If it is not continuous the leading 1/t behaviour is unchanged, but the 1/t 2 corrections might differ. For a generic continuous root density at large x, t and fixed α = x 2t , the two integrals over λ and µ in the correlation function (164) are both dominated by the saddle point at α . Applying (285) to the one particle-hole contribution gives The contribution due to two particle-hole excitations is more subtle and requires determining the asymptotic behaviour of oscillatory integrals with principal values, whose saddle point falls on the singularity. The general strategy is to write each singularity as and then to carry out a regular asymptotic analysis of the multiple oscillatory integrals successively. The sgn (ξ) factors introduce discontinuities which result in contributions on top of those from the saddle points. Let us treat the case of C x,t in detail. We write and then apply a saddle point approximation to the u and v integrals using (285) to obtain This can be simplified by performing an integration by parts on the ξ integral of the subleading term Similarly one finds for the D x,t term D 2α t,t = iπ sgn (t) ∞ −∞ ρ(ξ)ρ h (ξ)ξ sgn (α − ξ)dξ To deal with the A x,t term we use that in a distributional sense Putting everything together we arrive at (198).

C.1.4 GHD predictions
The GHD result for the asymptotics of the density-density correlator is [87,88] σ(x, t)σ(0, 0) = with ρ, ϑ defined in (27) and (28), and where the other functions are defined as follows: we find for z ≈ 0 We observe that for z close to zero we have Z − < Z + if and only if z < − η 2q (q +2Q) , in which case Z + − Z − = q |z| 1−z . Above the threshold we have 1 1 1 ω − <ω<ω + = 0 so that the last term in (201) vanishes. Of the remaining terms only the one proportional to Z + −Z − z 2 diverges near z = 0, so that The behaviour near the lower threshold is obtained through a similar analysis.

C.2.2 Behaviour at small q, ω
We start by writing the two particle-hole contribution as where Ψ 1,2 denote the first and second terms respectively. The integral for Ψ 1 with q, ω → 0 at fixed γ = ω 2q is well-defined and finite. In this limit we have q 3 = q 4 = γ + p(1 − z) and q 1 = q 2 = γ − pz. Changing variables to v = γ − pz and u = γ + p(1 − z) we have As for Ψ 2 , we first perform a change of variables from µ to v = q + λ − µ We now observe that the four ρ factors are invariant under the change of variable v = −q 2γ−2λ−q 2v . (324) We apply this change of variable to all the terms except and express the term as one half of itself plus one half of itself after the change of variables. We obtain for q > 0 Since there are no non-integrable divergences in the integrand at small v, in this representation one can set q = 0 in the ρ terms as well as in the integrand, at small q . It yields (329) We obtain then the claimed result.

C.2.3 High frequency tail
We start with the representation (184) for the two particle-hole contribution to the DSF expressed as a single double integral. We first decompose the double integral into the two regions |q + λ − µ| > and |q + λ − µ| < and focus on the latter part. Since we have assumed that ρ decays faster than any power law at infinity λ has to remain smaller than any power law of ω for the integral not to vanish at any order O(ω −n ). Since |q + λ − µ| < the same holds true for µ. But this implies that λ necessarily grows as ω 1/2 , which makes this contribution vanish at any order O(ω −n ). Hence at any given order O(ω −n ) we can impose |q + λ − µ| > . This removes all poles in the integrand of (184) and one can consider all contributions separately. Moreover, since ρ decays faster than a power law at infinity the term proportional to ρ( ω −q 2 2q )ρ h ( ω +q 2 2q ) is negligible at order O(ω −n ). We then split the µ integral into the sum of positive and negative µ parts and perform the change of variables This way the DSF can be brought to the form dz ρ(λ)ρ z + f + (z, λ, ω ) g + (z, λ, ω ) where the dots indicate subleading corrections that decay faster than any inverse power in ω and f ± (z, λ, ω) = ± 1 2 We now observe that any part of the integral where the argument of one of the two ρ's grows as a power-law in ω will give contributions that decay faster than any power-law, since ρ is assumed to decay faster than any power-law at infinity. From the expression of f ± one sees that z cannot grow faster than ω 1/4 . Consequently, with an error that goes to zero faster than any power law in ω one can replace the limits of the integrals and the arguments of the ρ h 's by ±∞. This gives We now expand f ± (z, λ, ω ) and g ± (z, λ, ω ) in Laurent series in λ, λ−z and ω −1/2 , and Taylor expand ρ(z + f ± (z, λ, ω )). This produces terms of the type ρ(λ)ρ (a) (z)λ b (λ − z) d ω −e/2 with a, b, e ≥ 0 integers and d a positive or negative integer. We integrate this by parts a times over z so that the integrand involves only ρ(z), and then write the full result S (2) (q, ω) as one half of itself plus one half of itself after swapping the dummy variables λ and z. We observe that there remain only positive powers d ≥ 0, and one obtains the first two terms of the expansion ω −7/2 , ω −9/2 stated in the text.