Resurgence for superconductors

TL;DR

This paper explores the connection between the energy gap in superconductors and the mathematical theory of resurgence, proposing that the perturbative series diverges factorially with a singularity linked to the energy gap.

Contribution

It conjectures a deep link between the superconducting energy gap and the divergence of perturbative series, supported by exact calculations in the Gaudin-Yang model.

Findings

Perturbative series for the ground state energy diverges factorially.

Leading Borel singularity is governed by the superconducting energy gap.

Identified a class of diagrams responsible for factorial growth.

Abstract

An important non-perturbative effect in quantum physics is the energy gap of superconductors, which is exponentially small in the coupling constant. A natural question is whether this effect can be incorporated in the theory of resurgence. In this paper we take some steps in this direction. We conjecture that the perturbative series for the ground state energy of a superconductor is factorially divergent, and that its leading Borel singularity is governed by the superconducting energy gap. We test this conjecture in detail in the attractive Gaudin-Yang model, an exactly solvable model in one dimension with a BCS-like ground state. In order to do this, we develop techniques to calculate the exact perturbative series of its ground state energy up to high order. We also argue that the Borel singularity is of the renormalon type, and we identify a class of diagrams leading to factorial growth. We give additional evidence for the conjecture in other models.

Equations

$${\Delta }^2\sim {\mathrm{e}}^{-A/\gamma }$$

$$V(x)=-2c\delta (x).$$

$$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),$$

$$N\to \infty ,L\to \infty ,n=\frac{N}{L}\text{fixed},$$

$$E=\underset{L\to \infty }{\lim }\frac{E_0(N,L)}{L}.$$

$$k_F=\frac{\pi n}{2}.$$

$$\gamma =\frac{c}{n}.$$

$$e(\gamma )=e_0+\sum _{n\ge 1}{c}_n{\gamma }^n,$$

$$e(\gamma )=\frac{E(\gamma )}{n^3}$$

$$e_0=\frac{{\pi }^2}{12}$$

$$c_1=-\frac{1}{2}.$$

$$c_2=-\frac{1}{12}.$$

$$c_3=-\frac{\zeta (3)}{{\pi }^4}.$$

$$\Pi (q)=-\mathrm{i}{\int }_{{\mathbb{R}}^{d+1}}\frac{{\mathrm{d}}^{d+1}k}{(2\pi {)}^{d+1}}G(k+q)G(k).$$

$$\frac{\mathrm{i}}{2\ell }{\int }_{{\mathbb{R}}^{d+1}}\frac{{\mathrm{d}}^{d+1}q}{(2\pi {)}^{d+1}}{\left(\kappa \Pi (q)\hat{V}(\mathbf{q})\right)}^{\ell },$$

$$-\frac{1}{\ell }{\int }_{{\mathbb{R}}^d}\frac{{\mathrm{d}}^d\mathbf{q}}{(2\pi {)}^d}{\int }_0^{\infty }\frac{\mathrm{d}\omega }{2\pi }{\left(\kappa \Pi (\mathbf{q},\mathrm{i}\omega )\hat{V}(\mathbf{q})\right)}^{\ell }.$$

$$E_{\mathrm{RPA}}={\int }_{{\mathbb{R}}^d}\frac{{\mathrm{d}}^d\mathbf{q}}{(2\pi {)}^d}{\int }_0^{\infty }\frac{\mathrm{d}\omega }{2\pi }\left\{\log \left(1-\kappa \Pi (\mathbf{q},\mathrm{i}\omega )\hat{V}(\mathbf{q})\right)+\kappa \Pi (\mathbf{q},\mathrm{i}\omega )\hat{V}(\mathbf{q})\right\}.$$

$$\Pi (q,\mathrm{i}\omega )=-\frac{m}{2\pi k_{F}y}\log \left(\frac{(y/2+1{)}^2+{\nu }^2}{(y/2-1{)}^2+{\nu }^2}\right),$$

$$y=\frac{q}{k_F},\nu =\frac{m\omega }{q{k}_F}.$$

$$\lambda ={\left(\frac{\kappa }{2}\right)}^2\gamma .$$

$$e(\lambda ;\kappa )=\frac{1}{4}\frac{E/\kappa }{(n/\kappa {)}^3}.$$

$$e(\lambda ;\kappa )=\sum _{n=0}^{\infty }{e}_n(\lambda ){\kappa }^{-n}.$$

$$e_0(\lambda )=\frac{{\pi }^2}{12}-\lambda .$$

$$e_1(\lambda )=\lambda -\frac{\pi }{4}\sum _{\ell =2}^{\infty }\frac{1}{\ell }{\left(\frac{4\lambda }{{\pi }^2}\right)}^{\ell }{\mathcal{I}}_{\ell },$$

$${\mathcal{I}}_{\ell }={\int }_0^{\infty }\mathrm{d}yy{\int }_0^{\infty }\mathrm{d}\nu {\left(\frac{1}{2y}\log \left(\frac{(y/2+1{)}^2+{\nu }^2}{(y/2-1{)}^2+{\nu }^2}\right)\right)}^{\ell }.$$

$${\Pi }_{\mathrm{pp}}(P)=\mathrm{i}{\int }_{{\mathbb{R}}^{d+1}}\frac{{\mathrm{d}}^{d+1}q}{(2\pi {)}^{d+1}}{G}_{\mathrm{p}}(P/2-q)G_{\mathrm{p}}(P/2+q),$$

$${\Pi }_{\mathrm{pp}}(\mathbf{P},\mathbf{k})=m{\int }_{{\mathbb{R}}^d}\frac{{\mathrm{d}}^d\mathbf{q}}{(2\pi {)}^d}\frac{\theta (|\mathbf{P}/2+\mathbf{q}|-k_{\text{F}})\theta (|\mathbf{P}/2-\mathbf{q}|-k_{\text{F}})}{{\mathbf{k}}^2-{\mathbf{q}}^2+\mathrm{i}ϵ}.$$

$$-2c{\int }_{{\mathbb{R}}^d}\frac{{\mathrm{d}}^d\mathbf{P}}{(2\pi {)}^d}{\int }_{{\mathbb{R}}^d}\frac{{\mathrm{d}}^d\mathbf{k}}{(2\pi {)}^d}\frac{\theta (k_{\text{F}}-|\mathbf{P}/2+\mathbf{k}|)\theta (k_{\text{F}}-|\mathbf{P}/2-\mathbf{k}|)}{1+2c{\Pi }_{\mathrm{pp}}(\mathbf{P},\mathbf{k})}.$$

$$s=\frac{|\mathbf{P}|}{2k_F},t=\frac{|\mathbf{k}|}{k_F}.$$

$$\mathrm{Re}{\Pi }_{\mathrm{pp}}(P,k)=\frac{m}{2\pi k_F}{F}_{\mathrm{pp}}(s,t),$$

Full text

1 Introduction

The energy gap of a superconductor is a robust property of Fermi systems with a weak attraction. It has a non-analytic, exponentially small dependence on the coupling constant bcs , and it should be therefore regarded as a non-perturbative effect. It is well known that many non-perturbative effects in physics, from quantum mechanics to string theory, can be detected by looking at the large-order behavior of the perturbative series. In fact, they provide a detailed characterization of the perturbative series at large orders, as first shown in the pioneering work of Bender and Wu bw ; bw2 . This relationship between perturbative series and non-perturbative effects has evolved into the theory of resurgent asymptotics, which can be considered as a general framework to understand non-perturbative phenomena in physics and mathematics (see mmlargen ; abs for reviews and references).

It is then natural to ask whether the theory of resurgence has something to say about the energy gap of a superconductor, providing in this way a characterization of the behavior of the perturbative series for weakly interacting fermion systems. In this work we take some steps in this direction. We argue that the perturbative series for the ground state energy of a superconductor is factorially divergent and non-Borel summable, and that the corresponding non-perturbative ambiguity is closely related the energy gap. More precisely, let

[TABLE]

be the square of the superconducting energy gap, where $\gamma$ is the coupling constant. Then, we conjecture that the first singularity in the Borel plane of the coupling constant is given by $A$. The fact that the non-perturbative ambiguity involves $\Delta^{2}$ can be justified heuristically by noting that in BCS theory the non-analytic correction to the ground state energy is proportional to $\Delta^{2}$.

Some aspects of our proposal have been hinted before in the literature, although the behavior of perturbative series in many-body theory has not been studied extensively. Previous works pointed out the possibility that the series diverges factorially baker-review ; jw , although instanton estimates in Fermi systems have suggested a milder growth of the coefficients parisi-fermi ; baker-pirner , due to cancellations among diagrams. The Cooper instability has been sometimes linked to the lack of analyticity of the perturbative series, but no precise connection seems to have been made between these two phenomena. One possible reason for the scarcity of studies of large order behavior in many-body theory is probably the lack of examples where the perturbative series can be calculated to sufficiently large order. In Quantum Mechanics, calculable models have been crucial in the development of the theory of resurgence. For example, the groundbreaking results of Bender and Wu in bw ; bw2 were made possible to a large extent by the data provided by the perturbative series for the ground state energy of the quartic oscillator.

