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

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

Guillaume Lang Université Grenoble Alpes, LPMMC, F-38000 Grenoble, France CNRS, LPMMC, F-38000 Grenoble, France    Frank Hekking Université Grenoble Alpes, LPMMC, F-38000 Grenoble, France CNRS, LPMMC, F-38000 Grenoble, France    Anna Minguzzi Université Grenoble Alpes, LPMMC, F-38000 Grenoble, France CNRS, LPMMC, F-38000 Grenoble, France
###### Abstract

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 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 repulsive contact interactions has been introduced in Lieb () as a generalization to finite interaction strengths of the Tonks-Girardeau (TG) gas of hard-core repulsive bosons Girardeau (), and is known in the literature as the Lieb-Liniger (LL) model or the -Bose gas. Its exact ground state is encoded in a set of coupled equations obtained in Lieb () and LiebII () by coordinate Bethe Ansatz (BA). The LL model is one of the simplest quantum integrable models McGuire (); YangC (); Baxter (); Thacker (). Its exact solution has helped understand various aspects of the many-body problem in one dimension, the most appealing features being the effective fermionization of bosons at large interaction strength Girardeau (); Yukalov () and the existence of two branches of excitations LiebII (), one of them reminiscent of the Bogoliubov dispersion Lieb (); Bogoliubov (), the other linked with a quantum analog of classical solitons Faddeev (). The equilibrium grand canonical description of the LL model has been developed by Yang and Yang who introduced the Thermodynamic Bethe Ansatz YangYang (), thereby opening new investigation lines, such as finite-temperature thermodynamics Wang (); Guan (); Vogler (); Rosi () or quantum statistics of the model Tomchenko (). Later on, Haldane used the Lieb-Liniger model as a testbed for the universal description of low-energy properties of gapless 1D systems within the bosonization technique Haldane (), known as the Tomonaga-Luttinger liquid (TL) framework Tomonaga (); Luttinger (); Mattis (); Luther (); Efetov (); Giamarchi (); Cazalilla ().

The calculation of the dynamical correlations of the LL model remained for a long time an open problem. The determination of its exact time-dependent density-density correlations and their Fourier transform, known as the dynamical structure factor, is easily tractable in the Tonks-Girardeau regime only VignoloMinguzziTosi (). Perturbation theory allows to tackle the strongly-interacting regime ChernyBrand (), and the TL predictions at arbitrary interaction strengths Haldane (); Haldanebis (); Pitaevskii1 () are accurate only within a small low-energy range Pereira (); Lang1 (). Finding the exact solution at arbitrary interaction strength actually required the development of fairly involved algebraic Bethe Ansatz (ABA) techniques relying on the quantum inverse scattering method Korepin (). The form factors were computed numerically using the ABACUS algorithm CauxA (), first at zero Calabrese (); Slavnov (), then at finite temperature CauxT (), while mathematically-oriented works focus on analytic and algebraic general considerations Terras (); Maillet (); Shashi (). All of them tend to validate a nonlinear extension of the Tomonaga-Luttinger liquid theory developed in parallel Glaz ().

The experimental progress in cooling and trapping ultracold atoms has led to a renewed interest for the LL model. Strong confinement in two transverse dimensions Olshanii (); Dunjko (); Seiringer () and low enough temperatures have led to realizations of (quasi-)1D systems of bosons Gorlitz (), some of them suitably described by the LL model Moritz (). The possibility to tune interaction strength in a controlled way Guarrera () has allowed to probe any interaction regime from weak to strong Paredes (); Kinoshita (), as well as attractive ones Haller (), yielding a strongly excited state called ’super-Tonks-Girardeau’ (sTG) gas Boronat (). For a review of experimental studies of the properties of the LL model we refer to Xi-Wen (); Zwerger (). In particular, the momentum distribution has been measured Paredes (), as well as local two- Wenger () and three-body correlations Tolra (); Jacqmin (); Nagerl (). The phase diagram at finite temperature predicted in Petrov () has been explored in Bouchoule (). More recently, measurements of the dynamical structure factor have validated the ABA predictions in a regime where the standard Tomonaga-Luttinger liquid approach is not applicable Clement (); Meinert (). In current experiments, it is also possible to investigate out-of-equilibrium dynamics Kinoshitacradle (); Atas (). From a theoretical point of view, the LL model allows to address a rich variety of topics, e.g. the effect of quantum quenches Kormos (); Cartarius (); Essler (); Andrei (); Holzmann (), few-body physics Muga (); Oelkers (); Sakmann (); Makin (); Davis (), extended versions of the model to anyonic statistics Patu (); Lundholm (), spatial exponential decay of the interaction Rincon () and a supersymmetric version Guilarte (), along with mappings onto other models, such as attractive fermions Cheon (), a BCS model Fuchs (); Batchelor (); Wada (), the Kardar-Parisi-Zhang (KPZ) model LeDoussal (), directed polymers Luca (), three-dimensional black holes Panchenko () or the Yang-Mills problem on a 2-sphere Flassig ().

