Mobile impurities in integrable models

We use a mobile impurity or depleton model to study elementary excitations in one-dimensional integrable systems. For Lieb-Liniger and bosonic Yang-Gaudin models we express two phenomenological parameters characterising renormalised inter- actions of mobile impurities with superfluid background: the number of depleted particles, $N$ and the superfluid phase drop $\pi J$ in terms of the corresponding Bethe Ansatz solution and demonstrate, in the leading order, the absence of two-phonon scattering resulting in vanishing rates of inelastic processes such as viscosity experienced by the mobile impurities


Introduction
Currently there is an ongoing interest in physics beyond Luttinger Liquid in one-dimensional systems [1]. An example of such behaviour is provided by systems featuring mobile impurities, representing one or few particles moving in a correlated background liquid. This situation was recently realised in ultra-cold atoms by either using atoms in different hyperfine states [2,3,4] or by creating a highly imbalanced mixture of different atoms [5,6]. By exploiting unprecedented control over parameters of the underlying Hamiltonian these experiments have uncovered rich dynamics of mobile impurities and have revived interest in their theoretical studies either in the context of liquid helium [7,8,9] or, more recently, in the context of one-dimensional quantum liquids [10,11].
There are two main types of questions one can ask about these systems. The first one concerns dissipative dynamics of mobile impurities resulting from the strong coupling to environment [12,13,14,15,16,17,18,19,20], while the second type of questions deals with correlation properties of quantum liquids affected by the presence of mobile impurities [21,22,23]. These questions can be extended to one-component quantum liquids as their high-energy nonlinear excitations, like solitons, can be described as effective mobile impurities [24,25].
Theoretically, the physics of mobile impurities is modelled by a single localised degree of freedom coupled to an extended environment. For sufficiently low temperatures and weak external fields the latter is only slightly perturbed so it can conveniently be represented by a collection of phonons in the form of Luttinger Liquid (LL) [26,27,28,29,30]. The coupling between the localised impurity and LL can be characterised by only two parameters: the number of background particles, N , expelled from the vicinity of impurity and the phase drop, πJ, of the local superfluid order parameter. The existence of these parameters, called collective charges below, is attributed to the existence of two local conservation laws: the conservation of number of particles and conservation of momentum. The length scale separation between localised impurities and extended phonons allows for extraction of the values of collective charges from the response of the equilibrium dispersion relation ε(k) of impuritylike excitations to small changes in background density and velocity fields and leads to a phenomenological depleton model [14,16] for which ε(k) is the only physical input.
This model was successful in predicting the rate of momentum exchange between the moving impurity and the phononic bath giving rise to viscous friction experiences by impurities [10,11,13]. A slightly different, but equivalent approach was used for determining relaxation rates toward thermal equilibrium in one dimensional systems due to mobile impurities [17,31].
On the other hand, the contribution of mobile impurities to correlation functions leads to power law singularities near the excitation threshold ε(k) [1]. It can therefore be expected that the corresponding edge exponents are determined by the same collective charges characterising the impurity-background interactions in the depleton model. The relation can indeed be established using the method of Kamenev and Glazman [23] based on Ref. [9] for extracting the edge exponents from the magnon dispersion of a general Galilei invariant spin liquid. The method was further developed in Ref. [32] for a variety of Galilei invariant bosonic or fermionic models. A similar calculation based on the equilibrium dispersion ε(k), as we show below, leads to the collective charges of Ref. [14], thus establishing their connection to the correlation properties of one-dimensional quantum liquids with mobile impurities.
For a general interacting system, the dispersion relation ε(k) which determines both dynamics and correlations in systems with mobile impurities can be either calculated analytically using an appropriate perturbation scheme or obtained numerically. There exist, nonetheless, a class of exactly solvable one-dimensional systems, where the dispersion of excitations is known exactly in terms of Bethe Ansatz solution [33]. Important examples of such systems include Lieb-Liniger model [34] of bosons interacting via a delta function potential and the closely related bosonic Yang-Gaudin model [35,36] of interacting spin one-half bosons.
The goal of this paper is to treat excitations of Lieb-Liniger and bosonic Yang-Gaudin models as mobile impurities and find their collective charges in terms of the corresponding Bethe Ansatz solutions. We were able to express collective charges directly in terms of shift functions (otherwise known as dressed phases in BA literature) of elementary excitations. This result reproduces the conjectured identity of chiral linear combinations of the collective charges, so called chiral phase shifts, with BA shift functions [22,32,37,38]. Recently, the calculation of correlation functions in Lieb-Liniger model by the form-factor approach by Kitanine et al. [39] confirmed the expression of the edge exponents in terms of BA shift functions obtained in [38].
The expression of the collective charges in terms of BA shift functions allows us to find a very useful and transparent representation for N and J in terms of derivatives of the excitation's momentum with respect to BA parameters, see Eq. (32) and (41) below. These identities lead straightforwardly to exact vanishing of the phononic back scattering amplitude in the leading, two-phonon, order. This result should not come as a surprise given the exact integrability of the models, but it is far from obvious in the phenomenological approach based on mobile impurities. Up to now it could only be verified in a number of particular cases of weak and strong interactions [13,25].
The paper is organised as follows. In Section 2 we review the depleton model containing collective charges N and J and give their expression in terms of the dispersion ε(k). We demonstrate the relation between the collective charges and chiral phase shifts and use it to reinterpret various edge exponents in terms of N and J. In Section 3 we review BA solution for elementary excitations of Lieb-Liniger and bosonic Yang-Gaudin models. Then we obtain exact relations between the collective charges N and J and solutions of BA equations for Lieb-Liniger model in Section 4 and for bosonic Yang-Gaudin model in Section 4.2. Finally, we discuss in Section 5 phononic back scattering amplitude and show that it vanishes identically for these models. We present conclusions and open questions in Section 6. Technical details needed for the proofs are delegated into four Appendices for clarity.