In order to make precise statements about the perturbative series in many-body fermionic systems, it is useful to have an example akin to the anharmonic oscillator in quantum mechanics, where the ground state energy can be explicitly computed to high order in perturbation theory. A natural candidate is the Gaudin–Yang model in one dimension gaudin ; yang . The ground state of the Gaudin–Yang model is a BCS-like state, and there is a superconducting gap which is non-perturbative in the coupling constant ko (see giamarchi for this and other results on one-dimensional models). At the same time, the energy of the ground state can be computed exactly with the Bethe ansatz, through a simple integral equation.

However, extracting the perturbative series from this integral equation turns out to be highly non-trivial. The same difficulty was observed in the closely related Lieb–Liniger model ll . In spite of many efforts iw-gy ; tw1 ; tw2 ; guan-review , and after fifty years, only the first three terms of the series have been obtained analytically (further coefficients have been obtained numerically in prolhac ). In this paper, building upon previous work hutson ; popov ; volin ; volin-thesis , we develop a method to calculate analytically the perturbative series from the Bethe ansatz equations, and we obtain explicit results for its coefficients up to order $50$ in the coupling constant (in principle, one can go further by using enough CPU time). This allows us to study the large order behavior of the series and to determine (numerically) the leading non-perturbative ambiguity. We find that the first singularity in the Borel plane corresponds indeed to the energy gap, in agreement with the conjecture above.

In general, there are two different sources for the factorial behavior of a perturbative series (see mmbook for a review): instanton singularities and renormalon singularities. Instanton singularities are due to non-trivial saddles in the path integral, and they encode the factorial growth in the number of diagrams. Renormalon singularities are typically associated to sequences of diagrams in which each diagram diverges factorially with the loop order, after integration over momenta (see e.g. beneke for a review and references). We argue that the large order behavior found in the Gaudin–Yang model should be attributed to renormalon-type singularities. To support this, we show explicitly that both ring and ladder diagrams lead to factorially divergent sequences.

We give further evidence for our conjecture by examining an even simpler example, the one-dimensional Hubbard model at half-filling lw . Here the series for the ground state energy is known in closed form, and we verify our conjecture analytically. We also examine the much more complicated case of the dilute Fermi gas in three dimensions. Although a systematic calculation of the ground state energy perturbative series at large order is currently out of reach, we find some encouraging indications that it might have the appropriate renormalon-type singularities.

This paper is organized as follows. In section 2 we present various approaches to the study of the ground state energy density of the Gaudin–Yang model: the perturbative approach, the BCS approach, and finally the exact Bethe ansatz solution. Most of section 2 is a review of known facts, but some of the material included here, like the explicit perturbative results for the ring and ladder diagrams, seem to be new. In section 3 we present in detail the method to extract the perturbative series from the Bethe ansatz equation. We apply it to the Gaudin–Yang model in the attractive and repulsive regimes, and to its generalization for $\kappa$ spin components. In section 4 we study the large order behavior of the series, we make the connection to the superconducting energy gap, and we argue that the Borel singularity is of the renormalon type. In section 5 we consider the one-dimensional Hubbard model and the dilute Fermi gas in three dimensions in the light of our conjecture. Finally, in section 6, we conclude and and list many open problems. There are also two Appendices. In the first one we evaluate the NNNLO diagram (of order $\gamma^{3}$) in the Gaudin–Yang model. In the second one we give some details of the matching procedure explained in section 3.

2 The Gaudin–Yang model

2.1 The perturbative approach

The Gaudin–Yang model describes a non-relativistic, one-dimensional gas of $N$ spin $1/2$ fermions with $\delta$-function interactions. The potential is given by

[TABLE]

The Hamiltonian of the system is

[TABLE]

i.e. we set $\hbar=2m=1$. We will mostly consider the attractive case. In our conventions, this corresponds to a positive coupling constant: $c>0$. Let $E_{0}(N,L)$ be the ground state energy for the $N$-particle system on a circle of length $L$. In the thermodynamic limit of fixed density $n$,

[TABLE]

we can define the ground state energy density as

[TABLE]

The Fermi momentum $k_{F}$ is related to the density by

[TABLE]

Since $c$ has the dimension of an inverse length, the dimensionless coupling constant is

[TABLE]

Note that weak coupling, $\gamma\ll 1$, corresponds here to high density, in contrast to what happens in the three dimensional Fermi gas with a delta function interaction. The ground state energy density can be written as a formal power series in $\gamma$ of the form

[TABLE]

where

[TABLE]

and

[TABLE]

corresponds to the free Fermi gas. The coefficients $c_{n}$, $n\geq 1$, can be calculated in terms of Feynman diagrams. The coefficient $c_{1}$ is given by the sum of Hartree and Fock terms, and its calculation is an elementary exercise:

[TABLE]

The coefficient $c_{2}$ is obtained by the sum of the two-bubble diagram in Fig. 1 and the diagram obtained by exchange from it. It has been calculated in magyar , and its value is

[TABLE]

The next coefficient was computed in magyar2 , and its value is:

[TABLE]

In Appendix A we rederive this result for an arbitrary number of spin components. Clearly, pushing this calculation to higher orders becomes increasingly hard, and as far as we know, no other results for higher coefficients have been obtained in the literature.

It is however possible to consider partial resummations of the perturbative series involving special types of diagrams. The most popular ones are the ring diagrams appearing in the RPA, and the ladder diagrams. In particular, one can obtain interesting weak-strong coupling interpolations for Fermi gases in three dimensions by resumming these diagrams steele ; hf ; schafer ; kaiser1 ; kaiser2 ; kaiser2d .

Let us then start by considering the resummation of ring diagrams, i.e. the random phase approximation (RPA). The building block of the RPA approximation is the so-called polarization loop or Lindhard function, which is given in $d$ dimensions by

[TABLE]

Here, $G(k)$ is the free fermion propagator at $T=0$. In terms of the Lindhard function, the contribution of a ring diagram with $\ell\geq 2$ bubbles to the ground state energy density is

[TABLE]

where $\kappa$ is the number of spin components and $\widehat{V}({\bf q})$ is the Fourier transform of the potential. By rotating the integration contour and taking into account that $\Pi({\bf q},{\rm i}\omega)$ is an even function of $\omega$, we can write (2.14) as

[TABLE]

It is possible to sum these terms in closed form to find the usual RPA contribution to the correlation energy,

[TABLE]

We note that in one dimension, the Lindhard function can be written as electron-liquid

[TABLE]

where the variables $\nu$, $y$ are defined as

[TABLE]

One important property of ring diagrams is that they give the first non-trivial contribution to the ground state energy density in the limit of a large number $\kappa$ of spin components (see e.g. hf ). In this limit, which is akin to a large $N$ limit, one has to scale the coupling constant in the potential and the density $n$ in such a way that every interaction gives a power of $\kappa^{-1}$ and that the Fermi momentum is fixed.

Let us now analyze the ring diagrams and the large $\kappa$ limit in the case of the Gaudin–Yang model. For a delta function potential (2.1) $\widehat{V}({\bf q})=-2c$. The Fermi momentum is given by $k_{F}=\pi n/\kappa$, therefore the appropriate scaled coupling constant is

[TABLE]

It plays the rôle of a ’t Hooft parameter and is kept fixed in the large $\kappa$ limit. We also define the rescaled ground energy density as

[TABLE]

It has the following large $\kappa$ expansion,

[TABLE]

Note that, when $\kappa=2$, $\lambda$ and $e(\lambda;\kappa)$ become the coupling constant $\gamma$ and the energy density $e(\gamma)$, respectively. The leading order function $e_{0}(\lambda)$ to (2.21) is given by the free gas result plus the Hartree term,

[TABLE]

The subleading function $e_{1}(\lambda)$ is the sum of the contributions of ring diagrams,

[TABLE]

where the first term corresponds to the Fock diagram, and

[TABLE]

We can represent this result schematically as in Fig. 2.

Another useful reorganization of the perturbative series is the particle-particle ladder expansion. In this approach, one considers diagrams with the shape of a ladder with $n$ steps, where each of the steps is an interaction. The ladder is then closed into a circle, as shown in Fig. 3, in order to have a vacuum diagram. We will follow the approach in steele ; schafer ; he-huang and consider ladders in which two of the parallel lines are “hole” lines (i.e. their momenta are inside the Fermi sphere), while the rest of the lines are “particle” lines (i.e. with momenta outside the Fermi sphere). We will refer to the resulting diagrams as pp-ladders (or particle-particle ladders), although a more appropriate name would be $n$pp-$1$hh bubbles he-huang . We will also include the diagrams obtained from the ladders by exchange of the hole lines.

Let us summarize the calculation of pp ladder diagrams, following the presentation in he-huang . The building block of the pp ladder approximation is the particle-particle bubble in $d$ dimensions,