In this work, we focus on theoretical and mathematical issues associated with the analytical description of the ground state of the LL model. After so many years since the discovery of the closed-form system of equations by Lieb and Liniger, the explicit analytical solution of the model is still lacking. In particular, a wide range of experimentally relevant, intermediate repulsive interactions is hardly accessed analytically by perturbation theory. A complete analytical understanding of the ground state, whose importance is comparable to the solution of the 2D Ising model, would be a benchmark towards solutions at finite temperature or generalizations to multicomponent systems, where explicit results from integrability are scarce to date. The goal of this article is to provide analytical estimates for the ground state energy, excitation spectrum and other related observables covering the experimentally relevant regime of intermediate interaction strengths as defined in Eq. (3) below, whose accuracy would be comparable to the ones accessed by numerical methods. For this purpose, we start from the strongly interacting regime. We use a systematic method to obtain corrections to the TG limit up to high orders. Since this method is intrinsically limited to high interaction strength and the convergence with the order of expansion quite slow, we transpose the problem to number theory and guess a general, partially resummed structure. We study its accuracy by computing several observables and comparing the results with numerics, and approximations available in the literature. We conclude that our conjecture is valid in a considerably wider range of interactions than high-order asymptotic expansions. Then, we study the excitation spectrum, harder to access analytically. By a careful analysis of its symmetry properties, we find new analytical approximate expressions that involve quantites computed at equilibrium and compare them with numerics. We then find that the most accurate approximation is actually valid in a wide range of intermediate to strong interactions. Moreover, the accuracy of our methods enable us to study quantitatively the regime of validity of Tomonaga-Luttinger liquid theory in terms of the interaction strength. We show that the range of validity in momentum and energy increases with the coupling constant.

The structure of the paper is as follows: in Section II 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 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. In Section III, 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. In Section IV, we give our conclusions and outlook. Several Appendices dwell deeper on the mathematical details; an exact mapping onto the classical physics problem of the circular plate capacitor is also discussed.

## Ii Ground-state energy and local correlations from analyticity: methods and illustrations

### ii.1 Model and Bethe Ansatz equations in the thermodynamic limit

The 1D quantum gas composed of identical spinless point-like bosons of mass with contact interactions, confined to a line of length , is described by the Lieb-Liniger Hamiltonian which reads Lieb ()

 HLL=N∑i=1⎡⎣−ℏ22m∂2∂x2i+g\scriptsize{1D}2∑j≠iδ(xi−xj)⎤⎦, (1)

where label the positions of the bosons, is the Planck constant divided by , is an effective one-dimensional coupling constant which can be deduced from experimental parameters Olshanii (); Dunjko () and is the Dirac function. We assume , corresponding to repulsive interactions. In second-quantized form, the Hamiltonian reads Lieb ()

 HLL[^ψ]=ℏ22m∫L0dx∂^ψ†∂x∂^ψ∂x+g\scriptsize{1D}2∫L0dx^ψ†^ψ†^ψ^ψ, (2)

where is a bosonic field operator satisfying the canonical commutation relations with its Hermitian conjugate: , . The dimensionless coupling constant for this model, known as Lieb’s parameter, reads

 γ=mg\scriptsize{1D}n0ℏ2, (3)

where is the mean linear density. For a given atomic species, the coupling constant is experimentally tunable by confinement-induced resonances, Feshbach resonances or control over the density . The limit yields the Tonks-Girardeau gas. In this regime, an exact mapping on a noninteracting spinless Fermi gas allows for full solution of the model, with arbitrary external potential Girardeau ().

For a uniform 1D Bose gas 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 at arbitrary interaction strength. The procedure, called coordinate Bethe Ansatz, is well known for this system and widely explained in the literature, we refer to Lieb (); Zvonarev (); Panfil (); Franchini () for details. At finite , this yields a system of transcendental equations. In the thermodynamic limit, the latter reduces to a set of three equations, namely

 g(z;α)−12π∫1−1dy2αg(y;α)α2+(y−z)2=12π, (4)

where denotes the quasi-momentum distribution in reduced units, is the pseudo-momentum in reduced unit such that its maximal value is and is a non-negative parameter, in a one-to-one correspondence with the Lieb parameter introduced above via a second equation,

 γ∫1−1dyg(y;α)=α. (5)

The third equation yields the dimensionless average ground-state energy per particle , linked to the total energy by , according to

 e(γ)=∫1−1dyg(y;α(γ))y2[∫1−1dyg(y;α(γ))]3. (6)