Depleton model, its collective charges and edge exponents of correlation functions near excitation thresholds
Depleton model is designed to describe mobile impurities within path intgral formalism and is defined by the action S = S ph + S d , where describes extended environment by using the slow hydrodynamic fields ϕ(x, t) and ϑ(x, t). The latter is related to density as n(x, t) = n + ρ(x, t) = n + ∂ x ϑ/π, where n is the equilibrium density. Correspondingly, the field ϕ(x, t) describes the superfluid phase and is canonically conjugated to ρ(x, t). The Luttinger Liquid action, Eq. (1), describes phonons with linear dispersion relation ω = v s |k|, where v s is speed of sound. The relative magnitude of fluctuations of the fields ϕ and ϑ are controlled by the Luttinger parameter K = πn/mv s , where n is the average density and m is the mass of the background atoms.
To describe a localised object with coordinate X(t) and the canonically conjugate momentum P (t) we consider action This action describes a sharp "kink" in the smooth configuration of the hydrodynamic fields located at X(t) and parametrised by collective charges N (t) and J(t). The depleton charge N (t) represents the number of background particles expelled from the vicinity of the depleton, therefore πN (t) is the magnitude of the "kink" in the otherwise smooth configuration of ϑ.
Similarly, πJ(t) is the size of the "kink" in ϕ. In the absence of phonons, the charges take their equilibrium values N (t) = N , J(t) = J obtained from the condition ∂H/∂N = ∂H/∂J = 0. The dispersion relation is the equilibrium value of the energy H(P, N, J) = ε(P ). It was shown in [14] that the equilibrium collective charges are given directly in terms of the dispersion 1 , Here v(k) = ∂ε(k)/∂k is the velocity of the excitation. In addition to describing dissipative dynamics of mobile impurities [14,16] the depleton action S can be used to calculate leading power law behaviour of various correlation functions the vicinity of the excitation threshold ω ∼ ε(k) as shown in Appendix A. In this case the path integral is dominated by the stationary value of the action S which becomes logarithmically large and leads to the power-law behaviour ∼ |ω − ε(k)| −µ , where the edge exponent is determined by the equilibrium values of the collective charges (3), (4). 1 Ref. [14] uses notations N, Φ, where N has the same meaning as in this paper and Φ = πJ.
The collective charges N and J can be connected with the parameters δ ± , called chiral phase shifts, of the unitary transformation in the standard calculation [1] of the correlation functions using the X-ray edge approach [40,41,42]. In the latter, the chiral phase shifts characterise the discontinuities in the chiral phononic fields χ ± = ϑ/ √ K ± √ Kϕ created by boundary condition changing operators [43,44]. The charges N, J, on the other hand, characterise the discontinuities in the fields ϑ and ϕ. The simplest way to establish connection between these two sets of parameters is by comparing Eqs. (3), (4) with those of the phase shifts δ ± first obtained in Ref. [23,32] in terms of the mobile impurity dispersion relation. It is easy to see that the phase shifts and the collective charges are related by The expression (5) reproduces known results for edge exponents if rewrite it in terms of the chiral phase shifts as The first line reproduces exactly the results for the edge exponents µ 1,2 of Dynamic Structure Factor (DSF) of interacting bosons, Eq. (16) of Ref. [38], while the second line is Eq. (13) of Ref. [23] for the exponent µ m of Dynamic Spin Structure Factor (DSSF) of a general spin one half quantum liquid. To reproduce edge exponents µ ± of the spectral function A(k, ω) in Ref. [38] one notices that in this case, in addition to creating a depleton, one removes exactly one particle from the system, so one have to replace the number of missing particles N by N + 1. In terms of the phase shifts δ ± this leads to which reproduces exactly Eq. (17) of [38]. Similarly, to calculate µ ± one replaces N by N − 1, which corresponds to adding one particle so that, again in full accordance with Eq. (18) of Ref. [38]. The same logic applies for calculating the edge exponents away from the fundamental zone 0 < k < 2πn. Consider the lower edge, the Lieb II mode ε(k) = ε 2 (k) as an example. A remarkable fact with far fetching physical consequences is that ε 2 (q) is a periodic function, ε 2 (q+2πn) = ε 2 (q). This is due to the fact that the momentum difference 2πn results from the presence of additional uniform background supercurrent and the corresponding phase winding of 2π. The energies of the state with the supercurrent and without it are the same in the thermodynamic limit. The additional phase winding 2πl (l is an integer) affects the value of J and can be taken into account by the substitution J → J − 2l in Eq (5) so that the resulting edge exponent becomes which is the bosonic version of Eq. (87) of Ref. [1]. Thus the expression (5) obtained within the depleton formalism nicely combines previous results on edge exponents of dynamical correlation functions obtained using the mobile impurity model. The collective charges N and J entering this expression can be obtained from phenomenological dispersion of the mobile impurity excitations by Eqs. (3), (4). Below we obtain the collective charges for integrable models using excitation dispersion from their Bethe Ansatz solutions.