[TABLE]

where $G_{\rm p}(k)$ is the particle part of the propagator (in which $|{\bf k}|>k_{F}$). In the calculation of the ground state energy, one needs the on-shell version of $\Pi_{\rm pp}(P)$, which is given by

[TABLE]

Here $\theta(x)$ is the Heaviside function. For a delta function interaction (2.1), the contribution of pp ladders to the ground state energy is a geometric series which can be resummed as

[TABLE]

Note that only the real part of $\Pi_{\rm pp}({{\bf P}},{\bf k})$ contributes to this integral. We introduce now the dimensionless variables

[TABLE]

Specializing now to one dimension, we find

[TABLE]

where

[TABLE]

We conclude that the pp ladder contribution to $e(\gamma)$, which we will denote by $e_{\rm pp}(\gamma)$, is given by

[TABLE]

This can be re-expanded in power series of $\gamma$ to obtain,

[TABLE]

where the very first coefficients can be computed explicitly in closed form,

[TABLE]

The first two coefficients correspond to the diagrams shown in Fig. 1, i.e. to the Hartree–Fock contributions and to the contributions at order $\gamma^{2}$. Note that the perturbative series up to order $\gamma^{2}$ is correctly reproduced by the ladder resummation.

2.2 The BCS approach

Perturbation theory has of course many limitations, in particular if one is interested in the strong coupling regime where $\gamma\gg 1$. Resummations of particular types of diagrams can be sometimes extrapolated to this regime schafer ; kaiser1 , but it is important to have other approaches to the problem. One possibility is to use BCS theory. Indeed, in the attractive case, the Hamiltonian (2.2) leads to Cooper pairing, and the ground state of the Gaudin–Yang model can be regarded as a superconductor. BCS theory leads to approximate expressions for the energy gap and for the ground state energy which interpolate between weak and strong coupling. The application of BCS theory to the Gaudin–Yang model was developed in montse ; quick , and we will complement the results in those papers by presenting explicit expressions for all quantities in terms of elliptic functions, as it was done in strinati in the three-dimensional case.

In the BCS approach, we assume that the ground-state wavefunction is of the BCS type and we use the variational or mean-field method (a detailed exposition can be found in section 37 of fw ). To write down the resulting equations, we introduce

[TABLE]

where $\mu$ is the chemical potential and $n$ is the density. The variational method leads to the famous gap equation for the energy gap $\Delta_{\rm BCS}$, which in this case reads

[TABLE]

In addition, we have to relate $\mu$ to the density, and this gives the equation

[TABLE]

The BCS expression for the ground state energy per particle is

[TABLE]

We note that this energy incorporates the Hartree–Fock contribution, which is sometimes included only in the fluctuations around the mean-field solution.

It turns out that all the integrals appearing above can be calculated in closed form in terms of elliptic integrals of the first and the second kind. To write explicit expressions, we introduce the variable

[TABLE]

as well as the square elliptic modulus

[TABLE]

The gap equation can now be written as

[TABLE]

where $K(m)$ is the elliptic integral of the first kind. The equation relating the density to $\mu$ is

[TABLE]

where $E(m)$ is the elliptic integral of the second kind. The energy per particle reads

[TABLE]

These expressions are valid for any value of the coupling $c$, but it interesting to analyze them in the weak coupling limit, where the gap is exponentially small. It follows from (2.39) that in this regime $m\approx 1$. Expanding the functions above around this limit, we find an explicit expression for the gap,

[TABLE]

where $E_{F}=k_{F}^{2}$ is the Fermi energy. Let us note that this is not the correct expression for the energy gap as $\gamma\rightarrow 0$, but rather the approximation to it provided by BCS theory. However, it captures already the correct exponentially small dependence on the coupling constant, as we will see. The ground state energy is given by

[TABLE]

where $E_{0}$ is the energy of the free Fermi gas. In terms of $e_{\rm BCS}=E_{\rm BCS}/n^{3}$, we find

[TABLE]

This expression includes the Hartree–Fock energy, together with exponentially small corrections in $\gamma$. We note that the non-analytic correction is proportional to the square of the gap, $\Delta_{\rm BCS}^{2}$.

2.3 The Bethe ansatz approach

One of the most important aspects of the Gaudin–Yang model is that its ground state energy can be computed exactly by using the Bethe ansatz. This exact solution was obtained in gaudin ; yang for fermions of spin $1/2$, and extended in taka to fermions with $\kappa$ spin components. An excellent presentation of the spin $1/2$ case can be found in sb-book . We will now summarize some ingredients of the solution which will be useful in the following.

We will first consider the attractive case. The basic ingredient is the density of Bethe roots $\phi(k)$ describing the ground state (which has zero polarization). This density has a compact support $[-K,K]$ and satisfies the so-called Gaudin’s integral equation

[TABLE]

The density of the gas is given by

[TABLE]

while the ground state energy density is

[TABLE]

We will now introduce a convenient normalization due to tw1 ; tw2 . We introduce the function

[TABLE]

so that Gaudin’s integral equation reads

[TABLE]

Here, $B=K/c$. The coupling (2.6) is given by

[TABLE]

while the normalized ground state energy density (2.8) is

[TABLE]

Let us note that, in the Bethe ansatz solution, the basic parameter is $B$, the endpoint of the interval where $f(x)$ is supported. The coupling constant $\gamma$ is a non-trivial function of $B$, and $e(\gamma)$ depends on $\gamma$ through $B$. In addition, the weak coupling limit $\gamma\rightarrow 0$ correspond to the limit of large $B$.

Unfortunately, the Gaudin integral equation can not be solved analytically. Therefore, although the Bethe ansatz solution makes it possible to calculate the ground state energy density exactly, one usually has to resort to numeric methods or to an asymptotic analysis in particular limits. It turns out that the study of the Gaudin equation in the weakly coupled regime is particularly difficult, since the integral equation becomes singular when $c\rightarrow 0$. The strong coupling regime is much easier to analyze, and one finds that as $\gamma\rightarrow\infty$,

[TABLE]

The Bethe ansatz solution presented above describes the ground state with zero polarization, but it can be extended to arbitrary polarization. This makes it possible to calculate exactly the energy gap (also known as spin gap). The energy gap is defined as (minus) the energy difference between the ground state, in which there are equal numbers of spin up and spin down particles, and the first excited state, in which one breaks a Cooper pair and obtains in this way a triplet state with total spin $S=1$. It is easy to see from the Bethe ansatz solution that the spin gap can be calculated by solving an integral equation (see e. g. ko ; zhou-exact ), and one finds, at leading order in a $\gamma$ expansion ko ; frz ,

[TABLE]

Note that the exponent is identical to the one obtained in the BCS approximation to the gap (2.43).

The case of a repulsive interaction with $c<0$ can be also solved with the Bethe ansatz, but one obtains a different integral equation for the ground state energy density (see sb-book for a useful presentation). This equation can be put in the form

[TABLE]

The kernel appearing in this equation is

[TABLE]

More explicitly, one has

[TABLE]

where $\psi(z)$ is the digamma function. Note that this kernel is no longer algebraic, and it is similar to the ones appearing in the Bethe equation for integrable field theories in two dimensions (like the $O(n)$ sigma model or the Gross–Neveu model). The coupling constant in the repulsive case is defined as

[TABLE]

and it is given by

[TABLE]

The ground state energy density is

[TABLE]

As in the attractive case, it is difficult to use the integral equation above to study the weakly coupled gas. However, one can easily study the strongly coupled limit, and one finds (see e.g. guan-review and references therein)

[TABLE]

Finally, the Bethe ansatz solution for the ground state energy density can be extended to fermions with $\kappa$ spin components and an attractive interaction taka (see e.g. guan-multi for additional information and references). The density of Bethe roots satisfies now the integral equation

[TABLE]

where the kernel is given by

[TABLE]

The density $n$ and energy $E$ are now given by

[TABLE]

where

[TABLE]

The ’t Hooft-like parameter $\lambda$ introduced in (2.19) is given by

[TABLE]

and the rescaled energy density defined in (2.20) can be obtained from

[TABLE]

One can use the Bethe ansatz exact solution to assess the accuracy of well-known approximations to the calculation of the ground state energy density. For example, in montse ; quick , the exact solution was compared to the BCS result in the attractive case. The BCS approximation provides a powerful weak/strong coupling interpolation procedure, and one can see that it is very successful, as illustrated in Fig. 4. As noted in montse ; quick , it matches exactly both the weak coupling and the strong coupling limits.

The resummation of diagrams also provides a mechanism to interpolate between weak and strong coupling. In schafer ; kaiser1 , the resummation of pp ladder diagrams in the three-dimensional Fermi gas with a contact interaction was studied, and it was used to find an estimate of the Bertsch parameter of the unitary gas. A similar calculation can be done in the repulsive Gaudin–Yang model. The ladder resummation (2.31) provides an estimate for the ground state energy, $e_{\rm pp}(-\gamma)$, which can be compared to the exact answer $e_{r}(\gamma)$ obtained from the Bethe ansatz solution. This is illustrated in Fig. 4. In this respect, it is interesting to define the analogue of the Bertsch parameter $\xi$ for the repulsive Gaudin–Yang model, as

