Exact perturbative results for the Lieb-Liniger and Gaudin-Yang models

TL;DR

This paper develops a systematic method to derive and analyze the divergent perturbative series for the ground state energies of the Lieb-Liniger and Gaudin-Yang models, revealing their factorial divergence and connection to non-perturbative effects.

Contribution

It introduces a procedure to explicitly compute perturbative series coefficients from Bethe ansatz solutions for these models, and studies their large order behavior and non-perturbative implications.

Findings

Both series diverge factorially and are not Borel summable.

The first Borel singularity in the Gaudin-Yang model relates to the non-perturbative energy gap.

Provides insights into the Cooper instability through perturbative series analysis.

Abstract

We present a systematic procedure to extract the perturbative series for the ground state energy density in the Lieb-Liniger and Gaudin-Yang models, starting from the Bethe ansatz solution. This makes it possible to calculate explicitly the coefficients of these series and to study their large order behavior. We find that both series diverge factorially and are not Borel summable. In the case of the Gaudin-Yang model, the first Borel singularity is determined by the non-perturbative energy gap. This provides a new perspective on the Cooper instability.

Equations

$$H_N=-\sum _{j=1}^N\frac{{\partial }^2}{\partial x_j^2}+2c\sum _{1\le i<j\le N}\delta (x_i-x_j).$$

$$\frac{f_{\mathrm{B}}(x)}{2}-\frac{1}{2\pi }{\int }_{-B}^B\frac{f_{\mathrm{B}}(y)\mathrm{d}y}{(x-y{)}^2+1}=1.$$

$$\frac{1}{\gamma }=\frac{1}{4\pi }{\int }_{-B}^B{f}_{\mathrm{B}}(x)\mathrm{d}x,$$

$$e_B(\gamma )=\frac{E}{n^3}=\frac{{\gamma }^3}{4\pi }{\int }_{-B}^B{x}^2{f}_{\mathrm{B}}(x)\mathrm{d}x.$$

$$\frac{f_{\mathrm{F}}(x)}{2}+\frac{1}{2\pi }{\int }_{-B}^B\frac{f_{\mathrm{F}}(y)\mathrm{d}y}{(x-y{)}^2+1}=1.$$

$$\frac{1}{\gamma }=\frac{1}{\pi }{\int }_{-B}^B{f}_{\mathrm{F}}(x)\mathrm{d}x,$$

$$e_{\mathrm{F}}(\gamma )=\frac{E}{n^3}=-\frac{{\gamma }^2}{4}+\frac{{\gamma }^3}{\pi }{\int }_{-B}^B{x}^2{f}_{\mathrm{F}}(x)\mathrm{d}x.$$

$$e_{\mathrm{B}}(\gamma )$$

$$e_{\mathrm{F}}(\gamma )$$

$$R(x)={\int }_{-B}^B\frac{f(x^{\prime })}{x-x^{\prime }}\mathrm{d}{x}^{\prime }.$$

$$f(x)=-\frac{1}{2\pi \mathrm{i}}\left(R^{+}(x)-R^{-}(x)\right),$$

$$R(x)=\sum _{k\ge 0}⟨x^k⟩x^{-k-1},⟨x^k⟩={\int }_{-B}^{B}f(x)x^k\mathrm{d}x.$$

$$\frac{1}{x^2+1}=-\frac{1}{2\mathrm{i}}\left(\mathsf{D}-{\mathsf{D}}^{-1}\right)\frac{1}{x},\mathsf{D}={\mathrm{e}}^{\mathrm{i}{\partial }_x},$$

$$(1-\mathsf{D}){\tilde{R}}_{\mathrm{B}}^{+}(x)-(1-{\mathsf{D}}^{-1}){\tilde{R}}_{\mathrm{B}}^{-}(x)$$

$$(1+\mathsf{D})R_{\mathrm{F}}^{+}(x)-(1+{\mathsf{D}}^{-1})R_{\mathrm{F}}^{-}(x)$$

$$B,x\to \infty ,u=\frac{x}{B}\text{fixed}.$$

$$B,x\to \infty ,z=2\left(x-B\right)\text{fixed}.$$

$${\tilde{R}}_{\mathrm{B}}(x)=-2\pi \sqrt{x^2-B^2}$$