Bethe Ansatz solution and elementary excitations of Lieb-Liniger and bosonic Yang-Gaudin models
Up to now our discussion of the effective mobile impurity with dispersion relation ε(k) was purely phenomenological. We now turn to the models where the dispersion relation and, consequently, the parameters N, J determining the edge exponents can be obtained analytically by Bethe Ansatz. We consider Lieb-Liniger model [34] describing one-dimensional bosons interacting with a short-range potential V (x) = cδ(x) which is arguably the most studied continuous model solved by Bethe Ansatz. Generalising this model to two species of bosonic particles of the same mass and equal interaction couplings leads to bosonic Yang-Gaudin model [35,36]. This model can be conveniently formulated using effective spin 1/2 particles. Its ground state is fully polarised [45] and is identical to that of Lieb-Liniger model. Below we present main equations which allow us to obtain dispersion relation ε(k) used for calculation of the collective charges. All results of this Section can be found in Refs. [33,46] and are reproduced here to make presentation self-contained. The ground state of Lieb-Liniger model is characterised by the density of quasi-momenta ρ(ν) found from the equation Here K(λ) = ∂θ/∂λ = 2c/(c 2 + λ 2 ), where θ(λ) = 2 arctan(λ/c) is the scattering phase shift. We use units in which the mass of the particles is m = 1/2 and = 1. In these units the coupling constant c and quasi-momenta λ have dimensions of velocity. The "Fermi momentum" q limits the support of the ground state density and is determined from the normalisation condition The excitations of Lieb-Liniger model are in one to one correspondence with those of free fermions and consist of either a particle excited above the Fermi sea (Lieb I), so its quasimomentum λ is constrained by |λ| > q or a hole (Lieb II) with −q < λ < q [47]. The difference from free fermions is that the positions of the ground state quasi-momenta are shifted as a result of the excitation and these shifts contribute collectively to momentum and energy. This effect is characterised by the shift function F (ν|λ) which obeys the equation Taking the limit λ → ±∞ and comparing with Eq. (12) it follows that F (ν| ± ∞) = ∓πρ(ν). The shift function represents the relative change of the ground state quasi-momenta δνρ(ν) = ∓(F (ν|λ) + πρ(ν)) as a result of a particle/hole with quasi-momentum λ. Here and below the upper sign corresponds to a particle-like (Lieb I) excitation while the lower sign corresponds to a hole-like (Lieb II) excitation. Some useful properties of the shift function F (λ|ν) are summarised in Appendix B.
The momentum and energy of Lieb I,II excitations are given by The integrals in the right hand side represent the collective contribution of the displaced momenta in the "Fermi sea". Eliminating λ from Eqs. (15) and Eq. (16) leads to the dispersion relation ε(k). Apart from momentum, the dispersion relation ε(k) depends on the density n via the limiting momentum q fixed by normalisation (13). The chemical potential h is then fixed from the condition ε(q) = ε(−q) = 0. As both signs in Eq. (15) should provide the same result for λ = q, it is clear that k(q) = 0. Since for a hole the momentum increases with decreasing λ one can show that k(0) = πn and k(−q) = 2πn. For a particle k(−q) = −2πn.
In the bosonic Yang-Gaudin model in addition to Lieb I,II excitations there is another type of excitations corresponding to a spin flip of one of the particles. The flipped spin particle is introduced with quasi-momentum λ causing the change in the quasi-momenta δνρ(ν) =F (ν|λ) + πρ(ν). The shift functionF (ν|λ) is found from the equatioñ This equation differs from Eq. (14) by the factor two in the argument of the bare phase shift. The momentum of the magnon excitation is obtained by the shift of the ground-state quasi-momenta Eqs. (15), (17) The last equality follows from Eq. (19). The corresponding energy of a magnon can be similarly expressed as where ǫ is the solution of Eq. (18) and the last equality is obtained by integrating by parts.