[TABLE]

In view of (2.61), the Bethe ansatz gives the exact answer

[TABLE]

while the pp ladder approximation gives

[TABLE]

On the other hand, this approximation is not suitable for the attractive regime, and it fails to capture even qualitatively the strong coupling behavior (2.53).

3 Exact results for the perturbative series

As we mentioned in the Introduction, the analysis of the Gaudin integral equation at weak coupling turns out to be highly non-trivial, since in this limit the integral equation becomes singular. A similar problem occurs in the Lieb–Liniger model, as it was already noted in ll . As a result, it is very difficult to determine analytically the perturbative series for the ground state energy density of these models. In spite of numerous efforts iw-gy ; tw1 ; tw2 , the series has only been derived up to second order in the coupling constant (even at the physical level of rigor). More recently, the combination of numerical analysis and inspired guesswork has led to conjectures for the exact form of some of the additional coefficients prolhac ; lang .

In this section, we will solve this longstanding problem and present a method that determines systematically the perturbative series for the ground state energy density of the Gaudin–Yang model from a direct analysis of the integral equation111This method can be also used to obtain the perturbative series for the Lieb–Liniger model, which is presented in a separate paper mr-ll .. The method goes back to the work of Popov on the Lieb–Liniger equation popov (and even earlier, to the work by Hutson on the circular plate condenser hutson ), but it was substantially improved by D. Volin in a remarkable tour de force volin ; volin-thesis . Volin’s work focused on Bethe ansatz equations appearing in the context of integrable field theories in two dimensions (in particular the $O(n)$ sigma model). Gaudin’s integral equation is an even simpler example where this method can be applied.

3.1 Resolvent in the bulk and edge regimes

The method we will use has three fundamental ingredients. The first one is to rewrite the integral equation as a difference equation for the resolvent of the density of Bethe roots. The second one is to consider two different regimes for the equation: the bulk regime near the origin, and the edge regime near the endpoints of the support. One then matches asymptotic expansions in the two regimes (this was already the main idea in hutson ; popov ). The third ingredient involves finding an explicit solution for the difference equation in the edge regime, after a Laplace transform.

We will now apply the method to Gaudin’s integral equation. As we will see in the next section, it is also possible to appeal to the conventional Wiener–Hopf method to obtain the wished-for solution for the resolvent in the edge regime without even writing the difference equation. This makes it possible to analyze a more general family of integral equations.

Let us first write the kernel of Gaudin’s integral equation (2.50) as

[TABLE]

where

[TABLE]

and

[TABLE]

is the shift operator which acts as $x\rightarrow x+{\rm i}$. Following volin ; volin-thesis (see also ksv ), we introduce the resolvent of the density of Bethe roots as

[TABLE]

This function is analytic in the complex $x$-plane except for a discontinuity in the interval $[-B,B]$, given by

[TABLE]

Let us note that $R(x)$ is an odd function on the real axis away from $[-B,B]$, while its discontinuity is an even function of $x$. From its definition we deduce that

[TABLE]

In particular, by using (2.51) we can compute the coupling constant $\gamma$ as a function of $B$ from the residue at infinity of the resolvent. It is now easy to see that Gaudin’s integral equation leads to the following difference equation for $R(x)$:

[TABLE]

The weak coupling limit of the Gaudin–Yang model corresponds to the limit of large $B$. Therefore, a natural strategy is to study the resolvent in a systematic expansion in $1/B$. To do this, we study the resolvent in two different regimes. The first one is the so-called bulk regime, in which we take the limit

[TABLE]

in such a way that

[TABLE]

is fixed. This is therefore appropriate to study $f(x)$ near $x=0$. The second regime is the so-called edge regime, in which we also have (3.78) but we keep fixed the variable

[TABLE]

This is therefore appropriate to study $f(x)$ near the edge of the distribution $x=B$.

For the bulk regime, and inspired by hutson ; popov ; iw-gy ; volin , we propose the following ansatz:

[TABLE]

where $p(k)$ is [math] or $1$ and

[TABLE]

Further justification for this ansatz will be clear when we study the edge regime. Note that $R(x)$ is indeed an odd function on ${\mathbb{R}}\backslash[-B,B]$. Its discontinuity at the interval $[-B,B]$ leads to the following ansatz for $f(x)$,

[TABLE]

where

[TABLE]

and $x\in[-B,B]$. The coefficients $c_{n,m,k}$ have all the information needed to compute the ground state energy. First of all, one has that

[TABLE]

We also have,

[TABLE]

The coefficients $c_{n,m,k}$ can be partially determined by using the difference equation (3.77), 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)=R(x(z))$ as a Laplace transform,

[TABLE]

As we will see in a moment, $R(s)$ decreases as $1/s$ at infinity, and the Bromwich inversion formula reads

[TABLE]

The inverse Laplace transform of (3.77) is obtained by replacing

[TABLE]

and it reads

[TABLE]

To find the correct solution to (3.90) we demand the following properties for $\hat{R}(s)$ at each order in a $1/B$ expansion:

$\hat{R}(s)$ is analytic everywhere except on the negative real axis. 2. 2.

$\hat{R}(s)$ has simple poles at $s/\pi=-n-1/2$, where $n$ is a non-negative integer. 3. 3.

$\hat{R}(s)$ has an expansion as $s\rightarrow\infty$ in integer, negative powers of $s$.

These analyticity properties can be justified as in volin ; volin-thesis . Note that the structure of zeroes in the resolvent follows from (3.90). The third property follows from the following consideration. The function $f(x)$ is regular at $x=B$, and therefore it has an expansion

[TABLE]

To reproduce this structure through the discontinuity formula we must have

[TABLE]

We now use the following (analytically continued) Laplace transform,

[TABLE]

In this way, the inverse Laplace transform of (3.92) gives

[TABLE]

which is precisely the third property above. The connection between (3.91) and (3.94) is very useful since it makes it possible to test the solution for $\hat{R}(s)$ against a numerical study of the distribution $f(x)$ at the edge.

The most general solution of (3.90) which satisfies the analyticity properties stated above is the following:

[TABLE]

where

[TABLE]

and $Q(s)$ is given by

[TABLE]

The dependence of $Q(s)$ on $B$ is a consequence of the structure of (3.81). We note that the coefficients $Q_{n,m}$ are undetermined functions of $\log B$. The overall coefficient in $\Phi(s)$ is fixed by studying (3.90) near $s=0$. In the next section we will give a different argument, based on the Wiener–Hopf method, to obtain the function $\Phi(s)$.

It turns out that the information contained in (3.96) is enough to compute the coefficients $c_{n,m,k}$ and $Q_{n,m}$, therefore the ground state energy density of the Gaudin–Yang model. This is done by “matching” the bulk and the edge regimes, as already proposed in hutson ; popov . In Volin’s reformulation of the method, we first consider the ansatz (3.81) for $R(x)$ in the bulk regime, by setting $x=z/2+B$ and expanding around $B=\infty$, $z/B=0$. We do the Laplace transform of the resulting expression, and we compare the result to the expansion of (3.95) around $s=0$. This turns out to fix all the coefficients. In Appendix B we work out explicitly the very first terms to illustrate the procedure. Note that the terms obtained in the expansion of (3.95) provide a further check of the bulk ansatz (3.81). The general expressions (3.85) and (3.86) give, after using the values for the coefficients obtained in (B.4),

[TABLE]

They are in agreement with the results in tw1 ; tw2 . We also obtain, for the distribution $f(x)$ in the bulk regime,

[TABLE]

where $u=x/B$. The leading order term in this expression reproduces the result of iw-gy .

As in volin ; volin-thesis , the recursive procedure to determine the coefficients $c_{n,m,k}$ can be fully automatized in a symbolic program, 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 $\rho$ and another for $\langle x^{2}\rangle$, whose coefficients are polynomials in $\log B$. This leads to a series of the same type for $\gamma$, which can be inverted to obtain $1/B$ as a power series in $\gamma$ (involving as well $\log\,\gamma$). One can then re-express the r.h.s. of (2.52) as a power series in $\gamma$ in order to calculate the perturbative series for the ground state energy density. It turns out that the dependence on $\log\,\gamma$ fully disappears in the final result (it was erroneously claimed in ko that this is not the case). This is similar to what has been observed in integrable field theories in fnw1 ; bbbkp ; volin ; volin-thesis .

In this way one finds, up to order $\gamma^{10}$,

[TABLE]

Note that, up to $\gamma^{3}$, this agrees with the perturbative results obtained above. The coefficients in (3.100) are in accord with the numerical calculations of S. Prolhac prolhac . We have pushed the calculation up to order $50$ in $\gamma$, which makes it possible to determine the large order behavior of the series with some precision.