$$+\sum _{n=0}^{\infty }\sum _{m=0}^{\infty }\sum _{k=0}^{n+m+1}\frac{c_{n,m,k}(x/B{)}^{p_{\mathrm{B}}(k)}}{B^{m-n-1}{\left(x^2-B^2\right)}^{n+\frac{1}{2}}}{\ell }^k(x,B),$$

$$\ell (x,B)=\log \left(\frac{x-B}{x+B}\right).$$

$$R_{\mathrm{F}}(x)=-\log \left(\frac{x-B}{x+B}\right)$$

$$+\sum _{n=1}^{\infty }\sum _{m=0}^{\infty }\sum _{k=0}^{n+m}\frac{c_{n,m,k}(x/B{)}^{p_{\mathrm{F}}(k)}}{B^{m-n}{\left(x^2-B^2\right)}^n}{\ell }^k(x,B),$$

$${\int }_{-B}^B{f}_B(x)\mathrm{d}x=\pi B^2+B\sum _{m=0}^{\infty }\frac{c_{0,m,0}-2c_{0,m,1}}{B^m}.$$

$$R(z)={\int }_0^{\infty }\hat{R}(s){\mathrm{e}}^{-sz}\mathrm{d}s.$$

$$\sin (s)\left({\mathrm{e}}^{-\mathrm{i}s}{\hat{R}}_{\mathrm{B}}(s^{-})+{\mathrm{e}}^{\mathrm{i}s}{\hat{R}}_{\mathrm{B}}(s^{+})\right)$$

$$\cos (s)\left({\mathrm{e}}^{-\mathrm{i}s}{\hat{R}}_{\mathrm{F}}(s^{-})-{\mathrm{e}}^{\mathrm{i}s}{\hat{R}}_{\mathrm{F}}(s^{+})\right)$$

$$\hat{R}(s)=\Phi (s)\left(\frac{1}{s}+Q(s)\right),$$

$${\Phi }_{\mathrm{B}}(s)$$

$${\Phi }_{\mathrm{F}}(s)$$

$$Q_{\mathrm{B}}(s)$$

$$Q_{\mathrm{F}}(s)$$

$${\tilde{R}}_{\mathrm{B}}(z)=-2\pi \sqrt{Bz}+\frac{\sqrt{B}}{\sqrt{z}}\left(c_{0,0,0}+c_{0,0,1}\log \left(\frac{z}{4B}\right)\right)+\cdots$$

$$\frac{{\Phi }_{\mathrm{B}}(s)}{s}=\sqrt{\pi B}\left(\frac{1}{s\sqrt{s}}+\frac{1-{\gamma }_E+\log \left(\frac{\pi }{s}\right)}{\pi \sqrt{s}}+\cdots \right),$$

$${\tilde{R}}_{\mathrm{B}}(z)=\sqrt{\pi B}\left(-2\sqrt{\pi z}+\frac{1+\log (4\pi z)}{\sqrt{\pi z}}+\cdots \right).$$

Full text

I Introduction

The Lieb–Liniger (LL) and Gaudin–Yang (GY) models describe a gas of one-dimensional bosons (respectively, spin $1/2$ fermions) with a delta function interaction. As shown in the pioneering papers Lieb and Liniger (1963); Gaudin (1967); Yang (1967), these models are exactly solvable with the Bethe ansatz, and for this reason they are excellent laboratories to test our understanding of interacting quantum gases. In addition, they can be realized in the laboratory, leading to a fascinating interplay between integrable quantum models and experiment (see Jiang et al. (2015); Guan et al. (2013); Cazalilla et al. (2011); Bloch et al. (2008) for reviews and references).

Typically, exact solutions can be used to generate series in the coupling constant of the model, providing in this way a powerful test of perturbative methods. However, as already noticed in Lieb and Liniger (1963), the weakly coupled regime of the LL and GY models turns out to be very difficult to analyze by using the Bethe ansatz equations. Although the exact solution of these models has been known for fifty years, and in spite of many efforts Popov (1977); Iida and Wadati (2005); Tracy and Widom (2016a, b, 2018), only the first three terms of the perturbative series for the ground state energy density of these models have been derived. More recently, precise numerical calculations have made it possible to guess the analytic form of some additional coefficients Prolhac (2017); Lang (2018); Ristivojevic (2019).

