Conjectures about the structure of strong-and weak-coupling expansions of a few ground-state observables in the Lieb-Liniger and Yang-Gaudin models

In this paper, we apply experimental number theory to two integrable quantum models in one dimension, the Lieb-Liniger Bose gas and the Yang-Gaudin Fermi gas with contact interactions. We identify patterns in weak- and strong-coupling series expansions of the ground-state energy, local correlation functions and pressure. Based on the most accurate data available in the literature, we make a few conjectures about their mathematical structure and extrapolate to higher orders.


I. INTRODUCTION
The Lieb-Liniger model describes spinless bosons with contact interactions, whose motion is confined to one dimension [1]. It is arguably the most conceptually simple, integrable model in the continuum. Each of its observables is amenable to exact calculations using Bethe Ansatz techniques [1, 2], providing valuable insights in the properties of strongly-correlated systems in one dimension [3]. The accuracy of the ground-state solution can be used to benchmark effective theories, such as the Luttinger liquid approximation [4][5][6] and its non-linear generalization [7,8], as well as numerical techniques [9][10][11][12]. Last, but not least, this model can be experimentally realized in ultracold gases experiments, both in the repulsive [13][14][15][16][17][18] and attractive regime [19].
A key observable is the ground-state energy, which is related to many other physical quantities through theorems and thermodynamic relations. For instance, in the Lieb-Liniger model, the local two-body correlation function can be extracted from the latter [36], as well as the Luttinger coefficients [5] that govern the long-range behaviour of non-local correlation functions, and Tan's contact, that yields the high-momentum tail of the momentum distribution [37][38][39]. Most remarkably, even the excitation spectrum can be obtained [40], by calculating the effective mass [41] and higher-order coefficients of its low-momentum expansion [42]. This spectrum, in turn, is related to the dynamical structure factor [43].
Although the exact ground-state energy can be obtained, in principle, from the exact Bethe Ansatz equations, only weak-and strong-coupling expansions are accessible to date [1, 40,[44][45][46][47][48][49][50][51]. As far as the Lieb-Liniger model is concerned, recently-developed algorithms allow to obtain these expansions to any order [40,49,50]. An important step forward would be to understand them well enough to predict their coefficients at arbitrary order without evaluating them algorithmically. This paper is organized as follows. In section II, we briefly describe the Lieb-Liniger and Yang-Gaudin models, introduce a few notations and the Bethe Ansatz equations that give access to the ground-state energy. In section III, we sum up the main known results about their strong-and weak-coupling expansions. We then propose new conjectures about the ground-state energy of the Yang-Gaudin model, and discuss the radius of convergence of the partial resummations we perform. In section IV, we sum up the main analytical developments related to the local correlation functions of the Lieb-Liniger model and propose a conjecture about their structure. Then, in section V, we link together two independent descriptions of the local two-body correlation function to evaluate an observable that we identify as the pressure of the Bose gas. We provide conjectures about its strong-coupling expansion. In section VI, we conclude and give a few outlooks.

II. MODELS AND EQUATIONS
The one-dimensional quantum gas composed of N identical spinless point-like bosons of mass m with contact interactions, is described by the Lieb-Liniger Hamiltonian H which reads [1] 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 and δ is the Dirac function. It is customary to introduce a dimensionless coupling constant, known as Lieb's parameter, that reads where n 0 is the average linear density.
The coordinate Bethe Ansatz procedure allowing to solve this model is well-known, see e.g. [52] for details. It yields a set of three equations, namely where g(z; α) denotes the distribution of quasi-momenta in reduced units, z is the pseudo-momentum in reduced units such that its maximal value is 1 and α is a real number in one-to-one correspondence with the Lieb parameter γ introduced above through a second equation, A third equation yields the dimensionless average ground-state energy per particle, e B (γ), linked to the total energy E 0 by according to A similar procedure allows to obtain the ground-state energy of the Yang-Gaudin model of spin-1/2 fermions with contact interactions, and even the κ-component model with an arbitrary number of spin components, see e.g. [53,54] for details. The Hamiltonian stays the same as in Eq. (1), but the statistics is modified to take into account the fermionic nature of the particles. In the case of a κ-component Fermi gas divided into [m 1 , . . . , m κ ] particles per species, a set of κ rapidities, distributions and integrations bounds denoted respectively k i , ρ i (k i ) and B i can be defined. Introducing the notation M i = κ j=i m j for all i ∈ {1, . . . , κ}, the κ coupled integral Bethe Ansatz equations for the ground state in the thermodynamic limit can be written as: together with the following normalization condition: The dimensionless average ground-state energy per particle reads:

III. GROUND-STATE ENERGY
Although solving the Lieb equations yields, in principle, the exact ground-state energy, only weak and strongcoupling expansions are known to date. In this section, we recap the most accurate results available and take a close look at their patterns to try and identify their generating series.

A. Strong-coupling expansion
Neumann series expansion of the kernel is the standard mathematical tool to solve integral equations like Eq. (3), see e.g. [46]. However, only the two first orders of this series expansion have been derived analytically [40], due to the increasing complexity of higher-order terms.
An algorithm that yields systematic series expansions at large Lieb parameter γ has been developed in [40]. It allows to evaluate analytically the ground-state energy of the Lieb-Liniger model in the vicinity of the Tonks-Girardeau regime γ → +∞ as a truncated series in the inverse Lieb parameter, We have pushed this procedure up to order n = 20 in [55], the explicit result is given here in appendix A. More generally, at order n the expression contains 1 + ⌊(n + 1)/2⌋⌊1 + n/2⌋ different terms, where ⌊. . . ⌋ is the floor function, such that ⌊x⌋ is the nearest lower integer to x, or x if the latter is integer. For this reason, explicit expressions of a k are rather uncanny for large n, all the more so as the coefficients become more and more complex. Moreover, the convergence of the series to the exact value at given γ is quite slow at intermediate coupling, and the procedure only applies when α > 2 in Eq. (3), corresponding to γ > 4.527.
The will to bypass these shortcomings led us to look at the problem from the perspective of experimental number theory, trying to identify patterns in the expansion to express it with more compact notations. This heuristic point of view turned out to be rather fruitful. We guessed a structure, seemingly valid at all orders, whose partial resummation led to the following pattern: where and for n ≥ 1, where E n (γ) = n−1 k=0 a k;n γ k is a polynomial of degree n−1 with rational coefficients. We also guessed a few properties of the latter. In particular, we identified the coefficient of the the highest-degree monomial of E n [55]. This conjecture generalizes a result of [56] to arbitrary order, and is in agreement with the exact series expansion to order 20. The role of the factor 1 + 2 γ that appears at the denominator when factorizing by γ is discussed in an appendix of [57]. Other lines could be followed to seek structures, see Appendix B for an example.
The validity of these assumptions has been checked at higher orders by independent numerical calculations [58]. Identifying the structure of coefficients in the polynomials E n allows to make more accurate calculations [55], compared to the mere series expansion in 1/γ. Indeed, a major advantage of this partially resummed expansion is that its radius of convergence is infinite. A closed-form formula for the coefficients that would lead to a full resummation would solve the ground-state energy problem explicitly.
In comparison, relatively few results have been obtained concerning the Yang-Gaudin and general κ-component model, due to the much more complex structure of the equations involved. In the balanced gas where all spin sectors contain the same number of fermions, the most accurate strong-coupling expansion of the average ground-state energy per particle available to date reads [60] where , with ψ the Euler psi function, C the Euler constant and ζ the Riemann zeta function.

Conjecture I:
The resummation above, Eq. (13), suggests a similar structure for e B (γ) and e F (γ; κ) through the Yang-You theorem as: where is inferred from the third-order strong-coupling expansion above. It might be that 2Z 1 (κ) + γ plays the role of 2 + γ at the denominator to all orders, and that linear combinations of 1 and terms of the form Z 2j+1 (κ) appear as prefactors at the numerator, but too few orders are available to date to make reliable conjectures about the higher-order structures e F,n when n ≥ 1.

B. Weak-coupling expansion
In the weakly-interacting regime, expansions of the ground-state energy of the Lieb-Liniger model have the following structure [46]: The zeroth-and first-order coefficients have been obtained in [1]. The second-order term was inferred in [44] and checked wih more rigor in [45,46]. Coefficients of order 3 to 5 were inferred from a very accurate numerical study [47,48], and involve the zeta function evaluated at odd arguments, ζ(3) and ζ(5). Very recently, further coefficients have been explicitly obtained using two independent methods. On the one hand, an efficient algorithmic expansion at low γ, combined to an integer coefficients algorithm to find the explicit form of the latter has been used in [49] and allowed to identify coefficients up to order 8. On the other hand, a systematic algorithmic procedure that yields directly the exact coefficients has been developed in [50], and has yielded their value up to order 34, though they have been published up to eighth order only due to their lengthy expressions. These studies revealed that from a certain order on, c B k involves products of the zeta function evaluated at odd arguments. The asymptotic expression of c B k for k ≫ 1 has even been inferred from numerical data [50]. As far as the Yang-Gaudin model is concerned, the weak-coupling expansion of the ground-state energy reads For a long time, only the first three coefficients were known [62]. A few more have been obtained numerically in [47,48], and the fourth has been computed explicitly [51]. A systematic procedure that yields the exact coefficients has been developed [50,51], and has given the explicit form of the coefficients up to order 34. They are given explicitly up to order 10 in [51].
We have carefully studied these new available data and found out structures that may be valid at arbitrary order. We thus make a few conjectures about subseries in the weak-coupling expansion of the Yang-Gaudin model.

Conjecture II:
We guess that ζ(3) appears as a global factor in a peculiar subsequence of the weak-coupling expansion of the ground-state energy of the Yang-Gaudin model, as: This result is in full agreement with [51] to all published orders. It has been confirmed up to order 50 by the authors of the aforementioned article [59]. This subseries converges provided that γ ≤ π 2 /2, and sums up to Conjecture III: We also point out that ζ(5) appears as a global factor in a peculiar subsequence of the weak-coupling expansion of the ground-state energy of the Yang-Gaudin model, that should read: This result is in full agreement with [51] to all published orders. As well as the former conjecture, this one has been confirmed up to order 50 by the authors of the aforementioned article [59]. The series converge provided that γ ≤ π 2 /2, and sums up to The limited radius of convergence of the series expansion in this conjecture and the previous one is due to the pole of the denominator in the functions involved in their resummation. At least one of the other subseries, most liketly the one that does not have any odd zeta prefactor, must have a radius of convergence equal to zero according to [51], where it was shown that the series on its whole has a null radius of convergence.

Conjecture IV:
We noticed that the weak-coupling expansions of the ground-state energy in the case of the Yang-Gaudin and Lieb-Liniger model, along with the Yang-You theorem, suggest a general structure for the weak-coupling expansion of the ground-state energy of the κ-component model that reads: wherec F k is a linear, rational combination of 1, and products of zeta function evaluated at odd arguments, except forc F 2 which also contains ζ(2). The fact that ζ(2) and more generally ζ(2m) is specific toc F 2 has been checked to low orders, and implemented by T. Reis and M. Marino in their calculations at very high orders to achieve workable computation times [59]. Equation (25) is a common generalization of Eqs. (19,20).
In the balanced case, i.e. when the number of particles in each component is linked to the total number of particles by N i = N/κ, it is known that [60] e F (γ; κ) = π 2 3κ 2 + κ − 1 κ γ + . . .

IV. LOCAL CORRELATION FUNCTIONS OF THE LIEB-LINIGER BOSE GAS
In the Lieb-Liniger model, the M-body local correlation functions are usually defined as where ψ(x) is the annihilation operator of the bosons. These observables indicate the probability of observing M particles at the same location at the same time. In a standard framework, their computation involves moments of the density of pseudo-momenta [36,63], , through the concept of connection, defined as a relationship between a local correlation function and coefficients of the short-distance expansion of the one-body correlation function [64]. This approach leads to rather complicated expressions for g M in terms of e 2k and their derivatives, already at order M = 3 [65]. A few properies of the moments are given in appendix C. Based on strong-coupling expansions of g 2 and g 3 to order 20 in the inverse coupling constant, and inspired by the resummation of e(γ) as in Eq. (13), we propose a conjecture about the general structure of the M -body correlation function.

Conjecture V:
We guess that the structure of the strong-coupling expansion of g M (γ) takes the form: It is in agreement with our expansions of g 2 and g 3 to order 20, and with the most accurate explicit result to date [66] up to order M 2 − M + 1 in 1/γ for all values of M , in particular through the known expression of G 0,M . It also agrees with the exact result g 1 (γ) = 1, on the condition that G n;0 = 0 for n ≥ 1.
Moreover, we have spotted out that G n;2 = n k=0 b k;n γ k with b 0;n = 4a 0;n and b n;n = −(2n + 1)a n−1;n with the notations involved in Eq. (13) and below.

V. PRESSURE OF THE LIEB-LINIGER BOSE GAS
Using an other formalism, the local correlation functions of the Lieb-Liniger model can be expressed in terms of solutions of other integral equations than the Lieb equation Eq. (3). For instance, the two-body local correlation function g 2 can be written as [67] where by definition The function g(z; α) has been defined above in Eq. (3), and f (z; α) is solution to the integral equation We combine several approaches to identify the physical meaning of p(γ). It is a well-known fact that the Hellmann-Feynman theorem yields e ′ (γ) = g 2 , where ′ denotes the derivation with respect to γ [36], thus On the other hand, the pressure P satisfies the thermodynamics relation [68] Therefore, we find out that p = e − γ 2 e ′ = P 2 n 2 0 /m is actually the dimensionless pressure of the Bose gas, an observable not quite studied so far in the case of the Lieb-Liniger model. More thermodynamic relations are derived in Appendix D. Our calculations of the strong-coupling expansion of the ground-state energy allowed us to obtain an expression for the pressure at the same order in 1/γ, whose more general structure can be identified as in the case of the ground-state energy.

Conjecture VI:
Based on our strong-coupling expansion of p(γ) at order 20 in the inverse coupling constant, we conjecture that it can be partially resummed as: where P n are polynomials of degree max(n−1, 0), with rational coefficients of alternating signs. We have identified the first few polynomials as: We found this structure using the property: and in analogy with the structure of the ground-state energy, Eq. (15). Moreover, the highest-degree monomial of P k , where k ≥ 1, seems to be affected of a coefficient According to Eq. (29), the ground-state energy e(γ) can be obtained from a differential equation knowing the pressure p(γ), and the coefficients of their series expansion at weak and strong coupling are linked together. At weak coupling, the coefficient of γ 2 in the series expansion of the pressure is null, so the corresponding coefficient in the series expansion of e B (γ), that involves ζ(2), is actually an integration constant. This might explain why ζ(2) is specific to this order and does not appear at higher orders anymore. The differential equation also gives an explanation to the fact that products of zeta functions, such as ζ(3) 2 , appear in the expansion.

VI. CONCLUSION AND OUTLOOK
In conclusion, in this article we used the most accurate available data to guess patterns in the weak-coupling series expansion of the ground-state energy of the Yang-Gaudin model. We proposed two conjectures about the mathematical structure of subseries whose terms share a common prefactor that involves ζ(3) and ζ(5) respectively. Although they encapsulate an infinite number of terms, they do not encompass multiples of zeta functions evaluated at odd arguments, that appear at higher orders. Lack of data did not allow us to generalize our conjectures to all odd arguments of the zeta function, but the expressions already obtained may serve as a guideline to identify more coefficients in the future. The subseries we have resummed have a finite radius of convergence, and are not responsible for the null radius of convergence of the full series, which should be understood from the subseries devoid of odd zeta terms.
All the results obtained for the ground-state energy of the Yang-Gaudin model, e F (γ; κ = 2), readily yield their equivalent in the super-Tonks-Girardeau regime of the Lieb-Liniger model, through the relation e ST G B (γ) = 16e F |γ| 4 . Surprisingly, the coefficients of the standard Lieb-Liniger model with repulsive interactions seem more difficult to guess than those of the Yang-Gaudin model.
We also made conjectures about the structures of the weak-coupling and strong-coupling series expansion of the ground-state energy of κ-component fermions. The fact that almost all weak-coupling coefficients are rational linear combinations of 1 and zeta function evaluated at odd arguments should be related to the model's integrability [59]. The only exception here concerns the coefficient of γ 2 , which might be a special feature of the κ-component model in general. Its role of an integration constant may be a clue to this special behavior. We note that a similar structure has been observed in the emptiness formation probability of the XXX spin chain [69], where the link between the appearance of the zeta function at odd arguments and integrability has been investigated and a similar conjecture has been stated. In this case, however, factors of ln(2) also appear, that are absent in the weak-coupling expansions of the ground-state energy of the Lieb-Liniger and Yang-Gaudin models. However, both factors ln(2) and ζ(2j + 1) can be encompassed in the more general framework of polylogarithms. Understanding at a deeper level the link between number theory and integrability is a thrilling problem in mathematical physics.
Then, we devoted our attention to the local correlation functions of the Lieb-Liniger model. Inspired by a conjecture about the strong-coupling expansion of the ground-state energy of the Lieb-Liniger model, and using series expansions of the two-and three-body correlation functions, we proposed yet another conjecture about the structure of general local correlation functions, that encompasses the most accurate analytical results available in the literature as special cases.
Finally, we combined two approaches to evaluate the two-body correlation function, to identify another groundstate observable as the pressure of the Lieb-Liniger Bose gas. We proposed a conjecture about the structure of its strong-coupling expansion.
As an outlook, it would be interesting to systematically compare several approaches to the local correlation functions. In [64], we pointed out the existence of a general framework through the notion of connection, that would deserve more investigations. In parallel, an other very interesting formalism exists, whose solution requires to solve several integral equations. It is thus more demanding in terms of computations, but better adapted to numerical studies as it does not involve derivatives. Its structure is simpler, and explicitly known at all orders M [70], it is valid at finite temperature and can be adapted to out-of-equilibrium situations, yet actual calculations of series expansions in terms of the coupling constant are still lacking.
As far as non-local correlation functions of the Lieb-Liniger model are concerned, they have been studied extensively in the classical textbook [71]. Working out explicitly their dependence on the coupling constant and temperature is still a partially-open problem.  To seek the structure of the ground-state energy, in this appendix we rely on the fact that any function can be decomposed in a unique way as a sum of an even and an odd function. The dimensionless ground-state energy is thus split into e B (γ) = e even (γ) + e odd (γ). (B1) We try and identify patterns of increasing complexity in e even and e odd , writing e even/odd (γ) = +∞ n=0 e even/odd,n (γ).
In doing so, we found involved in the calculation of the moments of the density of pseudo-momenta, Eq. (28).
Appendix D: More relationships for the pressure The dimensionless pressure p satisfies It is known that the dimensionless chemical potential satisfies [1] Combining these expressions yields, in real units, Also, introducing the operators and we notice that and µ = F (e).
Since F and G commute, and thus