3.2 Relation to the Wiener–Hopf method

One popular approach to study the integral equations appearing in the Bethe ansatz is to use the Wiener–Hopf method (see e.g. hmn ; hn ; fnw1 for examples in the context of relativistic field theories). In the cases we are interested in here, the integral equation is of the form

[TABLE]

where $\alpha$ is a real constant. In the Wiener–Hopf method, one extends this equation to the whole real axis in order to use Fourier transform techniques. This method seems however to be less efficient than the one presented here. For example, in the study of the Gaudin–Yang model in tw1 ; tw2 with the Wiener–Hopf method, only the terms up to order $\gamma^{2}$ are obtained, in spite of some effort. A key aspect of our approach is the explicit solution (3.96), since this provides the initial data that allow to fix all the constants. This function is a solution to the discontinuity equation satisfied by $\hat{R}(s)$. However, as noted already in volin , one can see in many examples that $\Phi(s)$ is closely related to the Wiener–Hopf decomposition of the kernel ${\cal K}(x)$. This makes it possible to obtain $\Phi(s)$ without the need to determine and solve the equations satisfied by the resolvent. We will now give a heuristic argument to establish this relation. The argument below is only valid when $G_{\pm}(0)$ is finite and provides information about the strict large $B$ limit, but it can be extended to more general situations.

Let us write the integral equation (3.101) in the edge regime, in terms of the variable (3.80), as

[TABLE]

Here, $f(z)=f(x(z))$ and $K(z)={\cal K}(z/2)/2$, in order to absorb the multiplicative factor from the change of variables. By extending the above equation to all real $z$, we obtain

[TABLE]

where $\xi(z)$ is a priori an unknown function. Let us define the following Fourier transforms

[TABLE]

The subscripts $\pm$ denote that something is analytic in the upper/lower half complex plane (including the real axis but possibly excluding the origin). We also introduce

[TABLE]

When $K(\omega)$ is an even function of $\omega$, we have,

[TABLE]

This decomposition can almost always be done provided $1-K(\omega)$ is well defined on the real axis, though some care might be necessary at $\omega=0$. The integral equation becomes,

[TABLE]

By using that

[TABLE]

we can rewrite (3.107) as

[TABLE]

From (3.109) it follows that $C(\omega)$ must be analytic in both the upper and the lower half complex planes, reducing it to an entire function. However, in their respective half planes (including the real axis), $F_{-}(\infty)=0$, $X_{\pm}(\infty)=0$ and $G_{\pm}(\infty)$ is a constant (the first two results follow from their definitions and (3.103), while the latter can be checked explicitly in examples). Both sides are thus bound at infinity, and by Liouville’s theorem, we conclude that $C(\omega)=0$. Thus, in the strict limit $B\rightarrow\infty$, we have

[TABLE]

We can now relate $F_{-}(\omega)$ to $\hat{R}(s)$ introduced in (3.87) (again in the strict large $B$ limit). We consider $\hat{R}({\rm i}s)$ for ${\rm Im}(s)<0$ and deform the Bromwich contour ${\cal B}$ in (3.88) around the negative real axis, without crossing it. We obtain in this way minus the contour ${\cal C}$ shown in Fig. 5. Therefore,

[TABLE]

We conclude that

[TABLE]

which is precisely of the form of the first term in the r.h.s. of (3.95). In other words,

[TABLE]

The results (3.112) and (3.113) show that we can skip the steps of finding and solving an equation for the resolvent. We can simply obtain the crucial leading term of $\hat{R}(s)$ from the Wiener–Hopf decomposition of the kernel, which is available in many situations. Note however that (3.110) gives the solution in the strict large $B$ limit. There are $1/B$ corrections to this result which lead to the term $Q(s)$ in (3.95).

As an example, let us reconsider Gaudin’s integral equation. In the notation of (3.101) we have

[TABLE]

and one has (see e.g. tw1 )

[TABLE]

It is easy to check that (3.113) reproduces the result (3.96) for $\Phi(s)/s$.

3.3 The repulsive case

The Gaudin–Yang model with a repulsive interaction is a good example where one can obtain the wished-for result for the resolvent by using (3.112): the kernel has a very simple expression after Fourier transformation, and one can also find easily its Wiener–Hopf decomposition. In the notations of the previous section, one has

[TABLE]

and

[TABLE]

From these results and (3.113) we obtain the function $\Phi(s)$. The ansatz for the resolvent in the bulk is the same as in (3.81), with the only difference that the leading logarithmic term has an additional factor of $4$:

[TABLE]

The ansatz for the edge is (3.95), albeit with the new $\Phi(s)$. Following the procedure of section 3.1, we find that the ground state energy density is given by the series (3.100) after we change $\gamma\rightarrow-\gamma$, as expected. It is also possible to write down a difference equation for the resolvent of this model, by using the representation of the Gamma function given in volin ; volin-thesis .

3.4 The Gaudin–Yang model with $\kappa$ spin components

The Gaudin–Yang model with $\kappa$ spin components is particularly interesting, since we can consider the limit of large $\kappa$ and make contact with the RPA solution at subleading order in $1/\kappa$.

To analyze the linear integral equation (2.62), we use the ingredients provided by the Wiener–Hopf method. We have, in the notation of (3.102),

[TABLE]

and

[TABLE]

We can now use (3.113) to obtain $\Phi(s)$. The ansatz for the bulk and for the edge are the same as for the Gaudin–Yang model, with the only difference that in (3.95) we have to use the new $\Phi(s)$. By using the method in section 3.1, we can now proceed to calculate the ground state energy, as a function of $\lambda$ and $\kappa$ (we recall that $\lambda$ was defined in (2.19)). It is convenient to present the results in terms of the large $\kappa$ expansion in (2.21). One finds, for the first two series $e_{n}(\lambda)$,

[TABLE]

The result for $e_{1}(\lambda)$ gives explicit predictions for the contributions of the first ring diagrams. In particular, we obtain a Bethe ansatz prediction for the values of the integrals (2.24). One finds, for example,

[TABLE]

The first two values can be verified analytically, since ${\cal I}_{2}$ is proportional to the two-bubble diagram evaluated in $c_{2}^{\rm pp}$ in (2.33), while the value of ${\cal I}_{3}$ follows from the calculation in Appendix A. The remaining equalities can be verified numerically (although the direct numerical evaluation of ${\cal I}_{\ell}$ becomes more and more difficult as $\ell$ is increased.) A detailed analysis of ring diagrams for a delta-function interaction in two dimensions has been made in kaiser2 . Interestingly, their evaluation is simpler than in the one-dimensional case, and it does not involve odd values of the zeta function. In the next section we will also examine the behavior of the ${\cal I}_{\ell}$ as $\ell$ becomes large.

4 Large order behavior, energy gap, and renormalons

4.1 Large order behavior of the perturbative series

The results obtained in the previous section make it possible to study the large order behavior of the perturbative series with numerical methods. We will first consider the series for $e(\gamma)$ in the original case of spin $1/2$, for which we know the exact form of the first $50$ terms.

It is very easy to see that the coefficients in (2.7) grow factorially, i.e. $c_{k}\sim k!$. Following standard practice in the asymptotic analysis of perturbative series (see e.g. mmbook ), we assume the following ansatz:

[TABLE]

The parameters $A$, $b$, and $\mu_{i}$, $i\geq 0$, can be extracted systematically from the sequence $c_{k}$, by a numerical analysis of appropriate auxiliary sequences (see e.g. msw ). For example, the parameter $A$ can be extracted from the sequence

[TABLE]

whose asymptotics is

[TABLE]

The precision in the determination of these parameters depends of course on the number of coefficients of the series which are available. In order to accelerate the convergence of the auxiliary sequences, we will use Richardson transforms bender-book ; msw . Our analysis gives

[TABLE]

with a precision of six decimal digits, once Richardson transforms are used. This is illustrated in the left graphic of Fig. 6. In addition, we have

[TABLE]

In other words, we find the asymptotic behavior

[TABLE]

Our numerical analysis also indicates that

[TABLE]

although this is somewhat less precise. In the graphic on the right in Fig. 6 we show the sequence

[TABLE]

together with its fifth Richardson transform. According to (4.6), this sequence should asymptote $1$, as we observe in the graphic on the right in Fig. 6.

These empirical results show that the perturbative series of the Gaudin–Yang model is factorially divergent. In the attractive case, the series is non-alternating, therefore not Borel summable. If we consider the Borel transform of $e(\gamma)$,

[TABLE]

the large order behavior (4.6) indicates that this function has its first Borel singularity at

[TABLE]

This is precisely twice the strength of the exponentially small term in the energy gap (2.54), and it confirms our conjecture in the case of the Gaudin–Yang model. As we noted in the Introduction, the factor of two can be justified heuristically by the fact that the non-analytic correction to the energy in the BCS approximation is proportional to $\Delta_{\rm BCS}^{2}$, as we saw in (2.44).