In this paper we solve this longstanding problem and give a systematic procedure to derive the perturbative expansion for the ground state energy density of these models, building on previous results in Volin (2010, 2011). The procedure is recursive and can be easily implemented in a symbolic program. In this way we are able to compute analytically the perturbative series up to order $50$ in the coupling constant. We use our results to determine their large order behavior, which is an important source of non-perturbative information. We find that, in both cases, the series are factorially divergent. In the LL and attractive GY models, the series are non-alternating and therefore non-Borel summable. In the case of the attractive Gaudin–Yang model, whose ground state is a BCS-like state, we show that the associated singularity in the Borel plane is closely related to the superconducting energy gap.

II The Lieb–Liniger and Gaudin–Yang models

The Hamiltonian for the LL and GY models is given by (we set $\hbar=2m=1$)

[TABLE]

In the LL case we have a repulsive interaction with $c>0$, while in the GY case we will consider the attractive case $c<0$. We will focus on the ground state energy density of the models: $E=\lim_{L\rightarrow\infty}E_{0}(N,L)/L$, where $E_{0}(N,L)$ is the ground state energy for $N$ particles on the interval $[0,L]$ with periodic boundary conditions. We consider the thermodynamic limit in which $L,N\rightarrow\infty$ while the density $n=N/L$ remains fixed. The value of $E$ can be obtained from the solution of a linear integral equation for the density of Bethe roots, which is supported on an interval $[-B,B]$ (see Samaj and Bajnok (2013) for an excellent presentation of the Bethe ansatz solution of these models). In the LL model, the integral equation reads

[TABLE]

The dimensionless coupling constant $\gamma=c/n$ is given by

[TABLE]

while

[TABLE]

In the GY model, the integral equation reads

[TABLE]

The dimensionless coupling constant $\gamma=|c|/n$ is now given by

[TABLE]

while

[TABLE]

Note that the ground state energy density depends on the coupling constant through $B$, the endpoint of the support of the distribution. The strong coupling limit of these equations, $\gamma\rightarrow\infty$, can be studied systematically in a rather straightforward way (see e.g. Guan et al. (2013); Ristivojevic (2014)). However, the weak coupling limit $\gamma\rightarrow 0$ corresponds to $B\rightarrow\infty$, which is a singular limit of the integral equations. For this reason, extracting systematically the perturbative expansion from these integral equations has been an open problem since Lieb and Liniger introduced their equation in 1963. The known analytic results are Tracy and Widom (2016a, b, 2018)

[TABLE]

III Exact perturbative results

The method we will use to determine these series goes back to the work of Hutson on the related mathematical problem of the circular plate condenser Hutson (1963), and to the work of Popov on the LL equation Popov (1977). It was substantially improved by D. Volin in Volin (2010, 2011), where he applied it to integrable field theories in two dimensions (see also Kostov et al. (2008)). We will explain now this method in some detail, since it is likely to have other applications. It has three ingredients. First, one rewrites the integral equation as a difference equation for the resolvent of the density of Bethe roots. Second, one studies two different regimes for the equation: the bulk regime near the origin, and the edge regime near the endpoint of the support of the distribution. One then matches asymptotic expansions in the two regimes (this was already the main idea in Hutson (1963); Popov (1977)). The third ingredient involves finding an explicit solution for the difference equation in the edge regime.

Let us then first introduce the resolvent of the densities:

[TABLE]

(When no subscript is introduced, it is assumed that the equations are valid for both LL and GY models.) The resolvent is analytic in the complex $x$-plane, except for a discontinuity in the interval $[-B,B]$, given by

[TABLE]

where $R^{\pm}(x)=R(x\pm{\rm i}\epsilon)$. Note that $R(x)$ is an odd function on the real axis away from $[-B,B]$ (its discontinuity (10) is an even function of $x$). From its definition we deduce that

[TABLE]

In particular, we can compute the coupling constant $\gamma$ as a function of $B$ from the residue at infinity of the resolvent. By writing

[TABLE]

is now easy to see that the LL and Gaudin’s integral equation can be transformed into difference equations for the corresponding resolvents:

[TABLE]

where $\tilde{R}_{\rm{B}}(x)=R_{\rm{B}}(x)-2\pi x$. The subtraction of the linear term takes into account the known leading order behavior $f_{\rm{B}}(x)\approx 2{\sqrt{B^{2}-x^{2}}}$ Lieb and Liniger (1963), and leads to a simpler procedure.

Since weak coupling is large $B$, we study the resolvent in a systematic expansion in $1/B$. To do this, we consider the resolvent in two different regimes. In the bulk regime we take the limit

[TABLE]

In the edge regime we consider instead

[TABLE]

For the bulk regime, and inspired by Hutson (1963); Popov (1977); Iida and Wadati (2005); Volin (2010), we propose the following ansatz in the LL case,

[TABLE]

where $p_{\rm{B}}(k)=k$ mod $2$, and

[TABLE]

In the GY case, it is given by

[TABLE]

where $p_{\rm{F}}(k)=k-1$ mod $2$. Note that they are both odd functions on ${\mathbb{R}}\backslash[-B,B]$. By computing their discontinuity we find an expansion for the density of roots similar to the one considered in Hutson (1963); Popov (1977); Iida and Wadati (2005). The coefficients $c_{n,m,k}$ have all the information needed to compute the integrals (3), (4), (6), (7) as power series in $B$. For example, in the LL model, one has

[TABLE]

The $c_{n,m,k}$ can be partially determined by using the difference equations (13), but this is not enough to fix their value. We then look at the edge regime.

To study the edge regime, we write $R(z)$ as a Laplace transform,

[TABLE]

The difference equations (13) become discontinuity equations along the negative $s$-axis:

[TABLE]

where $s^{\pm}=s\pm{\rm i}\epsilon$. We have simply denoted by $\hat{R}_{\rm{B}}(s)$ the inverse Laplace transform of $\tilde{R}_{\rm{B}}(z)$. Order by order in a $1/B$ expansion, the functions $\hat{R}(s)$ have to be analytic everywhere except on the negative real axis, and they must have an expansion as $s\rightarrow\infty$ in integer, negative powers of $s$ (this can be justified as in Volin (2010, 2011)). The structure of poles and zeroes can be read from (21). The most general solution of (21) which satisfies these properties is the following:

[TABLE]

where

[TABLE]

and

[TABLE]

The coefficients $Q^{{\rm{B}},{\rm{F}}}_{n,m}$ are undetermined functions of $\log B$. Finally, we match the bulk and the edge regimes as in Hutson (1963); Popov (1977); Volin (2010): we first consider the ansatz (16) and (18) for $R_{{\rm{B}},{\rm{F}}}(x)$ in the bulk regime, we set $x=z/2+B$, and we expand around $z=\infty$ and $B=\infty$. We do the Laplace transform of the resulting expression, and we compare the result to the expansion of (22) around $s=0$. This turns out to fix all the coefficients $c_{n,m,k}$ and $Q_{n,m}$, therefore the ground state energy density of the LL and the GY models.

An example will illustrate how this works. By expanding (16), we find

[TABLE]

The expansion of $\Phi_{\rm{B}}(s)/s$ at $s=0$ reads,

[TABLE]

and the Laplace transform of (22) gives

[TABLE]

By comparing the coefficients we obtain,

[TABLE]

The recursive procedure to determine the coefficients $c_{n,m,k}$ can be fully automatized, as in Volin (2010, 2011), and one can push the calculation to any desired order, the only limitation being CPU time. We obtain in this way two series in $1/B$, one for $\gamma$ and another for $\langle x^{2}\rangle$, whose coefficients are polynomials in $\log B$. When we re-express $e_{{\rm{B}},{\rm{F}}}(\gamma)$ in terms of $\gamma$, logarithms disappear in the final result, and one obtains a power series in $\gamma^{1/2}$ for the LL case, and in $\gamma$ for the GY case. In the LL case, we find

[TABLE]

while in the GY case we find,

[TABLE]

We have calculated the exact form of the first $50$ coefficients for both the LL and the GY series. Our results are in accord with the numerical calculations and conjectures of Prolhac (2017); Lang (2018); Ristivojevic (2019) for the very first coefficients. As a bonus, we note that (19) is (up to an overall factor) the asymptotic expansion of the capacity of the circular plate condenser, a problem which goes back to the XIXth century Hutson (1963). The method explained in this paper gives a systematic and explicit procedure to determine this expansion.