Interestingly, Eq. (4) is decoupled from Eqs. (5)-(6), which is a specificity of the ground state YangYang (). Equation (4) is a homogeneous type II Fredholm integral equation with Lorentzian kernel Polyanin (). In the following we will refer to it as the Lieb equation for simplicity. Prior to Lieb and Liniger’s work, it had been observed that its solution yields the exact capacitance of a capacitor formed of two circular coaxial parallel plates. For more details on this interesting link with a standard problem of electrostatics, and approximation methods that are not considered in the main text, we refer the reader to Appendix A. We also refer to Kuestler () for an historical review of attempts to exactly solve this problem, that has resisted more than a century of efforts from mathematicians and physicists alike.

In order to determine the ground-state energy of the Lieb-Liniger model, the most crucial step is to solve the Lieb equation (4). This was done numerically by Lieb and Liniger for a few values of spanning several decades of interaction strengths Lieb (). The solution procedure relies on the following steps. An arbitrary positive value is fixed for , and Eq. (4) is solved, i.e. is found with the required accuracy as a function of . Then, Eq. (5) yields , subsequently inverted to get . In doing so, one notices that is an increasing function, thus interaction regimes are defined the same way for both variables. The energy is then easily obtained from Eq. (6), as well as many interesting observables, that are combinations of its derivatives. They all depend on the sole Lieb parameter, which is the key of the conceptual simplicity of the model.

In the following, we briefly tackle the weakly-interacting regime , mostly for the sake of completeness, and dwell much deeper on the strongly-interacting regime .

### ii.2 Expansions in the weakly-interacting regime

Both numerically and analytically, finding accurate approximate solutions of Eq. (4) at small values of the parameter is a more involved task than at higher couplings, due to the singularity of the function at , whose physical interpretation is that noninteracting bosons are not stable in 1D. A guess function was proposed by Lieb and Liniger Lieb (), namely

 g(z;α)≃α≪1√1−z22πα, (7)

which is a semi-circle law, rigorously derived in Hutson (). Heuristic arguments have suggested the following correction far from the edges in the variable Hutson (); Popov ():

 g(z;α)≃α≪1,α≪|1±z|√1−z22πα+14π2√1−z2[zln(1−z1+z)+ln(16πα)+1], (8)

reproduced later by direct calculation in Wadati () after regularization of divergent series. A technical difficulty consists in evaluating precisely the validity range of both approximations Eq. (7) and Eq. (8) in , and whether the latter offers a significant improvement in accuracy. Systematic comparison to numerical solutions at is beyond the scope of this work. One can nevertheless, to a certain extent, discuss the relative accuracy of two solutions at a given value of , using the methods detailed in Appendix B.

The highest-order exact expansion proven to date for the dimensionless energy in the weakly interacting regime is

 eexact(γ)=γ−43πγ3/2+[16−1π2]γ2−Dγ5/2+O(γ3). (9)

The first order coefficient is readily obtained from Eq. (7), the second order one is found using Eq. (8) and coincides with the result from Bogoliubov’s approximation Lieb (). The prefactor of the term has been more controversial and required important efforts. Its expression was inferred on numerical grounds in Takahashi (), while in Wadatibis () a factor instead of is found analytically after lengthy calculations. It seems, though, that the factor is correct, as it was recently recovered from potential theory in Widom () in a (quasi-)rigorous way, and coincides with numerical independent calculations carried out in Lee (). The exact fourth term, given by multiple integrals, was numerically evaluated to Leebis (). A similar value was then found by fitting using very accurate numerics Kardar ().

Our contribution to the problem is purely numerical in this regime. We evaluated the function numerically by solving the Bethe Ansatz equations with a Monte-Carlo algorithm for values of the interaction strength spanning the interval . A fit containing the known analytical prefactors and four adjustable coefficients yields

 efit(γ)=γ−43πγ3/2+[16−1π2]γ2−0.002005γ5/2+0.000419γ3−0.000284γ7/2+0.000031γ4 (10)

with a relative error as low as a few per thousands in the interval considered. We stress that in this expansion the various coefficients are not supposed to (and do not) coincide with the exact ones in the Taylor expansion. The expression yields, however, a rather good estimate for the ground-state energy from the noninteracting to the strong-coupling regime for many practical purposes. In particular, it will allow us to check whether the results of next section, where this regime is specifically addressed, are also valid at intermediate or even weak interactions.

### ii.3 Expansions in the strongly-interacting regime

Close to the Tonks-Girardeau regime, accurate asymptotic solutions of the Lieb equation Eq. (4) can be found using a double series expansion of the function in and around . Our starting point is a slightly simplified version of a recently-developed procedure Ristivojevic (), we refer to Appendix C for a detailed derivation and mathematical discussions. This systematic expansion scheme is valid for and yields an approximate solution of the form

 g(z;α,M)=2M+2∑k=0M∑j=0gjkz2jαk, (11)

where is an integer cutoff, and , by construction, are polynomials of with rational coefficients. This expansion coincides with the Taylor expansion at the chosen order, and converges to for .