Let us now recall some basic facts of the theory of Borel–Écalle resummation (see e.g. mmbook ). When the Borel transform $\widehat{\varphi}(\zeta)$ of a perturbative series $\varphi(z)$ has singularities on the positive real axis, one can define lateral Borel resummations as

[TABLE]

where ${\cal C}_{\pm}$ are integration paths slightly above (respectively, below) the positive real axis. If the coefficients of the perturbative series are real, the two lateral Borel resummations have the same real part but imaginary parts with opposite sign. On the other hand, the large order behavior (4.1) is associated to a non-perturbative ambiguity or non-perturbative correction of the form:

[TABLE]

This ambiguity gives the leading contribution to $s_{+}(\varphi)(z)-s_{-}(\varphi)(z)$.

In the case of the Gauding–Yang model, it follows from (4.4), (4.5) and (4.12) that the non-perturbative ambiguity has the following form,

[TABLE]

at leading order in $\gamma$ and ${\rm e}^{-\pi^{2}/\gamma}$ as $\gamma\rightarrow 0$. This implies that the lateral Borel resummation of the perturbative series for $e(\gamma)$ has an imaginary piece given by

[TABLE]

Let us consider a numerical example to illustrate this result. It is convenient to parametrize $\gamma$ by $B$, following (2.51). Let us denote by $e_{\pm}^{\rm BP}(\gamma)$ the lateral Borel–Padé resummations of the perturbative series, in which we use a Padé approximant to extend analytically $\widehat{e}(\gamma)$ along the positive real axis (as usual in this technique, we obtain an approximate value for $e_{\pm}^{\rm BP}(\gamma)$ by keeping the digits which are stable as we increase the number of coefficients in the series). For $B=3$ one finds $\gamma(3)=0.4438...$, and

[TABLE]

On the other hand, the exact ground state energy, obtained by numerical solution of the Bethe ansatz, is

[TABLE]

Note that the real part of (4.15) agrees with the exact value $e(\gamma(3))$ with a precision of eighteen digits. In addition, the r.h.s. of (4.14) for this value of $\gamma$ is $0.98...\cdot 10^{-10}$, which compares well with the imaginary part of (4.15).

In order to recover the exact ground state energy, the perturbative series (3.100) should be extended to a full trans-series mmlargen ; abs , including exponentially small corrections in the coupling constant. The trans-series is expected to be of the form,

[TABLE]

The first term in the r.h.s. is the conventional perturbative series, while the second term is the infinite series of exponentially small corrections. $C_{\ell}$ and $b_{\ell}$ are parameters of the trans-series, and we note that the $C_{\ell}$ are generically complex. From the form of the ambiguity (4.13) it follows that $b_{1}=1$. We conjecture that the first coefficient has the value

[TABLE]

and plays the rôle of a Stokes parameter. As usual, the choice of sign is correlated with the choice of lateral resummation. To justify this, we note that with the value (4.18), the resummation of the first exponential correction to the trans-series cancels the imaginary part of the lateral Borel resummations of the perturbative series. In addition, we note that the real part of the Borel resummation of the perturbative series agrees with the exact result with a much higher precision than $\exp(-\pi^{2}/\gamma)$. Therefore, the first correction to the real part is likely to be of order $\exp(-2\pi^{2}/\gamma)$, and this implies that the real part of $C_{1}$ vanishes. We note that the structure of the trans-series appearing here is very similar to the Painlevé II example studied in mmnp .

In the repulsive case, the coefficients become alternating: $c_{k}\rightarrow(-1)^{k}c_{k}$. The leading singularity is now on the negative real axis, and the series is likely to be Borel summable. Let us denote by $e_{r}^{\rm BP}(\gamma)$ the Borel–Padé resummation of the perturbative series. Up to the numerical accuracy with which we can evaluate this resummation, $e_{r}^{\rm BP}(\gamma)$ is in perfect agreement with the exact $e_{r}(\gamma)$ calculated from the Bethe ansatz. Let us give an example. We parametrize again $\gamma$ by $B$, as dictated by (2.59). Then, one has, for $B=1$, $\gamma(1)=2.249..$, and

[TABLE]

The Borel summability of this series in the repulsive case is in agreement with general consideration on the fermion many-body perturbative expansion baker-review ; baker-pirner .

4.2 Rings, ladders and renormalons

Non-perturbative ambiguities in asymptotic series are the signal of non-perturbative effects that have to be taken into account. It is then important to identify the source of the ambiguity found in the previous section. In general, in quantum theory there are two sources for the factorial growth of the coefficients in the perturbative series (see mmlargen for a review). The first one is the growth in the number of Feynman diagrams as we increase the order. This is believed to be encoded in the instantons of the model, i.e. in non trivial saddle points of the Euclidean path integral. Another source of factorial growth are renormalons. Renormalons are special types of diagrams whose contribution grows factorially as one increases the order. This growth is not due to an increase in the number of diagrams, but to the integration over momenta. Typical examples of renormalons are bubble diagrams in QED or QCD beneke .

In theories with fermions, instanton methods suggest that the growth of the perturbative series is milder than the standard factorial growth, due to cancellations among diagrams parisi-fermi ; baker-pirner . For example, in some models in dimension three, instanton dominance in the large order behavior gives a growth of order $(n!)^{1/5}$, instead of $n!$ baker-pirner 222In rossi , the perturbative series for a Fermi gas in 3D with contact interaction was resummed in terms of an effective coupling, and a growth of $(n!)^{1/5}$ was predicted with instanton techniques. The setting is different from our since the effective coupling includes in particular the resummation of ladder diagrams. We would like to thank Félix Werner for clarifications on this issue.. This instanton-induced growth is expected to become even slower as we decrease the dimension. As we have seen, the ground state energy density in the Gaudin–Yang model has a fully factorial growth, so it seems unlikely that this is due to an instanton effect.

We turn then to the other possible source of factorial growth: renormalons. Are there sequences of diagrams in the Gaudin–Yang model which grow factorially after integration over momenta? If so, can they lead to the behavior we have found empirically? Before answering these questions, let us point out that the analysis of special classes of diagrams by itself does not allow to find the precise location of the Borel singularities. In the study of renormalons in quantum field theory, this is typically complemented by considerations based on the renormalization group and the OPE (see beneke ). However, the existence of these diagrams indicates that instanton estimates of the large order growth are likely to be inaccurate.

In the study of renormalons in gauge theories with fermions, it has been proved useful to organize the diagrams by considering the limit of a large number of flavours. We can follow the same procedure here, since as we have explained there is a similar limit in many-fermion theories: the limit of a large number of spin components. This limit is dominated by the ring diagrams appearing in (2.23). The first question to ask is then: what is the behavior of the integrals ${\cal I}_{\ell}$ as $\ell$ grows large? A little bit of numerical experimentation indicates that indeed, they grow factorially. This can be seen as follows. Let us consider the formal power series,

[TABLE]

which is essentially the Borel transform of the series appearing in (2.23). The factorial behavior

[TABLE]

is equivalent to the statement that the singularity of $\varphi(\zeta)$ which is closest to the origin occurs at $\zeta=C$. We can locate this singularity as follows. By using the explicit expression (2.24), we resum the series defining $\varphi(\zeta)$ as

[TABLE]

where

[TABLE]

The most important contribution to the above integral comes from the neighbourhood of the logarithmic singularity at $\nu=0$, $y=2$. If we introduce the variable

[TABLE]

the singularity occurs now at $z=\nu=0$. We can approximate the integral defining $\varphi(\zeta)$ by

[TABLE]

where $\Lambda_{z}$, $\Lambda_{\nu}$ define a region around the singularity. It is easy to see that this integral is finite for $0<\zeta<4$, but diverges as $\zeta\rightarrow 4$. The first singularity occurs at

[TABLE]

We conclude that

[TABLE]

and the ring diagram with $\ell$ bubbles diverges then as

[TABLE]

This corresponds to a non-perturbative ambiguity at large $\kappa$ of the form ${\rm e}^{-\pi^{2}/\lambda}$. When $\kappa=2$, the exponential dependence is precisely the one found in (4.13). As a side remark, let us point out that, as a consequence of this large order behavior, the series defining $e_{1}(\lambda)$ in (2.23) is factorially divergent. This illustrates the fact that, in the presence of renormalons, the large $N$ expansion fails to produce convergent series order by order in $1/N$, as noted in e.g. fkw .

The above analysis shows that there are renormalon-type singularities in the Gaudin–Yang model: ring diagrams grow factorially, due to integration over the momenta. Note that the source of the growth of the integral is the logarithmic divergence of the Lindhard function at $\omega=0$, $q=2k_{F}$, which is due to particle-hole excitations near the Fermi surface (this is in contrast to conventional renormalons in relativistic field theory, where the factorial growth comes from large or small momenta). The strength of the factorial growth is precisely what is needed to reproduce the observed large order behavior.