IV Large order behavior

One of the motivations for this work is the study of the large order behavior of the perturbative series in quantum many-body systems. Since the pioneering work of Bender and Wu Bender and Wu (1973), we know that this behavior provides crucial information on the non-perturbative structure of the model. This connection between perturbative and non-perturbative sectors has now evolved into the theory of resurgence (see Mariño (2014); Aniceto et al. for reviews). Although the large order behavior of perturbation theory has been investigated in Fermi systems in e.g. Baker (1971); Rossi et al. (2018), we know comparatively less than in quantum mechanics or even quantum field theory, due in part to the scarcity of examples where the perturbative series can be computed to high order. Our results on the LL and the GY models provide enough data to determine with some precision the large order behavior of the series for the ground state energy density. Let us denote,

[TABLE]

One finds, for $k\gg 1$,

[TABLE]

The values $8\pi$, $\pi^{2}$ have been guessed from the numerics with a precision of six decimal digits, see Fig. 1 for an illustration in the LL case. We see that both series are factorially divergent and non-alternating, and there are singularities in the positive real axis of the Borel plane at $8\pi$ (for LL) and $\pi^{2}$ (for GY). Therefore, the series are not Borel summable. However, one can perform lateral Borel resummations by deforming the integration contour slightly above or below the positive real axis. By standard arguments (see e.g. Mariño (2015)), the resulting resummations have an ambiguous imaginary part which is exponentially small. It follows from (32) that these non-perturbative ambiguities are of the form

[TABLE]

The perturbative series for the attractive GY model is obtained by changing $\gamma\rightarrow-\gamma$ (note however that the exact energies are not related by this transformation). The resulting series is therefore alternating and very likely it is Borel summable. In addition, its Borel resummation seems to agree with the exact result (see Mariño and Reis ).

The result for the non-perturbative ambiguity in the GY model is particularly interesting. Instanton estimates indicate that, in fermionic theories, the growth of the perturbative series should be slower than factorial Parisi (1977); Baker and Pirner (1983), due to cancellations among diagrams. Our explicit results show that, at least in the GY model, this does not occur, and suggest that the divergence of the perturbative series can not be explained by semiclassical effects. In addition, one observes that the non-perturbative ambiguity in the GY model has the same exponential dependence than the square of the superconducting or spin gap of the model, which goes like $\Delta\sim{\rm e}^{-\pi^{2}/(2\gamma)}$ Krivnov and Ovchinnikov (1975); Fuchs et al. (2004). This is consistent with the fact that, in BCS theory, the ground state energy receives corrections of order $\Delta^{2}$. In Mariño and Reis we argue that this behavior might be universal, namely, that the large order behavior of the perturbative series in weakly interacting Fermi systems is closely related to the superconducting gap. We also argue that this behavior is due to the factorial growth of certain classes of diagrams, akin to the renormalon singularity in quantum field theory.

V Conclusions

In conclusion, in this paper we have solved the longstanding problem of determining the perturbative series of the ground state energy density in the LL and GY models, from the exact Bethe ansatz solution. The method we have developed, building on Volin (2010, 2011), is likely to have many other applications (see Mariño and Reis for some illustrations). The results on the large order behavior of the perturbative series open the possibility to study their relation to non-perturbative phenomena, and more generally to apply the tools of the theory of resurgence in many-body physics. In the case of the GY model, we have argued that the non-perturbative ambiguity is closely related to the superconducting energy gap, providing in this way a new perspective on the Cooper instability. It would be very interesting to unravel the precise source of the non-perturbative ambiguity in the LL model, and to test recent resummation techniques like the one presented in Rossi et al. (2018) in these exactly solvable models.

Acknowledgements: We would like to thank Sylvain Prolhac, Wilhelm Zwerger and specially Thierry Giamarchi and Félix Werner for useful discussions and comments on the manuscript. This work is is supported in part by the Fonds National Suisse, subsidies 200021-175539, and by the NCCR 51NF40-182902 “The Mathematics of Physics” (SwissMAP).