Collective charges and shift functions in integrable models
The parametric expressions (15), (16) and (20), (21) of the previous section allows us to obtain the dispersion relations of excitation ε(k) in Lieb-Liniger andε(k) in bosonic Yang-Gaudin models in terms of shift functions F (ν|λ) andF (ν|λ) correspondingly. It is therefore expected that the same Bethe ansatz shift function determine the collective charges N, J or their chiral combinations δ ± (k). Indeed, the previous studies [38,22] suggested the following particularly simple and physically appealing relation for excitations of Lieb-Liniger model and a similar expression for magnon excitation of bosonic Yang-Gaudin model. The conjecture (22) is based either on identification of δ ± (k) with the parameters of a boundary changing operator of background particles [38] or on comparing the finite size spectra [48,49,37,22]. To the best of our knowledge, the direct proof of Eqs. (22) for integrable models is still lacking and the consistency of the thermodynamic definition Eq.

Lieb I and II excitations of Lieb-Liniger model
Our proof is based on the following expressions for the the derivatives of the dispersion relation with respect to the "natural" parameters λ and q. Defining for simplicity F ± = F (±q|λ) we show in Appendices C,D that and We also show there that and by differentiating both sides of Eq. (17) and comparing the result with Eq. (12) we have where we have used Eq. (58). We will also need the relation proven in Appendix D.1. We transform the derivatives (23), (24) of the dispersion relation with respect to the natural variables λ and q to those with respect to density n and momentum k. This can be achieved using Eqs. (25), (26) and (27). Re-introducing the particle's mass m and going back to ε = ±ǫ for particles/holes we rewrite Eqs. (23), (24) as If one identifies in these equations, accordingly to Eq. (22), the combinations and solves for δ ± using the upper sign for particles one obtains precisely Eqs. (6). The collective charges N, J can be obtained using Eq. (7) with the result For λ = ±q we have N = 2πρ(±q) = √ K = 1/J. In the case of a hole-like Lieb II excitation, one has k → −k, ε → −ε, but v(k) = ∂ε/∂k is unchanged. The expressions (3), (4) become Comparing Eqs. (28), (29) with Eqs. (25) and (26) allows one to rewrite the collective charges for both Lieb I (particles) and Lieb II (holes) via partial derivatives of momentum k(λ; q) as We are not aware of any previous studies discovering these nontrivial relations.