We stress that the lowest interaction strength attainable with the strong coupling expansion, i.e. , corresponds to an interaction strength Ristivojevic (). This is low enough to combine with the result in the weakly-interacting regime, and obtain an accurate description of the ground state of the model over the whole range of repulsive interactions. Furthermore, this method yields several orders of perturbation theory at each step and is automatically consistent at all orders. Nonetheless, it is crucial to capture the correct behavior of as a function of in the whole interval to obtain accurate expressions, whereas the expansion converges slowly to the exact value at the extremities since the Taylor expansion in is done around the origin. This reflects in the fact that the maximum exponent of varies more slowly with than the one of . All in all, this drawback leads to a strong limitation of the validity range in at a given order, and calls for high order expansions, far beyond the maximal order treated so far, Ristivojevic (). For a detailed study of the accuracy of the method in terms of the cutoff , we refer to Appendix C. The high number of corrections needed seems at first redhibitory, but we noticed that doing the average over two consecutive orders in , denoted by

 gm(z;α,M)=g(z;α,M)+g(z;α,M−1)2, (12)

dramatically increases the accuracy, yielding an excellent agreement with numerical calculations for all at , as illustrated in Fig. 1.

However, at increasing the method quickly yields too unhandy expressions for the function , as it generates terms. This motivated us to seek compact representations and resummations for the function , allowing to easily use them for further applications and to generate them up to large orders. We present below an analysis of the structure of the terms entering the expansion (11) and a conjecture for compact expressions, yielding a partial resummation. Our conjecture has been then verified and validated by a very recent numerical approach Prolhac (); Prolhacbis ().

### ii.4 Conjectural expansions and resummations in the strongly-interacting regime

#### ii.4.1 Conjectures for g(z;α)

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. Terms of the first kind already arise at low orders in and , while terms of the second kind appear at higher orders and are expected to play a crucial role in the crossover region of intermediate interactions. These structures are conjectural, we infered them on as few first terms as possible and systematically checked that their predictions coincide with all higher-order available terms. While the simplest patterns are trivial to figure out, others are far more difficult to find because of their increasing structural complexity. Denoting by the ’kind’ and by an index for the elusive notion of ’complexity’, we write

 g(z;α,M)=12π+∑m,nIMm,n(z;α) (13)

where each is itself a double sum over all terms of given kind and complexity that appear at order . As an illustration, we detail the terms found to all orders up to included in Eq. (11).

The simplest term is

 IM1,0=1π2αM∑j=0(−1)j(zα)2j2(M−j)+1∑k=0(2πα)k. (14)

Terms of the first kind with complexity sum as

 IM1,1=−13π2α3M−1∑j=0(−1)j(zα)2j2(M−j)−1∑k=0(2πα)k[2k+(j+1)(2j+1)]. (15)

Terms of the first kind and complexity are

 IM1,2=145π3α6M−2∑j=0(−1)j(zα)2j2(M−j)−4∑k=0(2πα)k(20k2+ajk+bj), (16)

where and . Terms of the first kind and complexity are

 IM1,3=−12835π3α8M−3∑j=0(−1)j(zα)2j2(M−j)−6∑k=0(2πα)k(280k3+cjk2+djk+ej) (17)

where , and . The only term of the first kind and complexity we have identified is

 (18)

which should correspond to the terms with index in . Note that, for terms of first kind, complexity actually corresponds to the degree of the polynomial in .

Terms of the second kind and lowest complexity are

 IM2,0=1π2α5M−2∑j=0(−1)j(2j+5)α2j. (19)

We also found

 IM2,1=−121π2α5(zα)2M−3∑j=0(zα)2j(−1)jj+1M−j−3∑k=0(−1)kα2k(2(k+j+3)2j+1). (20)

We note that some terms of the first kind could be interpreted as second kind, thus involving a shift of the summation index in one of them, so that the proposed classification may not be optimal. However, if one performs the summation of all terms explicited here, the resulting expansion is exact up to order included, thus the proposed structures are efficient to encode many terms of the series in a relatively compact way. By comparing with high-precision numerical calculations, we find that for the expansion has an error of the order of . The problem of the full resummation of the series expansion for remains open.

#### ii.4.2 Conjectures for e(γ)

In this section, we focus on resummation patterns directly on the dimensionless ground-state energy . Here, we shall write , where once again the index denotes a notion of complexity. Focusing on the large- asymptotic expansion, we identify the pattern of a first sequence of terms. We conjecture that they appear at all orders and resum the series, obtaining

 e0(γ)eTG=+∞∑k=0(−1)k2k(k+11)γk=γ2(2+γ)2, (21)

where is the value in the Tonks-Girardeau regime. Then, using the expansion of in up to (given in Appendix E with numerical coefficients) and guided by the property

 +∞∑k=0(−1)k2k(k+3n+13n+1)γk=(γγ+2)3n+2, (22)

