Ground-state energy and excitation spectrum of the Lieb-Liniger model : accurate analytical results and conjectures about the exact solution

We study the ground-state properties and excitation spectrum of the Lieb-Liniger model, i.e. the one-dimensional Bose gas with repulsive contact interactions. We solve the Bethe-Ansatz equations in the thermodynamic limit by using an analytic method based on a series expansion on orthogonal polynomials developed in \cite{Ristivojevic} and push the expansion to an unprecedented order. By a careful analysis of the mathematical structure of the series expansion, we make a conjecture for the analytic exact result at zero temperature and show that the partially resummed expressions thereby obtained compete with accurate numerical calculations. This allows us to evaluate the density of quasi-momenta, the ground-state energy, the local two-body correlation function and Tan's contact. Then, we study the two branches of the excitation spectrum. Using a general analysis of their properties and symmetries, we obtain novel analytical expressions at arbitrary interaction strength which are found to be extremely accurate in a wide range of intermediate to strong interactions.


I. INTRODUCTION AND MOTIVATION
The one-dimensional (1D) model of point-like bosons with arbitrary repulsive contact interactions was introduced more than half a century ago by Lieb and Liniger [2] as a generalization to finite interaction strengths of the Tonks-Girardeau gas (TG) [3] of hard-core repulsive bosons, and is known as the Lieb-Liniger (LL) model. In the thermodynamic limit, its exact ground state is given by a set of equations obtained in [2] and [4] by coordinate Bethe Ansatz (BA), a powerful technique introduced by Bethe to exactly solve the Heisenberg model [5]. Back in the 1960's, experimental realizations of 1D degenerate gases were not available, so the LL model was deemed as a toy model. Yet, it soon attracted attention in the context of quantum integrable models [6][7][8], being one of the simplest. In this respect, it helped understand various aspects of the many-body problem and why one-dimensional systems are so specific, the most appealing features being the effective fermionization of bosons at large interaction strength [3,9] and the existence of two branches of excitations [4], one of them reminiscent of the well-known Bogoliubov dispersion [10], the other linked with a quantum analog of classical solitons [11]. Meanwhile, it has served as a testbed for the universal description of low-energy properties of gapless 1D systems within the bosonization technique, known as the Tomonaga-Luttinger liquid framework [12][13][14][15][16][17][18].
The equilibrium grand canonical description of the LL model was provided by Yang and Yang who introduced the Thermodynamic Bethe Ansatz (TBA) [19], opening new investigation directions such as finite-temperature thermodynamics [20][21][22] and quantum statistics of the model [1,23]. The dynamical correlations of the LL model remained for a long time an open problem. The determination of its exact time-dependent density-density and phase-phase correlations and their Fourier transforms known as the dynamical structure factor and spectral function respectively, required the development of fairly involved algebraic Bethe Ansatz (ABA) techniques relying on the inverse scattering method [24]. The form factors were computed numerically using the ABACUS algorithm [25], first at zero [26,27], then at finite temperature [28], while parallel works focus on analytic and algebraic general considerations [29][30][31]. All of them contribute to validate a nonlinear extension of the Luttinger liquid theory developed approximately at the same time [32]. Very recently, new research lines are, among others, out of equilibrium dynamics with quantum quenches [33][34][35][36], few-body physics [37][38][39][40][41], and links with other models, such as attractive 1D fermions [42], a BCS model [43][44][45], the Kardar-Parisi-Zhang (KPZ) model [46], directed polymers [47], three-dimensional black holes [48] or the Yang-Mills problem on a 2-sphere [49], using mapping techniques, along with extended versions of the model to anyonic statistics [50], spatial exponential decay of the interaction [51] and a supersymmetric version [52].
As a matter of fact, the main reason of the renewed interest for the LL model in the early years of our century is experimental progress in cooling and trapping ultracold atoms. Strong confinement in two transverse dimensions [53][54][55] and low enough temperatures led to realizations of (quasi-)1D systems of bosons [56], some of them suitably described by the LL model [57]. Furthermore, the interaction strength is now tunable in a very controlled way [58], allowing to probe any interaction regime from weak to strong [59,60], as well as attractive ones [61], yielding a strongly excited state called 'super-Tonks-Girardeau' gas [62]. Experimental studies of the properties of the LL model are reviewed in [63,64]. In particular, the momentum distribution has been measured [59], as well as local two-body [65], three-body [66][67][68] and temporal correlations [58]. The phase diagram at finite temperature predicted in [69] has been explored in [70].
In this work, we focus on the theoretical and mathematical issues associated with the analytical description of the ground state. After so many years since the solution by Lieb and Liniger, the closed-form exact solution of the model is still lacking. More problematically, a wide range of experimentally relevant, intermediate repulsive interactions, is hardly accessed analytically by perturbation theory. This lack of accuracy affects the estimate of many physical quantities, such as the sound velocity or the Tan's contact of the gas. Moreover, approximation schemes and numerical solutions need benchmarks. A complete analytical understanding of the ground state, whose importance could be compared to the solution of the 2D Ising model, would be a benchmark towards solutions at finite temperature or generalizations to multicomponent systems, where the explicit results from integrability are even less accurate to date. From an analytical point of view, many properties of the LL model are easily studied in the TG limit of infinite interactions. At small enough interaction strengths, Bogoliubov expansion yields quite accurate results [2], and at arbitrary interactions the framework of Tomonaga-Luttinger liquid theory is appropriate to tackle many low-energy properties [71], but the validity range of both approximations is not known with accuracy.
The goal of this article is to provide analytical estimates for the ground state energy, excitation spectrum and other related observables in the experimentally relevant regime of interaction strengths γ ∈ [1, 10] as defined in Eq. (3) below. For this purpose, we will start from the strongly-interacting regime. We use a systematic method to obtain corrections from the TG limit and study its accuracy by computing several observables and comparing the results with the numerics and approximations available in the literature.
The structure of the paper is as follows: first, we introduce the main features of the model and briefly discuss the weakly-interacting regime. Then, focusing on the energy, we systematically evaluate corrections to the regime of infinite interactions up to order 20 and make conjectures about a possible resummation, using comparison with the numerics as an accuracy test. We compute various ground-state properties linked to this quantity and give results whose accuracy is unprecedented and already sufficient for all practical purposes. To finish with, we focus on the excitation spectra of the LL model and derive a simple expression which is very accurate in the whole range of repulsive interactions.