Magnon excitation in bosonic Yang-Gaudin model
For bosonic Yang-Gaudin model the phenomenological approach was used in Ref. [23] to express the phase shifts of the magnon excitation in terms of the magnon momentumk and derivatives of the magnon dispersion relationε(k) leading to the result (6) in which one replaces k →k, ε →ε and the group velocity v(k) →ṽ(k) = ∂kε. In later work [22] a relation between the phase shifts of the magnon excitation and Bethe Ansatz shift function, similar to Eq. (22) was proposed. In the notation of Ref. [23], and using a slightly different sign convention, this relation is given by To prove this relation by analogy to the case of the Lieb-Liniger model we need derivatives of the dispersion relation and derivatives of the momentum  (38) in full agreement with the results of Ref. [23]. The collective charges are as expected from Ref. [14]. Comparing Eqs. (30) and (39) we see that there are additional terms representing the bare momentum mṽ(k) of the impurity and its chemical potential πv s /K = mv 2 s /n = ∂µ/∂n. These additional terms cancel each other in Eq. (40) . In terms of the magnon shift functions, the collective charges are Again, comparing these expressions with Eqs. (36) and (37) allows to rewrite the collective charges via partial derivatives of momentum as

Phonon backscattering amplitude
It was shown in Ref. [14] that the viscous friction force acting on a moving impurity due to twophonon processes is proportional to the squared absolute value of the phonon backscattering amplitude Γ +− . It is expected that in integrable systems this amplitude vanishes due to the existence of infinitely many conservation laws. An example supporting this statement is provided by a dark soliton, which is a hole-like excitation of the Lieb-Liniger model in the weakly interacting regime and which was shown to have infinite life-time in Refs. [50,25]. Another example of a non-decaying excitation is a spin-flipped particle (magnon) in bosonic Yang-Gaudin model, as was shown in Ref. [13] in the limit of weak and strong interactions. Below we prove that to the leading two-phonon order the backscattering amplitude vanishes identically for Lieb-Liniger and bosonic Yang-Gaudin models for any value of interactions.
The expression for the back scattering amplitude was obtained in Ref. [14] as the following combination of partial derivatives of the collective charges Here M is the mass of the added impurity particle. As it stands, Eq. (42) is valid only for subsonic excitations. In our case they are the hole excitation of Lieb-Liniger model for which M = 0 since there is no additional particle and the magnon of bosonic Yang-Gaudin model for which M = m. For the supersonic particle-like excitation of Lieb-Liniger model the expression for the back scattering amplitude must be modified as explained below. We start with a hole-like excitation for which Γ +− is given by We convert the partial derivatives with respect to k and n into the partial derivatives with respect to the "natural" variables using where we have used Eqs. (28) and (29). Inverting the above expressions we get which can be used (with the lower sign) in Eq. (43). It leads to a particularly simply looking result as a direct consequence of the relations (32).