we conjecture that the structure of the term of complexity defined above is

 en(γ)eTG=π2nγ2Ln(γ)(2+γ)3n+2, (23)

where is a polynomial of degree , whose coefficients are rational, non-zero and of alternate signs. In this context, the notion of complexity is directly related to the power of the denominator. The first few polynomials are found as

 L1(γ)=3215, L2(γ)=−9635γ+848315, L3(γ)=512105γ2−4352525γ+131844725, L4(γ)=−102499γ3+1315845775γ2−4096275γ+117763465, L5(γ)=245761001γ4−2960506884729725γ3+4533678087882875γ2−2279444487882875γ+533377024212837625, L6(γ)=−409665γ5+614092835035γ4−469589196823648625γ3+371076300823648625γ2−1522810884729725γ+13433651242567525 (24)

by identification with the expansion to order . We conjecture that the coefficient of the highest-degree monomial of is .

Interestingly, contrary to the expansion, those partially resummed terms are not divergent at small , increasing the validity range. We also notice that corresponds to Lieb and Liniger’s approximate solution assuming a uniform density of pseudo-momenta Lieb (), and an equation equivalent to Eq. (21) appears in Astrabis (). The first correction was predicted rigorously in Wadatibis (), thus supporting our conjectures.

In the next section, we will evaluate the quality of our conjecture (23), (II.4.2) by comparing with high-precision numerical calculations.

### ii.5 High-precision ground-state properties

In this section we present our predictions for various physical quantities, using a combination of weak-coupling and strong-coupling expansion as well as the conjectures.

#### ii.5.1 Density of pseudo-momenta

The density of pseudo-momenta follows immediately from the solution of the Lieb equation (4). It is defined as

 ρ(k;γ)=g(Q(γ)z;α(γ)), (25)

where denote the pseudo-momenta, whose maximal value is , such that Zvonarev ()

 Q(γ)kF=1π∫1−1dzg(z;α(γ)), (26)

where is the Fermi wavevector in 1D. Within the method explained in Appendix C, we have access to analytical expressions for only. For a more general method to obtain this function, valid for any , we refer to Appendix D. Since the latter method is not appropriate to obtain analytical expressions of other quantities, we do not dwell further on it in the main text. Figure 2 shows our results for the density of pseudo-momenta. In the Tonks-Girardeau regime the density of pseudo-momenta coincides with the Fermi distribution at zero temperature, which is a manifestation of the effective fermionization. At decreasing interactions away from the Tonks-Girardeau gas, we find that the distribution remains quasi-uniform in a wide range of strong interactions. This shows the robustness of the effective Fermi-like structure, yet with Fermi wavevector which is progressively renormalized. Then, at lower interaction strengths, we witness an increase of the height of the peak of the distribution, which becomes progressively sharper and narrower around the origin in the weakly-interacting, quasi-condensate regime.

#### ii.5.2 Ground-state energy

We show in Fig. 3 the dimensionless ground-state energy per particle over a wide range of repulsive interactions. Our expressions in both the weakly-interacting regime as given by Eq. (10) and the strongly-interacting regime as given by Eq. (23) are compared to the numerics to emphasize their accuracy. As main result, we find that they have a wide overlap in the regime of intermediate interactions and are hard to distinguish from the numerics. In particular, extrapolating the conjecture in the strongly-interacting regime to low values of , we obtain a substantial improvement compared to the previously known approximations and in Eqs. (21) and (23) when .

We then show our results for the ratios of the mean kinetic energy and interaction energy per particle. According to Pauli’s theorem Lieb (), they are obtained as

 ep(γ)=γdedγ (27)

and

 ek(γ)=(e−γdedγ). (28)

These quantities are shown in the left panel of Fig. 4, normalized to the total energy in the Tonks-Girardeau regime as in Rigol (). 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 Girardeau (). The right panel of Fig. 4 shows the ratio of interaction to kinetic energy. This quantity scales as in the weakly-interacting regime and decreases monotonically at increasing , thus showing that does not represent this ratio, contrary to the mean-field prediction.

#### ii.5.3 Local correlation functions

In experiments, it is possible to access to the local -particle correlation functions of the Lieb-Liniger model, defined as

 gk=⟨[^ψ†(0)]k[^ψ(0)]k⟩nk0, (29)

where represents the ground-state average.

The local pair correlation (respectively three-body correlation ) is a measure of the probability of observing two (three) particles at the same position. In particular, governs the rates of inelastic processes, such as three-body recombination and photoassociation in pair collisions. The second-order correlation, , is easily obtained with our method using the Hellmann-Feynman theorem, that yields Gangardt (). Hence, the fact that is an increasing function is actually a direct consequence of the positiveness of . Higher-order correlation functions are related in a non-trivial way to the moments of the density of pseudo-momenta, defined as

 ϵk(γ)≡∫1−1dzzkg(z;α(γ))[∫1−1dzg(z;α(γ))]k+1. (30)