A. Model and main equations
The 1D quantum gas of N identical spinless point-like bosons of mass m with contact interactions confined to a line of length L is described by the Lieb-Liniger Hamiltonian H LL which reads in first-quantized form [2] where {x i } i∈{1,...,N } label the positions of the bosons, is the Planck constant divided by 2π, g 1D is an effective one-dimensional coupling constant which can be deduced from experimental parameters [53,54] and δ is the Dirac function. In the following we will focus on g 1D > 0, corresponding to repulsive interactions. Put into its second-quantized form, better suited to treat the thermodynamic limit, the Hamiltonian now reads [2] H where ψ is a bosonic field satisfying the canonical commutation relations with its Hermitian conjugate It is natural, then, to define a dimensionless interaction strength (also known as Lieb's parameter) as where n 0 ≡ N/L is the mean linear density of the gas. The interaction strength for a given atomic species is experimentally tunable by confinement-induced resonance, Feshbach resonance or dilution. The dependence on n 0 with negative power is rather counter-intuitive compared to higher dimensional bosons. For repulsive interactions, γ ranges in [0, +∞[. The case γ = 0 corresponds to noninteracting bosons, while the limit γ → +∞ corresponds to the Tonks-Girardeau gas. In this regime, an exact solution is available, based on a mapping on a noninteracting spinless Fermi gas, allowing for a simple determination of certain observables [3]. With periodic boundary conditions corresponding to a ring geometry and ensuring translational invariance, the Bethe hypothesis allows to derive exact expressions for the ground state energy, excitations and static correlations of the system. The procedure, called coordinate Bethe Ansatz, is well-known for this system and widely explained in the literature, we refer to [2,[72][73][74] for details and elegant physical discussions. At finite N this yields a system of N transcendental equations, which reduces to a much simpler set in the thermodynamic limit, and becomes analytic. The final result can be put in the form of a set of three integral equations, namely where g denotes the quasi-momentum distribution in reduced units and α is a non-negative parameter, in a one-to-one correspondence with the Lieb parameter γ introduced before via a second equation, The last equation of the set yields the dimensionless average ground-state energy per particle e, linked to the total energy E 0 by e(γ) ≡ 2m 2 E0(γ) N n 2 0 , according to Interestingly, the first equation of the system, Eq. (4), is decoupled from the others, which is a specificity of the ground state [19] and simplifies the problem to some extent. It is a homogeneous type II Fredholm integral equation with Lorentzian kernel [75], in the following we will refer to it as Lieb's equation. We mention in passing that it is also known in an other context as Love's equation [76][77][78], whose solution would yield the exact capacitance of a capacitor formed of two circular parallel plates, C(α) ∝ 1 −1 dzg(z; α) where α has then the physical meaning of an aspect ratio. In itself, it constitutes an extremely difficult open problem in engineering dating back to Clausius and Kirchhoff [79,80]. For a historical review of many attempts (-and failures) to solve this problem, we refer to [81]. Mathematicians were also inspired by this problem, in the context of convex analysis. Among their most important contributions to the problem are an inequality on C(α) [82] and a conjecture that its small α expansion should have a deep link with invariants of a Riemannian manifold [83].
To solve the system composed of Eqs. (4), (5) and (6) with high accuracy, the most crucial step is to solve the Lieb equation Eq. (4). It was done numerically by Lieb and Liniger for a few values of α spanning several decades of interaction strengths [2], and often again by other authors. Both numerically and analytically, the procedure based on the above equations to find the ground state properties of the LL model relies on the following steps. An arbitrary positive value is fixed for α, and Eq. (4) is solved, i.e. g is found with the required accuracy as a function of z ∈ [−1, 1]. Then, Eq. (5) yields γ(α), subsequently inverted to get α(γ), which may reveal thorny if one tries to do it analytically. Eventually, Eq. (6) is solved, yielding the ground-state energy per particle. Many interesting observables are then combinations of derivatives and moments of the latter. Note that they depend on a sole parameter (α or γ), which is the key of the conceptual simplicity of the model.
Both numerically and analytically, the case α 1 (corresponding to the weakly interacting regime γ 1) is known to be extremely difficult to solve compared to α 1 because of a singularity in the function g at α = 0 [2], whose physical interpretation is that noninteracting bosons are not stable in 1D. In the following, we will tackle quickly the weakly-interacting regime, mostly for the sake of completeness, and then dwell much deeper on the strongly-interacting regime.

B. Expansions in the weakly-interacting regime
Finding approximate analytical solutions of Eq. (4) at small values of the parameter α appears to be quite an involved task. A guess function was proposed by Lieb and Liniger [2], namely which is a semi-circle law and was derived rigorously in [78,84]. Heuristic arguments have suggested the following correction far from the edges in the variable z [84,85]: reproduced later by direct calculation in [86] after regularization of divergent series. We are not aware of other analytical studies claiming better accuracy in the weakly-interacting regime. A technical difficulty consists in evaluating precisely the domain of validity of both approximations Eq. (7) and Eq. (8) in α, and if the latter offers a significant improvement in accuracy. A systematic comparison to numerical solutions at α 1 is beyond the scope of this work. We can nevertheless, to a certain extent, discriminate two given analytical solutions in terms of their accuracy at a given value of α. To this end, we define a 'local error' functional as and a 'global error' functional If a proposed solution g(z; α) is exact for a given fixed parameter α, then trivially [g; α] as a function of z is uniformly zero. A sufficient condition for a proposed solution g to be more accurate than an other solutiong is that | [g; α]| ≤ | [g; α]| for all z in [−1, 1]. The global error functional yields a good accuracy criterion, by requiring that it is lower than a threshold. Both quantities can be used for numerical as well as analytical purposes. We used them to check that the correction in Eq. (8) improves locally the accuracy with respect to Eq. (7) close to z = 0 if α 1.
However, close to |z| = 1 we find that Eq. (8) is not necessarily more accurate. For an other criterion specifically designed to deal with the case α = 1, we refer the reader to Appendix A. From those approximate solutions, one can then derive the dimensionless energy. The highest-order exact expansion to date in the weakly-interacting regime is The structure of this expression suggests that the appropriate variable to write an expansion of the energy in the weakly-interacting regime is γ 1/2 . The first term of the right-hand side is readily obtained from Eq. (7), the second one is found when using next order in Eq. (8) and coincides with the result from Bogoliubov's approximation [2]. The prefactor of the γ 2 term is more controversial. Its expression was inferred on numerical grounds in [87], while in [88] a factor 1/8 instead of 1/6 is found analytically after lengthy calculations. It seems, though, that the factor 1/6 is correct, as it was recently recovered from potential theory in [89] in a quasi-rigorous way, and coincides with numerical independent calculations carried out in [90]. The same author evaluated numerically the exact fourth term, given by multiple integrals, and found D 0.0016 [91], in agreement with a fit on very accurate numerics obtained in [92]. In a graph from [91], we do not see any difference between Eq. (11) and the exact numerical result up to γ of the order of 4 − 4.5. It also shows that Bogoliubov's approach is numerically exact only up to γ 0.5 for the energy.
Our main contribution to the problem is purely numerical in this regime. We evaluated e numerically for 15 values of the interaction strength spanning the interval γ ∈ [1, 10], a fit containing the analytically known prefactors and four adjustable coefficients yields with a relative error as low as a few per thousands in the interval considered.

C. Expansions in the strongly-interacting regime
Since we aim at describing correctly the experimentally relevant regime of intermediate interactions, a complementary approach consists in focusing on the strongly-interacting case, then decreasing interactions as much as we can to reach this regime and cross-check the results found at low interactions. This approach is in principle easier to handle, both analytically or numerically, because one starts far from the singularity at α = 0, but on the other hand it requires higher-order expansions.
When interactions are strong, a good way to obtain accurate analytical solutions to the Lieb equation Eq. (4) consists in a double series expansion in z and 1/α around (0, 0). To do so, we use a procedure whose main ideas are found in [93], and fully developed in [1]. It has the advantage of yielding several orders of perturbation theory at each step and of being automatically consistent at all orders in z and 1/α. We refer to Appendix B for a quite detailed derivation and mathematical discussions. The expansion procedure yields an approximate solution of the form where M is an integer parameter, g jk are rational numbers and g(z; α, M ) converges simply to g(z; α) for M → +∞. In Eq. (13) one immediately sees that the maximum exponent of z 2 varies more slowly with M than the one of 1/α, which is the main drawback of the method. Indeed, as we will see later on, in order to obtain accurate predictions it is very important to capture the correct behavior of g as a function of z in the whole interval [−1, 1], but at the extremities it converges slowly to its true value since the Taylor expansion in z is done around the origin. It leads to a strong limitation of validity in 1/α at a given order, which we address in detail in the current work. We remark that, as a rule of thumb, the absolute value of the local error defined in Eq. (9) increases when α decreases and when z increases. Thus, very high order asymptotic expansions may be required to get an accurate description up to α = 2, which is the lowest attainable value in this approach. A good point, however, is to notice that α = 2 corresponds to an interaction strength γ ≡ γ c 4.527 [1], which is low enough to combine with the result at small interaction strengths, thus obtaining an accurate description of the ground state of the Lieb-Liniger model over the whole range of repulsive interactions.
In order to check the accuracy of the analytical results, at fixed α we compare the various approximations we find at increasing M with an accurate numerical solution obtained with a Monte-Carlo integration algorithm stopping when the condition on the global error functional E[g; α] < h is fulfilled, where h is a threshold value fixed at 10 −5 for α > 2. We only need to compare the values at z = 1, which are the most difficult to attain analytically, to get an idea of the accuracy of the expansion. A systematic comparison at increasing values of M is shown in Fig. 1. As explained in Appendix B, at fixed M one obtains simultaneously the expansions to orders k m = 2M +1 and 2M +2 in 1/α. Since we look for fully analytical solutions we are limited to relatively low values of M . We solved the Lieb equation up to M = 9 and checked that our result agrees with the Supplementary Material of [1], where it is given to order M = 3. In particular, again we study the worst case for our method, α = 2, where the convergence is the slowest with M because this value lies on the border of the convergence domain. In [93] it was estimated that an expansion to order k m 40 is needed to reach a 5-digit accuracy. We find that at k m = 20 the relative error is still as large as 8/1000. As a general fact, approximate solutions oscillate with M around the exact solution, and odd orders are more accurate, as pointed out in the same reference.
Quite remarkably, doing the arithmetic average of two consecutive even orders in k dramatically increases the accuracy, and this procedure can be repeated, allowing to treat lower and lower values of α. We use the notation g m (z; α, M ) ≡ [g(z; α, M ) + g(z; α, M −1)]/2 to denote the mean value obtained by averaging over two consecutive orders in M . As can be seen in Fig. 2, doing the mean of 18 th and 20 th orders is enough to cover the whole range [2, +∞[ with an accuracy comparable to the one of the numerics.

D. Conjectural resummations in the regime of intermediate interactions
A major drawback of the method used in this work is that it quickly yields too unhandy expressions for the function g to be given here. Indeed, at order M , there are i(i + 1) terms of order z 2(M +1−i) in the expansion for each i ∈ {1, . . . , M + 1} plus a common term 1/(2π) which is the value in the Tonks-Girardeau regime, yielding a total of 1 + (M + 1)(M + 2)(M + 3)/3 terms for a given value of M . This motivated us to seek compact representations for the function g(z; α, M ) allowing to easily use those expansions for further applications and to generate them up to large orders.
To do so, we shall consider all the mathematical operations yielding the expansion at a given order M as a black box, and try to recognize patterns in the coefficients g jk of the various terms appearing in the final result. In doing so, we were able to guess some expressions which, we conjecture, are valid at any order. Our work is a first step towards the exact solution of the Lieb equation (4), which would be achieved if one could find all the terms of the series expansion and resum them. As we will show below, our conjectures are in excellent agreement with the numerically exact solution and hence useful for practical applications.

Conjectures for g
By a careful analysis of the terms of the series expansion, we found two apparently distinct groups of patterns, which we arbitrarily call 'first kind' and 'second kind' respectively. The terms of first kind already arise at low orders in z and 1/α, while terms of second kind make their appearance at higher orders and should play a crucial role in the crossover region from strong to small interactions. While some are trivial to figure out, others are far more difficult to find because of their increasing complexity. Denoting by m the 'kind' and by n an index for the 'complexity', we write where each I M m,n is itself a double sum over all terms of given kind m and complexity n that appear at order M . To make things clearer, we illustrate our point by giving examples found at order M = 9.
The only term of first kind and complexity n = 4 we have identified is which should correspond to the terms with index j = 0 in I M 1,4 . Terms of second kind and lowest complexity are We also found . We note, however, that some terms of first kind could be interpreted as second kind, thus involving a shift of the summation index in one of them.
If one performs the summation of all terms considered here, the resulting expansion is exact up to order 1/α 11 included, thus the structures that we conjecture are efficient to encode many terms of the series in a relatively compact way. They are all easily resummed in the limit M → +∞. However, even at so high order, the expansion proposed here still has an error of the order of 1 − 2% at α = 2 as can be seen in the right panel of Fig. 1. Moreover, we were not able to find the generating term allowing to go from complexity n to n + 1, thus leaving open the question of the full resummation of the series.

Conjectures for e
It is also possible, and certainly even more fruitful for the physical problem we are concerned with, to look for patterns directly on the dimensionless ground-state energy e(γ). Here, we shall write it e(γ) = +∞ n=0 e n (γ), where once again the index n denotes a notion of 'complexity'. If we focus on the large γ expansion of e(γ) in 1/γ, we easily identify a sequence of 2M + 3 terms. We assume that they will appear at all orders and thus resum the series, obtaining in agreement with [94] where e T G = π 2 /3 is the value in the Tonks-Girardeau regime.
Guided by the expansion of e(γ) in 1/γ up to M = 9 (given in Appendix C with numerical coefficients) and the property we conjecture that the general form of the term of complexity n ≥ 1 defined above is where L n is a polynomial of degree n − 1, whose rational coefficients are all non-zero and of alternate signs. We also conjecture that the absolute value of the coefficient of the highest-degree monomial of L n is 3 * 2 2n+3 (n+2)(2n+1)(2n+3) . Interestingly, those partially resummed terms are not divergent at small values of γ, thus allowing to make predictions at higher 1/γ than the corresponding expansion at a given order, which converges very slowly to the exact result as M increases. We also noticed that e 0 corresponds to Lieb and Liniger's approximate solution [2], while the first correction e 1 was predicted in [88] with a different approach, giving more weight to our conjectures. In outlook, there are two steps left to exhibit the exact analytical solution of the Lieb-Liniger model. First, one should identify the generating function associated to the L n polynomials, then resum the series.

E. Illustrations and physical analysis
Now that we have found analytical expressions over the whole range of interaction strength, we illustrate them on a few observables of interest and give physical comments. First, in Fig. 3 we show the mean energy per particle over a wide range of repulsive interactions. Our expressions in both the weakly-interacting regime as given by Eq. (12) and the strongly-interacting regime as given by Eq. (17) are compared to the numerics to emphasize their accuracy. The analytical expressions have a wide numerical overlap in the regime of intermediate interactions and are quasiindiscernable from the numerics at this scale. In particular, extrapolating the expansion in the strongly-interacting regime to low values of γ, we obtain a substantial improvement compared to the previously known approximations e 0 and e 0 + e 1 in Eqs. (15) and (17) when γ > ∼ 1. It is then possible to evaluate what are the respective parts of the mean kinetic energy e k and interaction energy e p per particle, from Pauli's theorem [2] e p (γ) = γ de dγ (18) and e k (γ) = e − γ de dγ .
We show them in the left panel of Fig. 4, normalized to the total energy in the Tonks-Girardeau regime as in [95]. The kinetic energy is maximal in the Tonks-Girardeau regime of ultra-strong interactions. This can be seen as a manifestation of fermionization, since in several respects the particles behave as free fermions in this limit due to the Bose-Fermi mapping [3]. Actually, the ratio of interaction to kinetic energy scales as γ −1/2 in the weakly-interacting regime and decreases monotonically when interactions become larger. We remark that, contrary to common wisdom, γ does not represent this ratio. It would be true in mean-field theory only, and the curve in the right panel clearly shows that this picture is not appropriate.
Other interesting quantities are the local k-particle correlation functions g k defined as where . . . represents the ground-state average. They define the local many-body correlations, that is the expectation value of a product of 2k field operators (k = 1, 2, 3 . . . ) at zero space and time separation. The local pair correlation g 2 (respectively g 3 ) is a measure of the probability of observing two (three) particles at the same place. In particular, g 3 governs the rates of inelastic processes, such as three-body recombination and photoassociation in pair collisions. The function g 2 is easily obtained with our method using the Hellmann-Feynman theorem [96,97] yielding g 2 = de dγ [98]. This relation shows that the fact that e is an increasing function of γ is actually a direct consequence of the positiveness of g 2 . The three-body correlation function g 3 is obtained from the moments of the g function, defined as by [99] g 3 (γ) = 3 2γ We have already found the second moment, which is nothing else but the energy e. Upon computing the 20 th order expansion in 1/γ at high interaction, we found excellent numerical agreement with a very accurate fit function proposed in [99] down to γ 5 − 6. In principle, we could infer an expression for the fourth moment the same way as we did for e and propose a conjecture for the expansion, yet it would not be very useful in this context, since very accurate expressions for g 3 covering the whole interaction range are already given in [99]. We use those expressions in what follows. In Fig. 5 we plot g 3 , and our analytical expression of g 2 readily obtained from Eqs. (12) and (17), as well as their ratio as functions of γ. The fact that g 2 → γ→+∞ 0 (known as anti-bunching) is once again a consequence of fermionization. For very low values of γ, the ratio g 3 /g 2 can not be neglected, which means that the weakly-interacting 1D Bose gas is less stable than at higher interaction strengths. At high interactions fermionization and its associated Pauli principle preclude three-body processes in this system [98]. This property has been the key to realize the TG gas experimentally.
Finally, we study the first-order correlation g 1 (x) ≡ ψ † (x)ψ(0) n0 , describing the spatial correlations of the gas. Historically, it was first computed in the Tonks-Girardeau regime [100][101][102], where expansions at short and long distances are now known to high enough orders to match at intermediate distances, as can be seen in Fig. 6. Here we use the notation z = k F x, where k F = πn 0 is the Fermi wavevector in 1D. At large distances the expansion is [102](with sign mistakes corrected by Gangardt [103]) where G is the Barnes function defined as G(1) = 1 and the functional relation G(z +1) = Γ(z)G(z), Γ being the Euler Gamma function. Equivalently [101], G(3/2)/ √ 2 = π √ e2 −1/3 A −6 where A = 1.2824271 . . . is Glaisher's constant. At short distances, using the same technique as in [104] to solve the sixth Painlevé equation, we find the following expansion, where we added a few terms compared to [102]: The first sum is a truncation of the integer series defining the function sin(z) z , which is trivially g F 1 (z), the density matrix of free fermions. This is a typical example where a physical quantity in the Tonks-Girardeau regime is not the same as for free fermions, because it depends on the phase of the wavefunction and not only on its modulus.
At finite interactions however, analytical expansions are not accurate enough to match at intermediate distances in spite of some advances [105,106], so they are often treated numerically, e.g. with Monte-Carlo techniques [107], or with the ABACUS algorithm [27]. At short distances, one can write where the first coefficients are explicitly found as [106] c 1 = 0, c 2 = − 1 2 e k and c 3 = 1 12 γ 2 de dγ , while explicit expressions for higher order terms are still unknown. Once more, we readily have c 2 and c 3 with excellent accuracy. In Fig. 7 we show the short-distance asymptotics of g 1 to third order in z, for several values of γ. For all non-zero values of γ the g 1 function decays at long distances, showing the absence of true off-diagonal long-range order, especially for high values of γ. This is in agreement with the Mermin-Wagner theorem, a consequence of which is the impossibility of Bose-Einstein condensation in the LL model at finite values of γ.
The Fourier transform of the one-body correlation function g 1 is the momentum distribution of the gas, n(k) = n 0 +∞ −∞ dxg 1 (x)e −ikx . Even in the Tonks-Girardeau limit where g 1 (z) is known with excellent accuracy at any distance, the matched asymptotics are not suitable to perform the Fourier transform analytically. One can, however, construct the curve point by point for any value of k by numerical integration. Moreover, for any value of γ, the large momenta asymptotics is such that k 4 n(k) n 4 0 → k→+∞ C(γ) = γ 2 de dγ [106,108], where C is known as Tan's contact [109], shown in Fig. 8.

III. EXCITATION SPECTRUM, EXACT AND APPROXIMATE RESULTS
Now that we have found accurately the ground-state energy of the model, we can tackle a more complicated and partially open problem, namely the characterization of excitations above the ground state at zero temperature. A peculiarity of the LL model is that it features two excitation spectra as first pointed out by Lieb and Liniger themselves [4], historically called type I and type II respectively. in the Tonks-Girardeau regime as a function of the dimensionless distance z. Short-distance asymptotics given by Eq. (24) at ninth and fourteenth order (red and red, thick respectively) and long-distance asymptotics given by Eq. (23) (black, dashed). Adding a few corrections at small distances does not significantly improve the range of validity of the approximation.  In the Bethe Ansatz solution at finite N , the total momentum and energy of the system are given by P = N j=1 k j and E = 2 2m N j=1 k 2 j respectively [4], where the set of quasi-momenta {k j } j∈{1,...,N } satisfies the following system of N transcendental equations The I j 's, called Bethe quantum numbers, are integer for odd values of N and half-odd if N is even. Since we consider γ > 0, the quasi-momenta are real and can be ordered so that k 1 < k 2 < · · · < k N . Then, automatically I 1 < I 2 < · · · < I N [4]. The ground state corresponds to I j = − N +1 2 + j and has total momentum P GS = 0. In what follows, we will use the notations p ≡ P − P GS and ≡ E − E GS for the total momentum and energy of an excitation by reference to the ground state, so that the excitation spectrum is defined as (p). For symmetry reasons, we will only consider excitations such that p ≥ 0, those with −p having the same energy.
To give a physical meaning to the type I and type II excitations, we treat the simplest case, i.e. the Tonks-Girardeau regime. In the limit γ = +∞, the set of equations (26) decouples, and its solution becomes trivial, yielding k j = 2π L j. The simplest excitations are of two different natures. In a 'one-particle' excitation, the highest-energy particle labeled by j = N gains a momentum n 2π L , with n an arbitrary non-negative integer, compared to the ground state. This corresponds to the type I excitations. In the other type of elementary excitation, corresponding to type II, a set of n ∈ {1, . . . , N } rapidities starting at index j and belonging to adjacent momentum states are all shifted by 2π L , leaving a void at j in the set of Bethe numbers, interpreted as a hole, hence the name 'one-hole' excitation. By construction, to each of the N one-hole excitations corresponds a one-particle excitation, the two processes being identical for the lowest value of momentum gain. There is no theoretical limit to the number of particle excitations, and their energies are always larger than those of the hole ones. By construction, the one-hole excitation of highest momentum is such that p = 2π N L = 2p F , where p F denotes the Fermi momentum in the thermodynamic limit. This particular case can be interpreted as an umklapp process, where a particle on the left of the Fermi surface (which consists in two points in 1D) jumps to the other side of the latter with a small gain of energy. To finish with, we also mention that any combination of one-particle and one-hole excitation is possible, with intermediate energy. In the thermodynamic limit, they form a continuum between the type I and type II branches, respectively associated to one-particle and one-hole processes as we have just seen, known as the particle-hole continuum.
In the Tonks-Girardeau regime, one easily finds that the ground-state energy is , in agreement with the result in the thermodynamic limit, while the relative momentum and energy of the j th excitation with respect to the ground state are respectively p j = j 2π for holes and for particles. The variable j can be eliminated, one finds that all points lie on the curves and They cross exactly once, at the lowest allowed excitation mentioned before, whose energy is strictly zero in the thermodynamic limit, where we recover the well-known result for the thermodynamic limit of the excitation spectrum for noninteracting spinless fermions [3]. Meanwhile, the energy associated to the umklapp process vanishes, thus the system is gapless. Moreover, the type II dispersion then acquires the symmetry p ↔ 2p F −p. We illustrate all those points in Fig. 9, representing the type I and type II excitations for N = 4, 10 and 100 bosons, compared to the result in the thermodynamic limit.
In the general case of finite interaction strength, the system of equations (26) for the many-body problem can not be solved both analytically and exactly, although expansions in the interaction strength can be obtained in the weaklyand strongly-interacting regimes [38]. The solution is easily obtained numerically for a few atoms, yet to reach the thermodynamic limit with several digits accuracy the Tonks-Girardeau treatment suggests that N should be of the order of 100, moreover the interplay between interactions and finite N may slow down the convergence at finite γ [26]. Therefore, here we shall directly address the problem in the thermodynamic limit, where it reduces to two equations [1,110] : and Q(γ) ≡ n 0 /[ 1 −1 dyg(y; α(γ))], known as the Fermi rapidity, represents the radius of the quasi-Fermi sphere associated with the bosons with repulsive interactions, and the function f satisfies the integral equation that we will henceforth call 'second Lieb equation'. We solve it with similar techniques as we did for the Lieb equation, Eq. (4), a few details of which can be found in Appendix D.
The excitation spectra at a given interaction strength γ are given by the parametric plots (k; γ)[p(k; γ)], k ∈ [0, +∞[. With this representation, one can interpret the type I and type II spectra as a single curve, where the type I part corresponds to |k|/Q ≥ 1 and thus to quasi-particle excitations, while type II is obtained for |k|/Q ≤ 1, thus from processes taking place inside the quasi-Fermi sphere, which confirms that they correspond to quasi-hole excitations. Basic algebra on Eqs. (29) and (30) quickly yields the following interesting general results.
(i) The ground state (p = 0, = 0) trivially corresponds to k = Q(γ), showing that Q corresponds indeed to the edge of the Fermi surface.
(ii) The point k = −Q(γ) corresponds to the umklapp point (p = 2p F , = 0), always reached by the type II spectrum in the thermodynamic limit, regardless of the value of γ. In other words, the fermionic structure found in the Tonks-Girardeau regime is stable to finite interactions, and the Lieb-Liniger model is gapless.
(iii) The maximal excitation energy associated to the type II curve lies at k = 0, corresponding to p = p F .
, preserving the symmetry found in the Tonks-Girardeau regime.
(v) The type I curve I (p) repeats itself, starting from the umklapp point, shifted by 2p F in p. Thus, what is sometimes considered as a continuation of the type II branch can also be thought as a shifted replica of the type I branch.
(vi) Close to the origin in p, I (p) = − II (−p). This can be proven using the following sequence of equalities based on the previous properties: These properties will reveal most useful in the analysis of the spectra, they also provide stringent tests for numerical solutions. With the expansion method used before, we can obtain the type II curve parametrically, and even explicitly, with excellent accuracy provided α > 2. As far as the type I curve is concerned, however, we are not only limited by α > 2, but also by the fact that our approximate expressions for g(z; α) and f (z; α) are valid only if |z − y| ≤ α, ∀y ∈ [−1, 1], thus adding the validity condition, |k|/Q(α) ≤ α − 1. The latter is not very constraining as long as α 1, but for α 2 the best validity range we can hope is very narrow around p = 0. To bypass this problem, Ristivojevic used the iteration method to evaluate g and f [1]. However, in practice this method is applicable only for large interactions since it allows to recover only the first few terms of the exact 1/α expansion of (k, α) and p(k, α) (to order 2 in [1]). Moreover, the obtained expressions are not of polynomial type, it is then a huge challenge to substitute the variable k to express (p) explicitly, and one has to use approximate expressions at high and small momenta. Therefore, we rather considered the high-p regime of I (p) semi-numerically, first computing which allows us to obtain the type I spectrum with excellent accuracy for all values of p, as illustrated in Fig. 10. At small momenta, in the general case, due to the analytical properties of both g and f , for all values of γ the type I curve can be expressed as a series in p [111,112], In the Tonks-Girardeau regime, as follows from Eq. (27), one has v s = v F = k F /m, m * = m, and all other coefficients vanish. The parameters v s and m * , which can be seen as a renormalized Fermi velocity and a renormalized mass, have a deep physical meaning. Linear Luttinger liquid theory predicts that v s is the sound velocity associated with bosonic modes at very low p [17]. To all accessed orders in 1/γ, we have checked that our expression agrees with the thermodynamic equality [4] v Analytical expansions for the sound velocity are already known at large and small interaction strengths. The firstand second-order corrections to the Tonks-Girardeau regime in 1/γ are given in [71], they are calculated to fourth order in [72] and up to eighth order in [1]. In the weakly-interacting regime, expansions are found in [71] and [86]. Our expressions for e allow us to considerably increase the accuracy compared to previous works after easy algebra. Interestingly, there is actually no need to know e(γ) to find v s . It is enough to know the g function at z = 1 in the first Lieb equation due to the equality [24] v Reciprocally, knowing v s yields an excellent accuracy test for the g function, since it allows to check its value at the border of the interval [− 1,1] where it is most difficult to obtain. Here we use both approaches to find the sound velocity over the whole range of interactions with excellent accuracy. Our results are shown in Fig. 11. One can see that v s → γ→0 0, thus we shall find g(z; α) → z→1,α→0 +∞, which shows that the method of polynomial expansion fails at too low interaction strengths, as expected due to the presence of the singularity. Moreover, this argument automatically discards the approximate expression given in Eq. (7) close to the boundaries in z. Linear Luttinger liquid theory assumes that, for all values of γ, I/II (p) p/p F 1 v s p and II (p) |p−2p F |/p F 1 v s |p − 2p F |. This strictly linear spectrum is however a low-energy approximation, and nonlinearities cannot be neglected if one wants to deal with higher energies. Including the quadratic term and neglecting higher-order ones is actually a complete change of paradigm, from bosonic to fermionic excitations at low energy, at the basis of the Imambekov-Glazman theory of Nonlinear Luttinger liquids [113,114]. In this approach m * is interpreted as an effective mass, whose general expression is [1,115] We have verified that it is in agreement with our expression at all accessible orders in 1/γ. The dependence of the effective mass as obtained from Eqs. (12), (17), (34) and (36) is shown in Fig. 11. Comparing the two panels of this figure, one can see that at large interaction strengths the sound velocity and the inverse effective mass are of same order. As already pointed out in [1], this supports the fermionic picture rather than the bosonic one. For the type II spectrum, the properties that we have found suggest an other type of expansion, also used in [116]: Using the property (vi), on the other hand, allows us to write Equating both expressions to order p 2 , one finds that A similar result was recently inferred on purely numerical grounds in [117]. We can even go a step further and show that which is one of our main results. Systematic suppression of the variable k and higher-order expansions suggest that, at large enough values of γ at least, higher-order terms can be neglected in expansion (37). We have checked numerically the agreement between the exact solution and this truncated result in Fig. 12. At first sight the result at order one is satisfying, but the second order correction significantly improves the result at intermediate values of the Lieb parameter. Overall, we dispose of very accurate analytical predictions for the spectra, for γ ∈ [1, +∞[. Third and higher order corrections are negligeable in a wide range of strong interactions.

IV. CONCLUSIONS AND OUTLOOK
In conclusion, first we have solved with high accuracy the set of Bethe-Ansatz equations established by Lieb and Liniger for the ground state of a 1D gas of point-like bosons with contact repulsive interactions in the thermodynamic limit, thus obtaining the distribution of pseudo-momenta, the average energy per particle and all related quantities. Our main result in this part consists in two simple analytical expressions which describe to a good accuracy the ground state energy. Their combination spans the whole range of interaction strengths, thus providing an alternative to the use of tabulated values. We have also made several conjectures towards the exact closed-form solution of the model valid for all interaction strengths, and introduced various accuracy criteria to test analytical and numerical solutions. Then, we have studied the two branches of the excitation spectrum and expressed the main contributions in terms of the sound velocity and effective mass, which can be obtained from the ground-state energy through thermodynamic relations.
A first natural generalization of our study concerns finite-temperature effects. The correlations have already been studied extensively, g 2 is known non-locally and at finite temperature [118][119][120][121], and the 3-and more particle correlations too [122,123], yet analytical expressions are still limited to small ranges of temperature and interaction strength. Our results for the g function can also be taken as a first numerical input in Monte-Carlo programs to solve the Bethe Ansatz equations at low temperature [124].
The harmonic trap used in most of current experiments would destroy the integrability by breaking translational invariance, but its effect on correlation functions can be studied numerically [125,126] or within the local density approximation (LDA) [127]. In this respect, our improved results for the homogeneous gas can be used to increase the accuracy of theoretical predictions, and for comparison with exact results [128] to test the validity of the LDA.
The very accurate expressions we have obtained for the excitation spectra may reveal useful to better understand the link between type II excitations and quantum dark solitons. They have been mostly studied in the weakly-interacting regime so far [129][130][131][132][133][134][135], but they may be a general feature of the model [136][137][138][139]. Excitation spectra also yield the exponents governing the shape of the dynamical structure factor near the edges in the framework of 'beyond Luttinger liquid' theory [32,114].
All our techniques can be adapted easily to the metastable gaseous branch of bosons with attractive interactions [140][141][142], to study the super Tonks-Girardeau behavior of the model with increased accuracy. Extensions of the current method may be used to study the Yang-Gaudin model of spinful 1D fermions, which has attracted much attention recently [143][144][145][146][147][148][149]. Furthermore, since the Lieb-Liniger model can be seen as the special case of infinitely many different spin values [150], it allows to check the consistency with the general case, which is far less well understood, and to make approximate predictions at the highest experimentally relevant values of the number of spin components.
1 0 , where β is the beta function, defined as β(x) ≡ 1 2 ψ x+1 2 − ψ x 2 , with ψ the logarithmic derivative of the Euler Γ function, also known as the digamma function, one naturally defines an error functional such that if g is exact, For instance, in [76] the following expression was proposed: Then |Err[g]| 0.0032 g(0, 1). Our numerical result is very close to that function. Fitting it by an eighth-degree polynomial, we find exactly the same value for the error up to 4 th digit. This allows us to check once more the accuracy of our numerical algorithm. In this Appendix we explain in more details than in the original article [1] Ristivojevic's method to systematically find approximate solutions to Eq. (4) in the strongly-interacting regime. First, we recall some qualitative features of the function of interest, g. In [2] it was shown that at fixed α, g as a function of z is positive, bounded, unique and even. Moreover, it is analytic provided α > 0.
To make the system of equations finite, we cut off the series in n at an integer value M ≥ 0. The system Eq. (B8) truncated at order M can then be recast into a matrix form: where A is a (M + 1) × (M + 1) square matrix, which is inverted to find the set of coefficients a 2n (α). Actually, one only needs to compute (A −1 ) i1 , ∀i ∈ {1, . . . , M +1}, then combined with Eq. (B1) to obtain the final result at order M . For full consistency with higher orders, one shall expand the result in 1/α and truncate it at order 2M +2, while it goes to order 2M in z.
Appendix C: Expansion of the dimensionless energy per particle in the strongly-interacting regime Using the method of Appendix B up to M = 9 yields the expansion of e(γ) to order 20 in 1/γ. The expression is given here with numerical coefficients for the sake of compactness, and reads e(γ) 1.
It is illustrated in Fig. (13). To find the type I and type II dispersion relations, we need to solve Eq. (31) copied here for convenience: Adapting the techniques of Appendix B, using odd-degree Legendre polynomials because f is an odd function of z, we write f (z; α) = where B is a (M + 1) × (M + 1) square matrix. Actually, one only needs to compute B −1 i1 , ∀i ∈ {1, . . . , M +1}. For the same reasons as before, the expansion is valid for α > 2, i.e. in the strongly-interacting regime. At fixed α it converges faster than g to the exact value when M is increased. Once again, from the expansions obtained one can guess patterns and write f (z) = z + m,n J M m,n (z, α).