It should be noted however that the factorial growth of the ring diagrams is specific to one dimension. In dimensions two and three, the Lindhard function does not have a logarithmic singularity, and no factorial growth of the ring diagrams is observed for a contact interaction. Is there are another set of diagrams which might lead to a factorial growth of the renormalon type? A natural candidate is the set of ladder diagrams shown in Fig. 3. The factorial growth of this class of diagrams has been noted long ago baker-review , and they play an important rôle in the study of the Cooper instability (see e.g. as ).

Let us then look at the pp ladder diagrams considered in section 2.1. Like before, let us consider the series

[TABLE]

which is the Borel transform of (2.32). It can be written in closed form as

[TABLE]

The most important contribution to the integral comes from $s=0$, $t=1$. This corresponds to a pair of holes with opposite momentum near the Fermi surface. To extract this contribution, we can approximate the integral by

[TABLE]

which has its first singularity at

[TABLE]

We deduce that the coefficients $c_{n}^{\rm pp}$ behave as

[TABLE]

This can be also confirmed by a direct numerical evaluation of the coefficients. In this case, the singularity associated to this sequence of diagrams occurs at $\zeta=2\pi^{2}$, which is twice the value we have obtained in the previous section by studying the whole sequence. Therefore, the pp ladder diagrams underestimate the growth of the full sequence. This is not surprising, since we are not even taking into account the full contribution of ladder diagrams, but only those with two hole lines. In two kaiser2d and three dimensions kaiser1 , it is possible to include as well ladders with more than one pair of hole lines. This might change the location of the Borel singularity, and as we will see in section 5, in three dimensions it does. It would be very interesting to extend the calculations of kaiser1 ; kaiser2d to one dimension, and check whether the inclusion of additional ladder diagrams gives the correct singularity in the Borel plane.

In conclusion, there are factorially divergent subdiagrams in the Gaudin–Yang model leading to the correct singularity (in the case of ring diagrams) or to a closely related one (in the case of pp ladder diagrams). In our view, this supports the argument that the non-perturbative ambiguity we have found in the Gaudin–Yang model is of the renormalon type.

5 More examples

The connection between the non-perturbative ambiguity and the superconducting gap in weakly interacting Fermi systems is very natural, and as we have seen in the last section, there is clear numerical evidence that it holds in the Gaudin–Yang model. It would be interesting to investigate whether it holds in other models. In this section we look at two additional examples. Our first example, the one-dimensional Hubbard model at half-filling, is even simpler than the Gaudin–Yang model, and we can establish the conjecture analytically. The second example, the dilute Fermi gas in three dimensions, is much more difficult to analyze, but we can still give some partial evidence for our conjecture.

5.1 The Hubbard model

The Hubbard model in one dimension describes $N$ interacting fermions of spin $1/2$ in a lattice with $L$ sites. The Hamiltonian is

[TABLE]

where $c_{l\sigma}^{\dagger}$, $c_{l\sigma},$ $l=1,\cdots,L$, $\sigma=\uparrow,\,\downarrow$ are creation/annihilation operators for the fermions. This model can be regarded as a discretized version of the Gaudin–Yang model, and we will focus again on the attractive case with $c<0$.

The ground state of the system corresponds to the so-called half-filled band, in which $N=L$, and with a vanishing polarization. The model can be solved with Bethe ansatz techniques lw , and the solution leads to an exact expression for the energy of the ground state. In the thermodynamic limit, the ground state energy density can be written down in closed form in terms of Bessel functions:

[TABLE]

(We follow here the conventions in sb-book .) This has a asymptotic series expansion in powers of $c$ which was obtained in mis-ov . It reads,

[TABLE]

where

[TABLE]

The very first terms in this series can be tested against conventional perturbation theory metzner , as we have done for the Gaudin–Yang model in previous sections. Using that $\zeta(x)\rightarrow 1$ as $x\rightarrow\infty$, one can determine analytically the large order behavior of the $a_{k}$,

[TABLE]

The series diverges doubly-factorially and it is non-alternating, therefore non-Borel summable. The large order behavior (5.5) of the $a_{k}$ leads to a non-perturbative ambiguity of order

[TABLE]

On the other hand, the Hubbard model is known to have an energy gap which, in the half-filled case, is of the form (see e.g. marsiglio ; giamarchi )

[TABLE]

As before, the ambiguity (5.6) has the exponential dependence of the square of the gap. It would be very interesting to understand the origin of the Borel singularity by looking at the diagrammatic structure of the theory, as we have done for the Gaudin–Yang model.

5.2 The dilute Fermi gas in three dimensions

The two models discussed so far are the Gaudin–Yang model and the one-dimensional Hubbard model. They are both integrable and one can obtain explicit expressions for the coefficients of the perturbative series to large order (in the case of the Hubbard model with a half-filled band, to all orders). It would be interesting to explore other models where our conjecture can be tested, in particular models in higher dimensions, although an important limitation in this study is the lack of concrete results on the coefficients of the perturbative series.

As a modest first step in this direction, we will make a preliminary analysis of a dilute Fermi gas in three dimensions with a contact interaction. The Lagrangian of the system is written as

[TABLE]

where $M$ is the mass of the fermions. One can include additional interactions in order to match the scattering data. In this matching, the constant $C_{0}$ is related to the scattering length $a$ by $C=4\pi a/M$. The natural dimensionless coupling is

[TABLE]

where $k_{F}$ is the Fermi momentum. We follow the convention in kaiser1 , namely, an attractive interaction corresponds to $g>0$. One can expand the ground state energy per particle in a power series in $g$. Note that, in contrast to the Gaudin–Yang model, the weak coupling limit of small $g$ is here a limit of low density. Let us define

[TABLE]

where $N$ is the number of particles. Then, after including the diagrams shown in Fig. 1, one finds the well-known result (see e.g. bishop ; hf-effective ),

[TABLE]

where we have set for simplicity $\kappa=2$. This series is the three-dimensional counterpart to (3.100), but it is only known up to order $g^{4}$ fourth-order . An additional subtlety as compared to (3.100) is that $e^{\rm 3D}(g)$ also includes logarithmic corrections (i.e. terms involving $\log(g)$), so one should enlarge the formal framework we have used in the Gaudin–Yang model to take into account this additional expansion parameter. However, we can have a first diagnose on the large order behavior of (5.11) by looking at renormalon effects similar to those unveiled in the Gaudin–Yang model. In three dimensions, ring diagrams grow only exponentially, so we should look at ladder diagrams. As in section 2.1, we first consider pp ladder diagrams, which contain one hole-hole line and $n$ particle-particle lines. They can be resummed in a compact form (see e.g. steele ; schafer ; kaiser1 ; he-huang ), and the result reads

[TABLE]

Here, the function $F_{\rm pp}(s,t)$ is the three-dimensional analogue of (2.30), and it reads

[TABLE]

The function $I(s,t)$ was introduced in kaiser2 and it is a kinematic function to perform the integration over hole momenta. Its explicit expression is

[TABLE]

We can now re-expand $e^{\rm 3D}_{\rm pp}(g)$ in powers of $g$, to obtain

[TABLE]

where

[TABLE]

It is also possible to include and resum ladders with more than one pair of hole lines kaiser1 ; kaiser2d (see also he ). One finds,

[TABLE]

where the function $R(s,t)$ is given by

[TABLE]

After re-expanding, we find

[TABLE]

We can now ask what is the large order behavior of the coefficients $c_{n}^{\rm 3D,\,pp}$, $c_{n}^{\rm 3D,\,lad}$. This question can be addressed either with the Borel resummation techniques we used for the ring and ladder diagrams in the Gaudin–Yang model, or by inspecting numerically the behavior of the coefficients at large $n$. One finds,

[TABLE]

This can be clearly seen in the plots in Fig. 7, although the actual large order behavior sets in relatively late. In the pp case, for example, the series seems to grow exponentially (and to be alternating) until $n\sim 25$.

The first conclusion to extract from (5.20) is that ladder diagrams indeed grow factorially and lead to non-Borel summable series. The second conclusion is that, as we enlarge the type of diagrams included in the resummation, the location of the first Borel singularity changes as well: the pp ladders lead to a singularity at $4\pi$, while the improved resummation in kaiser1 moves it to $2\pi$. The associated non-perturbative ambiguities are

[TABLE]

We can now compare this to the result for the superconducting gap, which is of the form (see e.g. pitaevski-rev )

[TABLE]

According to our conjecture, the perturbative series should lead to a Borel singularity at $\pi$, which is half the singularity obtained with the improved ladder resummation of kaiser1 . This apparent mismatch shouldn’t be discouraging, since our conjecture concerns the full series, while the results above concern only a subclass of diagrams. Indeed, the location of the singularity changes by a factor of two as we pass from pp ladders to ladders including at least two holes lines. It is conceivable that including more diagrams would confirm our conjecture. In addition, one should remember that the perturbative series for the three-dimensional gas includes logarithms, so probably more care should be exercised. We find it encouraging that well-known resummations of the ground state energy of the dilute Fermi gas lead to Borel singularities which are closely related to the superconducting energy gap. In any case, a more detailed investigation is needed in order to determine the large order behavior of the ground state energy for this system.