In particular, is related to the two first non-zero moments by the relation Cheianovqboson ()

 g3(γ)=32γdϵ4dγ−5ϵ4γ2+(1+γ2)dϵ2dγ−2ϵ2γ−3ϵ2γdϵ2dγ+9ϵ22γ2. (31)

By definition, the second moment coincides with the dimensionless energy . In Fig. 5 we plot accurate expressions for obtained in Cheianov (), and our analytical expression of readily obtained from Eqs. (10) and (23). The fact that vanishes in the Tonks-Girardeau regime is once again a consequence of fermionization, as interactions induce a kind of Pauli principle and preclude that two bosons come in contact. This property has been the key to realize the TG gas experimentally Gangardt (). At very small interaction strengths, however, the ratio can not be neglected; this means that the weakly-interacting 1D Bose gas is less stable with respect to three-body losses than the strongly-interacting one.

#### ii.5.4 Non-local, one-body correlation function and Tan’s contact

Finally, we study the one-body non-local correlation function . We focus first on the Tonks-Girardeau regime Lenard (); Vaidya (); Jimbo (), 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. We use the notation , where is the Fermi wavevector in 1D. We recall first the large-distance expansion derived in Jimbo ()(with signs of the coefficients corrected in Gangardtbis ())

 gTG1(z)=G(3/2)4√2|z|[1−132z2−cos(2z)8z2−316sin(2z)z3+3320481z4+93256cos(2z)z4+O(1z5)], (32)

where is the Barnes function defined as and the functional relation , being the Euler Gamma function.

At short distances, using the same technique as in Forrester () to solve the sixth Painlevé equation, we find the following expansion, where we added six orders compared to Jimbo ():

 gTG1(z)=8∑k=0(−1)kz2k(2k+1)!+|z|39π−11|z|51350π+61|z|7264600π+z824300π2−253|z|971442000π−163z1059535000π2 +7141|z|11207467568000π+589z126429780000π2−113623|z|13490868265888000π−2447503z141143664968600000π2 +(140186125000π3+3366129452095953280000π)|z|15+5597693140566821595200000π2z16+O(|z|17). (33)

The first sum is a truncation of the integer series defining the function , which corresponds to the one-body correlation function of noninteracting fermions. The additional terms appearing in this expansion, and in particular the odd ones, are peculiar of bosons with contact interactions. Actually, the one-body correlation function for Tonks-Girardeau bosons differs from the one of a Fermi gas due to the fact that it is a nonlocal observable, depending also on the phase of the wavefunction and not only on its modulus.

The same structure is valid at finite interaction strength, where the short-distance expansion reads

 g1(z)=1++∞∑i=1ciπi|z|i (34)

and the first coefficients are explicitly found as Dunn ()

 c1=0, c2=−12ek, c3=112γ2dedγ, (35)
 c4=γ12dϵ4dγ−3ϵ48+2γ2+γ324dedγ−γe6−γ4ededγ+34e2. (36)

Our solution of the Lieb equation Eq. (4) hence allows to estimate the first terms of this expansion. Further progress on analytic expressions has been obtained in Jimbobis (); Dunn (), as well as numerically Astra (); Slavnov ().

The Fourier transform of the one-body correlation function is the momentum distribution of the gas,

 n(k)=n0∫+∞−∞dxg1(x)e−ikx. (37)

The short-distance behavior of in Eq. (34) allows to obtain the large momenta asymptotic behavior,

 k4n(k)n40→k→+∞C, (38)

where Dunn (); Olsh () is Tan’s contact Tan (). Figure 7 shows the value of Tan’s contact obtained from Eqs. (10) and (23)-(II.4.2).

## Iii Excitation spectrum, exact and approximate results

Now that the ground-state energy of the model is known with good accuracy, we proceed with a more complicated and partially open problem, namely the analytical characterization of excitations above the ground state at zero temperature. In particular, we are interested in the accuracy of field-theoretical approximations such as Luttinger liquid theory. In order to introduce these excitations, we consider the Bethe Ansatz solution at finite . The total momentum and energy of the system are given by and respectively LiebII (), where the set of quasi-momenta satisfies the following system of transcendental equations

 {kj=2πLIj−2LN∑l=1arctan(kj−klγn0)}j∈{−N−12,…,N−12}. (39)

The ’s, called Bethe quantum numbers, are integer for odd values of and half-odd if is even. Since we consider , the quasi-momenta are real and can be ordered in such a way that . Then, automatically LiebII (). The ground state corresponds to and has total momentum at arbitrary interaction strength. In what follows, we use the notations and for the total momentum and energy of an excitation with respect to the ground state, so that the excitation spectrum is defined as . For symmetry reasons, we will only consider excitations such that , those with having the same energy.