Summary and conclusions
To summarise, we have obtained the collective charges N and J of mobile impurities in integrable models directly in terms of the corresponding BA shift functions. Our method relies on the Bethe Ansatz solution which provides a unique parametrisation of the elementary excitations in these models by their rapidity as well as the parameter which limits the extent of the ground state rapidity distribution. This parametrisation is expressed via Bethe Ansatz shift functions and allows for exact calculation of derivatives of the excitation energy with respect to momentum and density which enter the phenomenological expressions for the scattering phase shifts and collective charges found in earlier works.
As a byproduct we have found a novel expression for the collective charges of the effective mobile impurity model in terms of the partial derivatives of the impurity momentum with respect to the Bethe Ansatz parameters mentioned above.
A straightforward consequence of these relations is the absence of phonon scattering off mobile impurities in the leading two-phonon order for all values of interaction parameters. This absence was previously conjectured based on approximate calculation in the limiting cases of weak and strong interactions. The proof of the expected absence of inelastic processes beyond two-phonon rate is an obvious extension of our work.
We expect that our methods can be generalised to study excitations in other models soluble by nested BA, such as, the fermionic Hubbard model and integrable spin chains. The lack of Galilean invariance in lattice models can be dealt with following techniques of Ref. [17].
A Semiclassical calculation of power-law edge exponents using depleton model We are interested in dynamical correlation functions of one-dimensional bosons. In particular we consider the zero-temperature Dynamic Structure Factor (DSF), and spectral function A(k, ω) = − 1 π Im G(k, ω) sgn ω, where is the Green's function. Here Ψ(x, t) and ρ(x, t) = Ψ † (x, t)Ψ(x, t) are boson annihilation and density operators, and T denotes time ordering. In the case where bosons have two internal states a, b we can define spin operators s + = s † − = Ψ † a (x, t)Ψ b (x, t) and the corresponding Dynamic Spin Structure factor (DSSF), It is well known (see Ref. [1]) that for various one-dimensional models the dynamical correlation functions have a rich structure in the (k, ω) plane. In particular, they exhibit power law singularities in vicinity of dispersion curve ε(k) of elementary excitations, This power law behaviour can be obtained by semiclassically evaluating the path integral with the action S = S ph +S d , see Eqs.
(1), (2). As we shall see below the semiclassical method is justified by the smallness of the energy difference ω − ε(k) leading to a logarithmically large action. Following Iordanskii and Pitaevskii [51] we divide the time countour into three intervals, (−∞, t 1 ), [t 1 , t 2 ], (t 2 , +∞). The kinematic constraints dictate that in this vicinity of the excitation energy ε(k) it is enough to consider configurations consisting of only one depleton [14,16] propagating in the time interval [t 1 , t 2 ]. Then the stationary configuration of the fields dominating the functional integral can now be described as follows: the depleton appears at t 1 at x 1 , propagates as a point-like particle with a constant velocity V = (x 2 − x 1 )/(t 2 − t 1 ) and disappears at t 2 at x 2 . Its trajectory is given by X(t) = x 1 + V t. This is reflected in the following behaviour of the collective chargeṡ Here the values N and J correspond to their corresponding equilibrium values calculated from Eqs. (3), (4) in which the momentum k is such that v(k) = V . The corresponding stationary configuration of the phononic variables is determined by solving the wave equation with the source terms Solving these equations leads to logarithmic behaviour of the action where t 21 = t 2 − t 1 , x 21 = x 2 − x 1 and we have introduced symmetric and antisymmetric combinations of the depleton collective charges.
To proceed with the correlation function we need the Fourier transform of the type where an infinitesimal imaginary part was added to avoid power-law singularities. The action of depleton S d (x, t) is considered as function of time and coordinate, and its full differential obeys The integral over coordinate x is performed by stationary phase method, which locks the momentum P of the depleton to the externally imposed value P = ∂S d /∂x = k and energy H = ε(k). The velocity of the depleton is now a function of momentum and its stationary trajectory is x = v(k)t. The collective charges N and J and their combinations Λ ± are functions of momentum k given by Eqs. (3), (4). The time integral is performed by contour integration and depends on the position of branch points in the complex plane of t. For a supersonic impurity v 2 (k) > v 2 s there is one branch point in the upper and one in the lower half-plane, which leads to a double-sided edge singularity in correlation functions. For a subsonic impurity v 2 (k) < v 2 s both branch points are in the upper half-plane, which leads to vanishing of the integral for ω < ε(k) and the result In both cases the edge exponent is Substituting the relations (54) into this expression one obtains Eq. (5).

B Properties of shift functions in Lieb-Liniger model
Similarly to the scattering phase θ(λ − µ), the shift function obeys F (λ|µ) = −F (−λ| − µ). We can also interchange the arguments using the following non-linear identity obtained by Slavnov in Ref. [52]. In addition, there are useful relation between the density of quasimomenta and the shift function. The first identity can be proven by the direct substitution of the left hand side into Eq. (12) and integrating the kernel K(λ − µ) to generate the phase shifts in the right hand side of Eq. (14). Taking λ = q one has 2πρ(q) = 1 + F (q|−q) − F (q|q) .
Multiplying both sides by 1 − F (q|q) − F (q|−q) and using Slavnov's identity Eq. (57) with λ = µ = q we arrive at another useful relation 1 2πρ(q) which was first established in Ref. [53]. Finally, using Slavnov's identity together with Eqs. (59), (60) and expression (64) for ρ(q) in terms of the Luttinger parameter, we can show that which is convenient for interchanging indices in shift functions.