6 Conclusions

The main goal of this paper has been to start the study of perturbative series in many-body physics with the tools of the theory of resurgence. History teaches us that a useful starting point is a non-trivial, yet quasi-solvable model where one can obtain extensive data on its perturbative structure. This is what Bender and Wu did in quantum mechanics: they first obtained data for the perturbative series of the ground state energy of the quartic oscillator bw , and then they discovered its precise connection to tunneling bw2 .

In this paper we have tried to adapt this program to many-body fermionic systems. Our toy model has been the Gaudin–Yang model, which can be solved exactly with the Bethe ansatz. However, extracting the perturbative series for the ground state energy from the exact solution is already non-trivial. Fortunately, we were able to rely on the powerful techniques of volin ; volin-thesis , and we found a systematic method to obtain the perturbative information. The series turns out to be factorially divergent and non-Borel summable, and we have shown that its large order behavior is controlled by the superconducting energy gap of the model.

In retrospect, the validity of our conjecture is very natural in view of what happens in the Gross–Neveu model gross-neveu , which can be regarded as a relativistic cousin of the attractive Gaudin–Yang model. At large $N$, which corresponds to the BCS mean-field theory, the mass gap of the Gross–Neveu model can be calculated exactly and is of the form fnw2

[TABLE]

where $\beta_{0}$ is the first coefficient of the beta function. On the other hand, this non-perturbative effect controls the large order behavior of the conventional perturbative series in the model, and it is a renormalon effect. Our conjecture is a transposition of these well-known facts to the non-relativistic case.

The energy gap is a robust non-perturbative effect in weakly coupled Fermi systems with an attractive interaction, and it sets the non-perturbative scale of the problem. Therefore it is natural to conjecture that this connection is general, and that the resurgent structure of superconducting systems is governed by the superconducting gap. We have also stated a quantitatively precise connection, in (1.1), which holds in the Gaudin–Yang model and in the closely related example of the one-dimensional Hubbard model. Our conjecture should be regarded as an exploration tool in a research program, so it might need qualification or modification in more general cases. Our preliminary analysis of the dilute Fermi gas in three dimensions indicates that the large order behavior of special types of diagrams goes into the right direction, but more work is clearly necessary. Recent progress on high-order calculations (as in e.g. montecarlo ) might shed light on these issues.

Our results seem to provide an interesting picture of the Cooper instability. It is often stated that, due to this instability, the perturbative expansion around the conventional vacuum is incorrect, and one should expand around the superconducting vacuum with a non-zero gap. Our findings make this statement more precise: the conventional perturbative series is a non-Borel summable series, which lacks non-perturbative information. In quantum theory, non-Borel summable series are often a signal that one is not expanding around the right ground state, or that one is not taken into account additional vacua. At the same time, the theory of resurgence teaches us that the shortcomings of conventional perturbation theory can be cured in many cases by extending the perturbative series to a trans-series. Therefore, once the full trans-series for the superconductor is obtained, it should be in principle possible to extract its physical ground state energy.

We have also argued that the non-perturbative ambiguity associated to the lack of Borel summability does not have a semiclassical interpretation, and is due to a renormalon-like effect. To this end, we have identified sequences of diagrams which display the typical factorial growth of renormalon physics, coming from integration over momenta. In the case of ladder diagrams, which are known to encode information about the Cooper instability, the kinematic region that leads to the factorial growth corresponds indeed to Cooper pairs near the Fermi surface.

The connection between the Cooper instability and renormalon singularities leads to a very intriguing interplay, and should be understood further333Some parallelisms between superconductivity and the large order behavior of perturbative series in QCD were pointed out in peris . Renormalon-like diagrams were found in the $\epsilon$ expansion of the unitary Fermi gas in nishida-son .. It would be interesting to see whether standard tools in renormalon physics, like the use of RG equations and the OPE, can be adapted to many-body physics and shed light on the location of the Borel singularities. This interplay might also lead to a more tractable arena to understand the phenomenon of renormalons in quantum field theory. It has been argued in recent years that renormalons have semiclassical realizations once the theory is suitable compactified argyres-unsal ; a-unsal-long ; dunne-unsal-cpn ; cherman-dorigoni-dunne-unsal ; cdu ; misumi-1 ; du-on ; misumi-2 . It is not clear whether this scenario holds in purely fermionic theories. For example, the Gross–Neveu model, which is historically the first model displaying renormalons gross-neveu , has not been given so far a semiclassical picture. It would be then interesting to see whether the Borel singularity we have found in the Gaudin–Yang model can be interpreted in semiclassical terms along the lines proposed in those works.

In addition to these more general issues, it would be fruitful to apply the techniques used in this paper to other models. One could for example study the Gaudin–Yang model with arbitrary polarization, for which some terms in the series are known iw-pol ; guan-review . More generally, one can apply the method developed in this paper to other integrable theories in two dimensions mr-ift . The techniques we have developed can be also applied to the Lieb–Liniger model mr-ll , and one obtains in this way explicit results for the perturbative series of its ground state energy density. The series obtained in that case is also factorially divergent and non-Borel summable. Given the importance of the Lieb–Liniger model, both theoretically and experimentally, it would be very interesting to have a precise physical understanding of the corresponding non-perturbative ambiguity. It is likely that this is due to an instanton in the non-relativistic version of the $O(2)$, two-dimensional $\lambda\phi^{4}$ theory describing the Lieb–Liniger model.

Another important question is to obtain the precise form of the trans-series for the Gaudin–Yang model. This trans-series is encoded in the exact Bethe ansatz solution, but it is not obvious how to extract it. Perhaps one should first address this question in the half-filled, one-dimensional Hubbard model, where the energy can be written down in closed form, as we have seen in (5.2). A closely related question would be to find the all-orders expansion of the spin gap of the Gaudin–Yang model from the Bethe ansatz. The leading order term quoted in (2.54) was determined in ko up to an overall constant, which was fixed in frz by using results from the Hubbard model444We would like to thank W. Zwerger for detailed explanations of this point.. It would be very interesting to see whether the techniques used in this paper to study the ground state energy can be also used to calculate the subleading corrections to the energy gap.

More generally, it would be important to explore with the tools of resurgence the vast landscape of models of interacting fermions, and to study the connection conjectured in this paper between Borel singularities and the non-perturbative energy gap. This might lead to a deeper understanding of non-perturbative effects in these theories.

Acknowledgements

We would like to thank Lianyi He, Norbert Kaiser, Sylvain Prolhac, Ricardo Schiappa, and Wilhelm Zwerger for useful discussions and correspondence. We are specially indebted to Thierry Giamarchi for his patient explanations of one-dimensional models, and to Félix Werner for a detailed reading of 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).

Appendix A Perturbative evaluation of the NNNLO contribution

The contribution to the ground state energy of the Gaudin–Yang model at order $\gamma^{3}$ was obtained in magyar2 by using standard perturbative techniques. Here we rederive this result and we consider an arbitrary number of spin components $\kappa$. We can adapt the results of bishop , who studied the perturbative expansion of the dilute Fermi gas in three dimensions. At order $\gamma^{3}$ there are three types of Goldstone diagrams that contribute: a purely pp ladder diagram (together with its exchange), a ladder diagram with 2 particle lines and 4 hole lines, and finally diagrams with 3 particle lines and 3 hole lines (which include the three-bubble ring diagram). The total contribution is of the form

[TABLE]

where $I_{j}$, $j=1,2,3$ correspond to the types of diagrams that we have just listed. $I_{1}$ comes from the pp ladder, which according to (2.31) is given by

[TABLE]

$I_{2}$ is given by

[TABLE]

where $n(k)=\theta(1-|k|)$. After some manipulations, we obtain

[TABLE]

Finally, we have

[TABLE]

It is convenient so split this integral in four parts, as follows,

[TABLE]

This can be evaluated as

[TABLE]

Since

[TABLE]

for $\kappa=2$ we obtain the result (2.12). This confirms the Bethe ansatz calculation of the term of order $\gamma^{3}$ in (3.100). In addition, since we also have the explicit spin factors, we can also verify that the contribution to $e(\lambda;\kappa)$ of order $\lambda^{3}$ is

[TABLE]

in agreement with the results for $e_{1,2}(\lambda)$ in (3.121).

Appendix B Matching the bulk to the edge

Let us illustrate the method developed in section 3.1 by working out in detail how the matching of bulk and edge is done. By expanding (3.81) around the edge regime, we find

[TABLE]

On the other hand, the expansion of (3.95) at $s=0$ reads,

[TABLE]

where $\gamma_{E}$ is the Euler–Mascheroni constant. Its Laplace transform gives

[TABLE]

By comparing the coefficients we obtain,

[TABLE]