The Lieb-Liniger model features two excitation branches, denoted as type I and type II LiebII (). In order to explain their features, we derive them in the Tonks-Girardeau limit GirardeauMinguzzi (), where the set of equations (39) decouples and . The type-I excitations corresponding to the highest energy excitation, occur when the highest-energy particle with gains a momentum and an energy . The corresponding dispersion relation is

 ϵI(p)=12m[2pFp(1−1N)+p2] (40)

where is the Fermi momentum.

The type-II excitations occur when a particle inside the Fermi sphere is excited to occupy the lowest energy state available, carrying a momentum . This type of excitation corresponds to shifting all the rapidities with by , thus leaving a hole in the Fermi sea. This corresponds to an excitation energy , yielding the excitation branch

 ϵII(p)=12m[2pFp(1+1N)−p2]. (41)

Any combination of one-particle and one-hole excitation is possible, giving rise to an intermediate energy between and , forming a continuum in the thermodynamic limit. Figure 8 shows the type-I and type-II excitation spectrum for bosons in the Tonks-Girardeau limit. We notice the symmetry , valid at large boson number for the type-II branch.

In the general case of finite interaction strength, the system of equations (39) for the many-body problem can not be solved analytically, although expansions in the interaction strength can be obtained in the weakly- and strongly-interacting regimes Oelkers (). The solution is easily obtained numerically for a few bosons. To reach the thermodynamic limit with several digits accuracy the Tonks-Girardeau treatment suggests that should be of the order of , and the interplay between interactions and finite may slow down the convergence at finite Calabrese (), but a numerical treatment is still possible. Here, we shall directly address the problem in the thermodynamic limit, where it reduces to two equations Ristivojevic (); Pustilnik () :

 p(k;γ)=2πℏQ(γ)∣∣∣∫k/Q(γ)1dyg(y;α(γ))∣∣∣ (42)

and

 ϵ(k;γ)=ℏ2Q2(γ)m∣∣∣∫k/Q(γ)1dyf(y;α(γ))∣∣∣, (43)

where, according to Eq. (26), . It is known as the Fermi rapidity, represents the radius of the quasi-Fermi sphere and equals in the Tonks-Girardeau regime. The function satisfies the integral equation

 f(z;α)−1π∫1−1dyαα2+(y−z)2f(y;α)=z, (44)

referred to as the second Lieb equation in what follows. We solve it with similar techniques as for the Lieb equation (4). Details are given in Appendix F.

The excitation spectra at a given interaction strength are obtained in a parametric way as . Within this representation, one can interpret the type I and type II spectra as a single curve, where the type I part corresponds to and thus to quasi-particle excitations, while type II is obtained for , thus from processes taking place inside the quasi-Fermi sphere, which confirms that they correspond to quasi-hole excitations. Using basic algebra on Eqs. (42) and (43) we obtain the following interesting general results:

(i) The ground state ) trivially corresponds to , showing that represents the edge of the Fermi surface.

(ii) The quasimomentum corresponds to the umklapp point , always reached by the type II spectrum in the thermodynamic limit, regardless of the value of .

(iii) The maximal excitation energy associated with the type II curve lies at , corresponding to .

(iv) If , and , hence , generalizing to finite interactions the symmetry found in the Tonks-Girardeau regime.

(v) The type I curve repeats itself, starting from the umklapp point, shifted by in . Thus, what is usually 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, . 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 explicitly, with excellent accuracy provided . As far as the type I curve is concerned, however, we are not only limited by , but also by the fact that our approximate expressions for and are valid only if , , thus adding the validity condition, . The latter is not very constraining as long as , but for the best validity range we can get is very narrow around . To bypass this problem, Ristivojevic used an iteration method to evaluate and Ristivojevic (). However, in practice this method is applicable only for large interactions since it allows to recover only the first few terms of the exact expansion of and (to order in Ristivojevic ()). Moreover, the obtained expressions are not of polynomial type, it is then a huge challenge to substitute the variable to express explicitly, and one has to use approximate expressions at high and small momenta.

In the regime , we first compute the M-th order mean approximant for the function , defined in Eq. (12),

 gm(z>1;α,M)=12π+απ∫1−1dygm(y;α,M)α2+(y−z)2, (45)

which then allows us to obtain the type I spectrum with excellent accuracy for all values of from Eqs. (42) and (43). Both excitation spectra are shown in Fig. 9 for several values of . We note that the area below the type II spectrum, vanishing in the noninteracting Bose gas, is an increasing function of .

At small momenta, in the general case, due to the analytical properties of both and , for all values of the type I curve can be expressed as a series in Furusaki (); RistiMat (),

 ϵI(p;γ)=vs(γ)p+p22m∗(γ)+λ∗(γ)6p3+…. (46)