C Relation of thermodynamic quantities of Lieb-Liniger model and its Bethe Ansatz solution
It is well known that solutions of Bethe Ansatz equations are related to thermodynamic quantities. These relations are demonstrated in Ref. [33]. For completeness we prove certain thermodynamic relations which we used in Section 4. We start by taking λ = q in Eq. (24) and using Eq. (59) to show that which is identical to Eq. (A.3.6) of Ref. [33]. For small momenta the dispersion relation of elementary excitations becomes linear so one can define the sound velocity where we used Eqs. (26). This leads to a remarkable exact value of the ground state density of quasimomenta at the edge, We obtain another useful relation, by differentiation of both sides of Eq. (18) with respect to h. Consider now the condition ǫ(q) = 0 which determines the dependence h(q). Differentiating both sides of this condition and taking into account Eq. (65) we arrive at This relation can also be obtained from the thermodynamic definition of sound velocity [33]. Eqs. (65), (66) can be used to calculate the derivative of the dispersion relation with respect to q. Combining them with Eq. (58) for the density of ground state quasi-momenta we obtain

D.1 The Resolvent
The equivalence of expressions (15), (16) for momentum and energy of an excitation and the corresponding expressions (17), (18) was proven in Ref. [33]. Below we use a similar method for proving this equivalence for reader's convenience. This will allow us to set up useful notations. We introduce the resolvent R(µ, ν) which solves the integral equation We rewrite this equation in the operator form which can be formally inverted leading tô or, equivalentlyÎ This identity leads to the formal solution of Eq. (12) for the ground state density of quasimomenta Here the vector ρ has elements ρ(λ) and, similarly, 1 has unity elements. The ground state density of rapidities in terms of the resolvent becomes Consider now the shift function F (ν|λ) which depends on λ as a parameter. Taking derivative of both sides of Eq. (14) with respect to λ and comparing with Eq. (68) one gets immediately ∂ λ F (ν|λ) = −R(ν, λ), so that Here we have used the fact that F (ν| − ∞) = πρ(ν). Substituting Eq. (74) into Eqs. (15) gives It is consequence of the theory of linear integral equations that R(ν, σ) = R(σ, ν). Using this fact together with Eq. (72) leads to Alternatively, the equivalence of Eqs. (15) and (17) can be proven by expressing the scattering phase in the right hand side of Eq. (15) in terms of the shift function using Eq. (14). We now show the equivalence of (16) and (18). Using Eq. (71) and Eq. (74) we can rewrite the latter as The integral in r.h.s. is similar to the one in Eq. (16) but has wrong sign and the arguments of the shift function F are interchanged. To manipulate them into the right order we use Slavnov's identity (57). This produces the result (16) plus an extra term, One recognises ǫ(±q) in the parenthesis of this expression which must be zero by the right choice of chemical potential h. We have therefore established the equivalence of Eqs. (16) and (18). The resolvent operator allows one to prove other useful identities. Let us start with differentiating with respect to q both sides of Eq. (12), Using the fact that ρ(q) = ρ(−q) and Eq. (68) the solution can be found at once: Consider now Interchanging the arguments of the resolvent and using Eq. (73) as well as Eq. (64) we get ∂ q n = 2ρ(q) + ρ(q) (2πρ(q) − 1 + 2πρ(−q) − 1) = 4πρ 2 (q) = K π .

D.2 Derivatives of energy and momentum of excitations in Lieb-Liniger model
By differentiating both sides of Eq. (14) with respect to q and using the expression (69) for the resolvent it is easy to show that so that Using this identity in Eq. (15) we can show that which upon substituting the result (64) becomes Eq. (25). Consider now the derivative of the dispersion relation with respect to quasimomentum λ at fixed q. Differentiating both sides of Eq. (18) we get Using the property ∂ λ K(λ − µ) = −∂ µ K(λ − µ) and integrating by parts leads to where we have used the fact ǫ(±q) = 0. The identity (71) allow to write the solution ∂ λ ǫ(λ) = 2λ + 2 q −q µR(λ, µ) dµ .