In the Tonks-Girardeau regime, as follows from Eq. (40), one has , , and all other coefficients vanish. At finite interaction strength, the parameters and can be seen as a renormalized Fermi velocity and mass respectively. Linear Luttinger liquid theory predicts that is the sound velocity associated with bosonic modes at very low Haldane (). To all accessed orders in , we have checked that our expression agrees with the exact thermodynamic equality LiebII ()

 vs(γ)=vFπ[3e(γ)−2γdedγ(γ)+12γ2d2edγ2(γ)]1/2. (47)

Analytical expansions for the sound velocity are already known at large and small interaction strengths. The first- and second-order corrections to the Tonks-Girardeau regime in are given in Cazalilla (), they are calculated to fourth order in Zvonarev () and up to eighth order in Ristivojevic (). In the weakly-interacting regime, expansions are found in Cazalilla () and Wadati (). Our expressions Eqs. (23, II.4.2) for the dimensionless energy density allow us to considerably increase the accuracy compared to previous works after easy algebra. Interestingly, one does not need to compute the ground-state energy to find . It is sufficient to know the function at in the first Lieb equation due to the useful equality Korepin ()

 vs(γ)vF=1[2πg(1;α(γ))]2. (48)

Reciprocally, knowing yields an excellent accuracy test for the function, since it allows to check its value at the border of the interval where it is the 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 the left panel of Fig. 10. One can see that , thus we shall find , which shows that the method of polynomial expansion must fail 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 .

Linear Luttinger liquid theory assumes that, for all values of , and . This strictly linear spectrum is however a low-energy approximation, and nonlinearities cannot be neglected if one wants to deal with higher energies. Here, we provide the first quantitative study of the regime of validity of the linear approximation at finite interaction. We denote by and respectively the half-width of momentum, and maximum energy range around the umklapp point such that the linearized spectrum is exact up to percent. These quantities should be considered as upper bounds of validity for dynamical observables such as the dynamic structure factor Lang1 (). Our results for and are shown in Fig. 11. One sees that Luttinger liquid theory works well at large interaction strengths. However, its range of validity decreases when interactions are decreased from the Tonks-Girardeau regime.

Including the quadratic term and neglecting higher-order ones is actually a complete change of paradigm, from massless bosonic to massive fermionic excitations at low energy, at the basis of the Imambekov-Glazman theory of nonlinear Luttinger liquids Glazmanbis (); Glazmanter (). In this approach is interpreted as an effective mass, whose general expression is Ristivojevic (); MatPustilnik ()

 mm∗=(1−γddγ)√vsvF. (49)

We have verified that the effective mass and the sound velocity obtained from the excitation spectrum satisfy Eq. (49) to all accessed orders in . Our results for the effective mass as obtained from a combination of the weak-coupling expression Eq. (10), the conjecture (23), and the use of Eqs. (47) and (49), is shown in the right panel of Fig. 10. We notice that the inverse effective mass vanishes for . This is also predicted by the Bogoliubov theory, where the small- expansion of the excitation dispersion reads . Hence, the non-vanishing inverse effective mass is a beyond-mean field effect.

For the type II spectrum, the properties (i)-(v) that we have detailed above suggest another type of expansion, also used in Brand ():

 ϵII(p;γ)=12m[f1(γ)p(2pF−p)+f2(γ)p2(2pF−p)2p2F+…]. (50)

Using the property (vi), on the other hand, allows us to write

 ϵII(p;γ)=vs(γ)p−p22m∗(γ)+… (51)

Equating both expressions to order , one finds that

 f1(γ)=vs(γ)vF. (52)

A similar result was recently inferred from Monte-Carlo treatment of 1D He and proved by Bethe Ansatz applied to the hard-rods model in Bertaina (); Bertainabis (). By the same approach we then show that

 f2(γ)=14(vs(γ)vF−mm∗(γ)), (53)

which is one of our main new results.

Systematic suppression of the variable and higher-order expansions suggest that, at large enough values of at least, higher-order terms in expansion (50) can be neglected. Figure 12 shows the value of the Lieb-II excitation spectrum at its local maximum value , as obtained from a numerical calculation as well as the expansion (50). We find that the result to order one is satisfying at large , but the second order correction significantly improves the result at intermediate values of the Lieb parameter. Our numerical calculations show that third and higher order corrections are negligeable in a wide range of strong interactions. Overall, we dispose of very accurate analytical predictions for the excitation branches, for .

## Iv Conclusions and outlook

In conclusion, first we have solved with high accuracy the set of Bethe-Ansatz equations Eqs. (4), (5) and (6) 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 of two simple analytical expressions which describe to a good accuracy the ground state energy, namely, a weak coupling expansion, Eq. (10), valid for interaction strengths , and a strong-coupling expansion to order , whose coefficients are given numerically in Eq. (E), valid for interaction strengths . Their combination spans the whole range of coupling constants, thus providing an alternative to the use of tabulated values and an opportunity to accurately benchmark numerical methods Cirac (); Ganahl ().

<