Testing the Bethe ansatz with large N renormalons

The ground-state energy of integrable asymptotically free theories can be conjecturally computed using the Bethe ansatz once the theory has been coupled to an external potential through a conserved charge. This leads to a precise prediction for the perturbative expansion of the energy. We provide a non-trivial test of this prediction in the non-linear sigma model and its supersymmetric extension, by calculating analytically the associated Feynman diagrams at next-to-leading order in the 1/N expansion, and at all loops. By investigating the large order behavior of the diagrams, we locate the position of the renormalons of the theory and we obtain an analytic expression for the large N trans-series associated to each. As a spin-off of our calculation, we provide a direct derivation of the beta function of these theories, at next-to-leading order in the 1/N expansion.


Introduction
In quantum field theory, the conventional perturbative series explores just the fluctuations around the trivial vacuum. It is believed that a full picture should involve additional sectors, corresponding to non-trivial vacua, and leading to exponentially suppressed contributions to physical observables. The origin of these sectors might be additional semi-classical saddles, also called instantons, but in many cases, like in Yang-Mills theory, the most important non-perturbative effects are associated to elusive objects called renormalons [1].
The contribution of instanton sectors can be obtained in principle by expanding the path integral around these non-trivial saddle-points. In the case of renormalons, we do not have such a semi-classical picture. It is known, however, that renormalons manifest themselves in the large order behavior of perturbation theory [2,3]. Therefore, one way to obtain precise quantitative information about renormalon sectors and their exponentially small corrections is to use conventional perturbation theory at large number of loops. This connection, in which a non-perturbative correction "resurges" in the perturbative series, is part of the general theory of resurgence, which gives a systematic framework to understand non-perturbative sectors (see [4][5][6][7][8] for introductions). In spite of some impressive lattice calculations [9] (see also [10]), unveiling the resurgent structure of realistic quantum field theories remains a difficult problem. One could however try to address these issues in more tractable quantum field theories. Among these, integrable field theories in two dimensions play an important rôle. On the one hand, they display many of the physical phenomena of interest, like asymptotic freedom and the presence of renormalon sectors. On the other hand, they are much more tractable analytically. In particular, exact expressions for their S matrices have been conjectured in [11,12].
It was noted by Polyakov and Wiegmann in [13,14] that, in integrable field theories, one can calculate exactly the dependence of the ground-state energy on a chemical potential coupled to a global conserved charge. This is done by combining the exact S-matrix with the Bethe ansatz, and the answer can be elegantly expressed in terms of a set of integral equations. One beautiful application of this observation was the determination of the exact mass gap of these theories, by comparing the Bethe ansatz to conventional perturbation theory (see [15][16][17][18][19] for various case studies, and [20] for a review). The ground-state energy as a function of the chemical potential, which we will call in short the free energy of the theory, is then an ideal observable to understand quantitatively the different sectors, perturbative and non-perturbative, of the theory.
However, even extracting the perturbative expansion of the free energy from the Bethe ansatz answer remained a challenge for a long time (this is a generic difficulty of the Bethe ansatz, which goes back to the Lieb-Liniger solution of the interacting Bose gas in one dimension [21]). An algorithm to do this was finally found by Volin in [22,23]. This opened the way to a resurgent analysis of quantum integrable systems. The presence of renormalons in the non-linear sigma model was tested numerically in [22]. In [24], this analysis was extended to many two-dimensional asymptotically free theories, and some classical predictions of renormalon asymptotics were verified in detail. A similar resurgent study of various non-relativistic models was carried out in [25][26][27]. Very precise results on the resurgent structure of the free energy in the O(4) sigma model have been recently obtained in [28,29], where Volin's method was used to generate the perturbative series up to very large order.
In this paper, we study the free energy of the O(N ) non-linear sigma model and its supersymmetric extension. Our first goal is a detailed test of the Bethe ansatz result against conventional perturbation theory. So far, this has been verified up to two-loops in [30]. To do a comprehensive test, we calculate the free energy to all loops, at next-to-leading order in the 1/N expansion. Our result matches the Bethe ansatz answer obtained in [22,24] to all available orders. In this way, it provides very convincing evidence for the validity of the Bethe ansatz, as well as for the techniques of [22,24] to extract perturbative series from it. It also provides additional evidence for the conjectural S matrices (although these have been tested directly, up to next-to-leading order in the 1/N expansion, in [31,32]). We should note that, to make contact with perturbation theory, our 1/N expansion is performed around the perturbative vacuum, and not around the non-trivial large N vacuum where particles get a non-perturbative mass already at tree level. In this sense, our calculation is very similar to what was done in [33]: we re-organize perturbation theory in powers of 1/N and keep the first two non-trivial contributions.
An interesting spin-off of our calculation is a new derivation of the beta function of the non-linear sigma model, at next-to-leading order in the 1/N expansion. This is a known result, going back to [34]. However, in [34] the beta function was derived from the epsilon expansion of the critical exponents (see [35] for a review). Our derivation is a direct one, similar in spirit and in the details to the computation of Palanques-Mestre and Pascual of the QED beta function in the limit of large number of fermions N f [36]. In the case of the supersymmetric non-linear sigma model, we verify the vanishing of the beta function at next-to-leading order in 1/N , which was also derived from the epsilon expansion of critical exponents in [37].
In some cases, one can calculate the 1/N expansion of the free energy directly from the Bethe ansatz equations [18,[38][39][40][41][42]. However, even if the Bethe ansatz result contains complete information about the groundstate energy, including both the perturbative expansion and non-perturbative corrections, at present no known method exists to extract analytically the nonperturbative sectors, not even at large N . In contrast, our calculation of the perturbative series at next-toleading order in 1/N is fully analytic. Using techniques developed in [27,33,43], we can extract the trans-series at large N from our all-loop results. This trans-series signals the presence of an isolated infrared (IR) renormalon singularity, 1 and an infinite sequence of ultraviolet (UV) renormalon singularities.
In the supersymmetric extension of the non-linear sigma model, an important consequence of our perturbative computation is an analytic proof, at next-toleading order in 1/N , that the IR singularity arising from bosonic diagrams does not cancel against the IR singularity arising from fermionic diagrams. Therefore, the supersymmetric model exhibits IR renormalons. 2 This is in contrast to the disappearance of leading IR renormalons in some supersymmetric theories, pointed out [44,45].
The organization of this paper is as follows. In Sect. 2, we review basic aspects of the non-linear sigma model, its 1/N expansion, and its Bethe ansatz solution. Section 3 presents the calculation of the effective potential at next-to-leading order in 1/N . This result is then used in Sect. 4 to extract the large N renormalons and their trans-series. In Sect. 5, we extend all these results to the supersymmetric non-linear sigma model. Finally, in Sect. 6, we present our conclusions and prospects for future developments. There are three Appendices with additional details and clarifications on our calculations.

The non-linear sigma model and its Bethe ansatz solution
The O(N ) non-linear sigma model is a quantum field theory in two dimensions for a vector field S(x) = (S 1 (x), . . . , S N (x)), satisfying the constraint The Lagrangian density is where g 0 is the bare coupling constant. The non-linear sigma model is asymptotically free [46], and it can be regarded as a toy model for gauge theories. It also has many different applications in condensed matter physics, where it is used to model the low-energy dynamics of statistical systems with a global O(N ) symmetry. We will write the beta function for the cou-pling constant g as With this convention, asymptotically free theories have β 0 > 0. All of our perturbative calculations will be done in the MS scheme. For the non-linear sigma model, the first two coefficients of the beta function are [47] β 0 = 1 4πΔ , The beta function is known up to four loops [34,48]. The non-linear sigma model can be also studied in the limit in which the number of components of S is large and relevant quantities are computed in a large N expansion. In this setting, one introduces the 't Hooft coupling which is kept fixed in the large N limit. As we will see, it will be more natural to make the expansion in powers of Δ, rather than in 1/N . The large N expansion makes it possible to obtain non-perturbative results for this model (see, e.g. [5]). In perturbation theory, one expands around an ordered vacuum with S = 0, in which the global O(N ) symmetry is spontaneously broken. This leads to a perturbative spectrum consisting of N − 1 Goldstone bosons. This can not be the case at the non-perturbative level, due to the Coleman-Mermin-Wagner theorem. Indeed, at large N , one finds a spectrum consisting of N massive particles in the fundamental representation of O(N ), which is thought to be the true spectrum of the theory. It is important, however, to keep in mind that perturbation theory around the "false vacuum" gives the correct asymptotic expansion of O(N ) invariant observables, as pointed out in [49,50].
Many quantities in the non-linear sigma model can be calculated systematically in a 1/N expansion, like for example critical exponents. Using this method, one finds that the beta function for the 't Hooft coupling, defined as has the following 1/N expansion [34]: where 10) and is related to the number of dimensions of the theory by As a spin-off of the results in this paper, we will rederive the result for β (1) (λ) by a direct calculation in perturbation theory. A conjectural expression for the exact S-matrix of the two-dimensional non-linear sigma model was put forward in [11,12]. This makes it possible the following exact computation [13]. Let H be the Hamiltonian of the model, and let Q be a conserved charge, associated with a global conserved current. Let h be an external field coupled to Q. h can be regarded as a chemical potential, and as usual in statistical mechanics we can consider the ensemble defined by the operator H − hQ. (2.12) The corresponding free energy per unit volume is then defined by where V is the volume of space and β is the total length of Euclidean time. This is the ground-state energy of the model in the presence of the additional coupling.
As pointed out in [13], we can compute using the exact S matrix and the Bethe ansatz. We will refer to (2.14) as the free energy. After turning on the chemical potential h beyond an appropriate threshold, there will be a density ρ of particles, charged under Q, and with an energy per unit volume given by e(ρ). These two quantities can be obtained from the density of Bethe roots χ(θ). This density is supported on an interval [−B, B] and satisfies the integral equation In this equation, m is the mass of the charged particles, and with a clever choice of Q, it is directly related to the mass gap of the theory. The kernel of the integral equation is given by (2.17) Let us note that B is fixed by the value of ρ, and this leads implicitly to a function e(ρ). Finally, the free energy can be obtained as a Legendre transform of e(ρ): (2.18) Note that the first equation defines ρ as a function of h.
The above program can be implemented in a number of models. In the case of the non-linear sigma model, one considers the conserved currents associated to the global O(N ) symmetry, We will denote by Q ij the corresponding charges. Usually [15,16], one considers in (2.12) the quantum version of Q 12 . The exact S matrix of the O(N ) non-linear sigma model, for particles charged under Q 12 , is given by and Δ is given in (2.5).
The perturbative series can be extracted from the Bethe ansatz solution with a method developed by Volin [22,23]. In his original work, the method was applied to the non-linear sigma model, but it was later extended to other quantum integrable models in [24][25][26][27]. It is convenient to use an expansion variable which can be connected directly to the perturbative answer. Such a variable was introduced in [30] and is defined by and Λ is the dynamically generated scale 3 A scheme for QCD closely related to (2.21) was studied in [51].
This scale is proportional to the mass m appearing in the Bethe ansatz. In the case of the non-linear sigma model, one has [15,16]: Λ. (2.24) In this way, we obtain a power series in α for the normalized energy density, The result, up to order α 4 , is The free energy δF(h) can also be calculated in perturbation theory from the effective potential. The coupling to the conserved charge Q 12 leads to the modified Lagrangian [15,16,30] (2.27) δF(h) was obtained at one-loop in [15,16], and at twoloops in [30]. The result of the perturbative calculation can be expressed in various convenient ways, depending on an appropriate choice of coupling. We can use the renormalization group (RG) to re-express the perturbative series in terms of the couplingḡ 2 (μ/h, g), defined by , (2.28) or in terms of the dynamically generated scale defined in (2.23). In terms ofḡ, the free energy is To compare with the result of the Bethe ansatz, we can use the Legendre transform of δF(h) to calculate e(ρ), and then re-express the result in terms of the coupling α introduced in (2.21). We note that, at leading order in the coupling expansion, One obtains in this way from perturbation theory, and up to two-loops [30], This is in agreement with the result of the Bethe ansatz calculation (2.26), as verified in [22,30]. It is convenient to organize the free energy as a 1/N expansion: (2.32) Similarly, the normalized energy density in (2.25) can be organized as a double power series expansion in α, Δ: e(ρ) where E (α) are power series in α. One finds, at leading order, and at subleading order in Δ we have, for the very first terms [22,24], (2.35) The series E (1) (α) has been computed analytically up to order 44 in [24]. In the next section we will compute E (1) (α) analytically and at all loops, directly in perturbation theory, and we will match the result (2.35) and up to order 44.

The 1/N expansion of the effective potential
To evaluate δF(h) in perturbation theory we have to calculate the effective potential in the theory with Lagrangian (2.27). As in other problems with an O(N ) symmetry, it is convenient to reformulate the model in terms of a linear sigma model, by including an additional field X which implements the constraint (2.1).
In this way, we consider We expand around the following classical vacuum: where σ and χ are constants that minimize the potential. Neglecting linear terms, the resulting Lagrangian can be organized as The tree-level Lagrangian is while the quadratic or Gaussian part is given by with η = (η 1 , . . . , η N −2 ). The interaction part is In the theory with h = 0, there is an ordered phase with χ = 0, σ = 0 in whichσ 2 , η are Goldstone bosons. This phase is not realized quantum-mechanically, as we explained in the previous section. However, once h is turned on, these η bosons acquire a mass. After writing the quadratic terms L G in momentum space and inverting the matrix for the fields (σ 1 ,σ 2 ,χ), we obtain the propagators

Effective potential at one loop
The effective potential is a function of the vacuum expectation values (vev) σ and χ, and the parameter h. It has a 1/N expansion in powers of Δ given by (3.8) The vevs σ and χ are obtained by extremizing the potential, and they can be also expanded in 1/N : The leading order vevs σ (0) and χ (0) , which are the only ones needed in our calculation, are obtained as At next-to-leading order in Δ, we have Let us first compute V (0) (σ, χ; h). It has contributions from the tree level Lagrangian L tree and from the 1/Δ = N − 2 fields η at one-loop: Using dimensional regularization to evaluate the integral, we find From here we can calculate σ (0) and χ (0) . There is a "disordered" non-perturbative vacuum with χ (0) = 0 and an "ordered" perturbative vacuum with χ (0) = 0. As in [33], to make contact with conventional perturbation theory, we choose the perturbative vacuum, where χ (0) = 0. Imposing (3.10), we find (3.14) We can now introduce the renormalized coupling λ by the usual equation, where parametrizes the scale choice μ and features the MS scheme, and γ is the Euler-Mascheroni constant. Z is the renormalization constant, for which we consider the 1/N expansion Cancelation of singular terms in the r.h.s. of (3.14) fixes Now, the r.h.s. of (3.14) is manifestly finite as → 0, and we find By evaluating the leading order effective potential (3.13) at the critical point σ (0) , χ (0) = 0, we obtain the leading order free energy δF (0) (h), defined in (2.32).
After writing the resulting expression in terms of the renormalized coupling λ and in the limit → 0, one finds (3.20)

Ring diagrams
As explained in [52] (see also [33]), the next-to-leading correction to the effective potential in the 1/N expansion is given by a sum of ring diagrams. A ring diagram at m loops is constructed with m bubbles of η particles successively connected by m propagators ofχ particles until the diagram closes on itself. Each bubble comes with a factor 1/Δ = N −2 (the number of particles contributing to the bubble), which cancels with the factor Δ coming from the pair of vertices that connect with theχ propagators. The sum of such ring diagrams is shown in Fig. 3. Following the structure used in [33], the contribution of these diagrams is the infinite sum (3.22) is the scalar polarization loop arising from the η bubbles.
The momentum integrals in (3.21) are divergent. However, going back to (3.13), after renormalization of the term with 1/λ 0 , the renormalization constant at next-to-leading order Z −1 (1) provides additional divergent terms that should cancel the divergences of our momentum integrals. By evaluating at the critical point with χ (0) = 0, we find the renormalized free energy at next-to-leading order (3.23) The integrals I m are defined as By first expanding (k 4 + 4h 2 k 2 0 ) m with the binomial theorem, the angular integral of the momentum can be computed term by term. I m can then be expressed as a finite sum of rotational invariant integrals. Finally, to compute the integral over k 2 , it will be crucial to write the scalar polarization function in terms of a hypergeometric function: In this equation, The integrals I m then can be expressed as where (x) n is the Pochhammer symbol, S d = 2π d/2 / Γ (d/2) is the volume of a d-dimensional sphere, and which we have appropriately expressed in terms of the integration variable z = x/(x+4). It is easy to see that, when ≥ 2, the integrals I m, are finite as → 0. Thus, for now we focus on the integrals with = 0, 1 and their singular part. It is convenient to re-express these integrals in yet another form, specially to extract the singular part. We use fractional linear transformations of the hypergeometric function to write Similar manipulations of the polarization loop were done in [53,54] to calculate critical exponents in the 1/N expansion. We can now write the integral I m, as

30)
In writing this expression, we have changed the integration variable from z to 1 − z. The integral (3.31) can be written as an infinite sum which will be useful for our analysis. Let us define the expansion coefficients Then, by expanding the binomial in (3.31) and using (3.32) for each of the hypergeometric functions arising in the binomial sum, we can integrate term by term in z, and we find: . (3.34) In Appendix A, we prove that the terms k = 0, 1 in (3.34) contain the Laurent expansion of B m, up to order m−1 . From the factor 1/ m in (3.30), we note that the expansion of B m, up to this order precisely corresponds to the singular part of I m, . Thus, the two terms B m, ,0 and B m, ,1 fully contain the singularities of I m, .

Divergences and beta function
The renormalized free energy (3.23) must be finite. By imposing cancellation of divergences we should be able to obtain an explicit expression for Z −1 (1) , and thus, for the next-to-leading order of the beta function in the 1/N expansion. This result, which we quoted in (2.10), has been known for some time [34], based on the 1/N calculation of critical exponents [53,54]. Our calculation provides a direct derivation of the beta function, very similar to the calculation by Palanques-Mestre and Pascual in [36], where they studied the beta function of QED in the 1/N f expansion.
Let us write the next-to-leading correction to Z −1 in (3.17) as a Laurent expansion with yet undetermined coefficients B (n) i : Firstly we will relate the coefficients B (n) i directly with the coefficients of the beta function and, secondly, we will find an explicit expression for B (n) i by imposing cancellation of divergences in (3.23) and using our result in (3.34).
The beta function for the 't Hooft coupling can be obtained from the renormalization constant Z −1 as (3.36) Using the Laurent expansion (3.35) and the leading order result (3.18) for the renormalized coupling, the above equality can be written at order Δ in the 1/N expansion as In addition, finiteness of the β function as → 0 requires that Thus, the computation of β (1) (λ) now reduces to determining the coefficients B (n) i . As we prove in Appendix A, the divergent part of the integrals I m, comes from a finite number of terms k in the sum of (3.33). To be more specific, B m, ,k leads to singularities only for k = 0, 1 when = 0, and for k = 0 when = 1. Thus, the divergent part of I m will be contained completely in the sum where the factor in front of the = 1 term arises from the binomial coefficient and Pochhammer symbols in (3.27). By combining appropriately the Γ factors in (3.34), one finds (3.41) By Taylor expanding f m (s , ) in powers of the first variable and using properties of the binomial coefficient, we can compute the sum in (3.40). We obtain We remark that the expression we have obtained is only valid up to order m , while the first term in the r.h.s. is equal to S m up to order m−1 (so it correctly encodes the singular part of the integrals I m ). We will, however, need the second term to extract the finite part of the free energy.
We are now ready to extract the singular part of the sum over ring diagrams in (3.23). First, let us introduce some notation. Let A( ) be a Laurent series in . We will denote by divA( ) the singular or principal part of the series. Then, it is easy to show that div and we have denoted The expression in (3.43) is obtained using (3.14) in place of λ 0 /σ 2 (0) and reexpanding in powers of λ 0 . Replacing S m in (3.44) by the first term of (3.42), we find div Π r (h; ) = div 1 (r + 1) r P ( , (r + 1) ) , (3.46) where We can now use a similar argument to the one in [36]. In the r.h.s. of (3.43), we replace λ 0 by the renormalized coupling at leading order in 1/N : We then obtain where we used that and we expanded Because the expansion functions P j (x) are regular at x = 0, it is obvious from (3.49) that only P 0 ( ) contributes to the singular part of the sum over ring diagrams. Using the reflection formula for the gamma function, we obtain the explicit expression which has a power series expansion at = 0 of the form: Requiring cancellation of divergences in (3.23) determines the expansion of (3.35), and we find the values which satisfy the constraint (3.38). The beta function at next-to-leading order in the 1/N expansion can be now computed by going back to (3.37) and using our result for the coefficients B (n) i : This result of course coincides with (2.10).

Finite part of the free energy and comparison with the Bethe ansatz
Once the singularities have been canceled, we can focus on the finite part of the renormalized free energy. The finite part arises from the integrals I m, with ≥ 2 and from the finite part of the integrals I m, with = 0, 1. This last contribution comes from three types of terms. Two of them are already written down in the previous section, since they have their origin in S m : the last term in the r.h.s. of (3.42), and the finite part in the r.h.s. of (3.49). In addition, we have to take into account the contribution from the terms which are not included in (3.39). It will be convenient to rewrite the series (3.57) in the following integral representation, which we derive in Appendix A: The derivation of this result relies on the same tricks that we used to obtain (3.42) (similar manipulations can also be found in [55]). We can now write the next-to-leading correction to the renormalized free energy as Many of the ingredients appearing in these formulae have been already spelled out in detail. The coefficients P 0,m can be read from (3.52). The coefficients P m (0) follow from (3.47) and (3.51): The function f m (y, x) is given in (3.41). It remains to compute the integrals B m,0 , B m,1 and I m, , = 2, . . . , m, for arbitrary m, . This can be done analytically, and the results are presented in Appendix B. This allows us to compute δF (1) (h) at any given order. The very first terms read (3) (3)) log (h/μ) 32 (3.63) To compare with the Bethe ansatz solution, we have to re-express this result in terms of the coupling constant α, defined in (2.21). The first step is to set λ to the running coupling constant at the scale μ = h, which defines λ = λ (μ = h). λ is related to the dynamically generated scale Λ and h by where ξ is defined in (2.22). At this scale, δF (1) (h) simplifies greatly to (3.65) This defines the coefficients v m , m ≥ 1. We can now use the Legendre transform (2.18) to obtain the normalized energy density (2.25). To do that, it is useful to introduce yet another coupling which was first considered in [30] and is related toλ bȳ (3.67) We can useα to write the free energy as (3.69) The Legendre transform gives (3.70) The final step is to relateα to the coupling α defined in (2.21): which leads to a remarkably simple expression for the normalized energy density: (3.72) Expanding β (1) (u) with (3.56), we can read the following result for the series E (1) (α), defined in (2.33), in terms of the perturbative coefficients v m : (3.73) It follows from this expression that where the functions in the r.h.s. were defined in (3.61).
One can now compare the expression (3.73) with the result of the Bethe ansatz (2.35). We find perfect agreement up to order α 44 .
As a side remark, up to order α 44 , we notice that every coefficient of E (1) (α) is the sum of a rational number plus a linear combinations of odd Riemann zeta functions. 4 In Appendix B we prove to all orders that Z(α; 1) can be written as linear combination of ζ(2k + 1) (in fact, with no rational term). On the other hand, X(α; 1) and Y (α; 1) do not satisfy this transcendentality property when alone, but the combination X(α; 1) + Y (α; 1) indeed does, up to order α 44 (in this case, a rational term has to be included with the linear combination of odd zetas), although we do not have a proof of this statement to all orders.

Large N renormalons and their trans-series
In [22,24], numerical evidence was given for the factorial growth of the perturbative series for (2.31), at fixed N . This was interpreted as a signature of renormalons [1][2][3]. In [24], the contribution of UV and IR renormalons was disentangled, and detailed evidence was given that the large order behavior of the perturbative series is in agreement with the predictions of renormalon physics. In particular, it was shown that the next-to-leading behavior of the asymptotics involves the first two coefficients β 0 and β 1 of the beta function. The evidence for these effects was based on a numerical study of the perturbative series and it focused on the leading singularities in the Borel plane.
At large N , the ring diagrams studied in the previous section should give the leading renormalon behavior. One advantage of having explicit results for these diagrams is that we can obtain from them analytic results on the large order behavior of the perturbative series. Equivalently, we can find explicit results for the exponentially suppressed trans-series associated to each Borel singularity. These can be obtained without even calculating the loop integrals. As shown in [27,33,43], it is enough to write the generating functions (3.61) in integral form and study their imaginary parts (or, equivalently, their imaginary discontinuities).
We will now present the integral forms for the series X(λ; 1), Y (λ; 1) and Z(λ; 1) appearing in (3.74). We note that the coefficients of the series W (λ; 1) do not grow factorially. This is easily observed from (3.62) and the fact that P 0,m are the Taylor coefficients of an analytic function at = 0.
Let us start with X(λ; 1). The integral form in this case is easily obtained using the explicit expression (3.59) and Laplace transforms. We find In the case of Y (λ; 1) we use the explicit expression (3.41) and write the Euler beta functions that appear in the resulting expression as integrals over z. In this way, we find We can use again Laplace transforms to sum up this series and we eventually find and E 1 (z) is the exponential integral. Finally, after using the identity the last series can be written as where Z(z, λ) = 1 + F (z)λ 1 + zF (z)λ , The advantage of the representations (4.1), (4.4) and (4.7) is that they lead to explicit, exponentially small imaginary terms. These are precisely the trans-series associated to the renormalon singularities.
Let us first consider the function X(λ; 1). It has discontinuities when λ < 0, due to the poles and logarithmic branch cut in the integrands of (4.1). The singularities of the integrand are located at z 1 and z 2 , which are defined by We find disc X(λ; 1) = 2πi 1 π z1 z2 dz z 2 − 1 π This discontinuity can be computed term by term as a trans-series in λ, i.e. as a power series in both e 2/λ and λ. See Appendix C for details of this computation. We find  Let us now consider the function Y (λ; 1). By investigating the function (4.5), we see that Y (λ; 1) has discontinuities both for positive and negative λ. When λ < 0, there is a pole at z = z 1 , where z 1 was defined in (4.9), and a discontinuity due to the exponential integral. For λ > 0, we have a discontinuity due again to the exponential integral, and one finds Finally, we consider the discontinuity of Z(λ; 1), which arises from two sources. The first one is a pole of the integrand, which appears for λ < 0. This occurs at a z 3 satisfying There is another source of discontinuity due to the square root inside the logarithm, which occurs when Z(z, λ) < 0. The discontinuity of this source is given by where (z 4 , z 3 ) is the subinterval of (0, 1) where Z(z, λ) is negative. The value z 3 is the pole previously discussed in (4.13), while z 4 satisfies The trans-series obtained in this way is  We can now put all these results together and calculate the trans-series associated to E (1) (α). It is given by 1/(2i) times the discontinuity, and reads  Let us analyze this result. The first term in the r.h.s. of (4.17) corresponds to a Borel singularity on the positive real axis at ζ = 2. It is an IR renormalon, which has been identified in [22,24] at finite N and more recently in [28,29] in the O(4) model. At this order in the 1/N expansion there are no additional IR renormalons. The next terms correspond to singularities in the Borel plane on the negative real axis, and they are UV renormalons. We conjecture that they are located at The derivation follows from the general theory of resurgence (see e.g. [5,7]), which relates the trans-series of a function to the singularities of its Borel transform. In turn, one can extract the large order behavior of the perturbative coefficients directly from the Borel singularities. We find (4.20) The first term in the r.h.s. gives the leading order asymptotics, which is due to the first UV renormalon singularity at ζ = −2. The second term, with fixed sign, is due to the IR renormalon singularity at ζ = 2; while the third term, which is exponentially subleading with respect to the first two, is due to the next UV renormalon at ζ = −6.
We can now compare the large order behavior (4.20) with the results of [24] for finite N . In general, one has to exercise care in this comparison, since some of the contributions to the asymptotic behavior detected at finite N might be suppressed at large N . In [24], it was found that the series (2.25) at finite N has the following asymptotics: where the first and second term in the r.h.s. correspond respectively to the IR and UV singularities at ζ = ±2, and C IR , C UV are in principle functions of N . After taking into account the difference in the labelling, m = n+1, we find that (4.21) leads to the behavior (4.20) at large N , since Δ → 0 in the arguments of the Gamma functions. This also implies that C IR , C UV are of order Δ.
In [24] it was argued that the IR renormalon at ζ = 2 is associated to the condensate of the operator O = ∂ μ S · ∂ μ S, and that the large order behavior (4.21) is compatible with the expectations of renormalon physics. It follows that our large N result for this IR renormalon might be also explained by the contribution of this condensate.
The calculation above determines the functional form of the trans-series associated to the different singularities. An additional, interesting question is the "semiclassical decoding" of the normalized density (2.33) at large N , i.e. its expression as a Borel-Ecalle resummation of these trans-series and the perturbative series. At finite N = 4, this has been recently done in [28,29] for the full energy density, using trans-series at finite N . Results along this direction in the large N expansion will be reported in [42].

The supersymmetric non-linear sigma model
We can extend all of the above results to the supersymmetric version of the non-linear sigma model considered in [56]. This model consists of the vector field S(x) of the purely bosonic version, satisfying as well the constraint (2.1), and an N -tuple of two-component Majorana spinors Ψ = (Ψ 1 , . . . , Ψ N ), satisfying the constraint The Lagrangian density is where we follow the conventions of [56] for the gamma matrices. This model is asymptotically free and its beta function is of the form (2.3), with and Δ is again given by (2.5). The model can also be studied in the large N expansion, where one finds a non-perturbative mass gap and dynamical breaking of the discrete chiral symmetry Ψ → γ 5 Ψ [57]. The beta function in the 1/N expansion has the structure (2.8), where β (0) (λ) is given again by the expression in (2.9), while β (1) (λ) = 0 due to cancellations between bosons and fermions [37]. We will rederive this result in Sect.

5.2.
As in the non-linear sigma model, we couple the present model to an external potential using again the conserved charge Q 12 associated to the global O(N ) symmetry. The dependence of the ground-state energy on the external potential can be obtained from the Bethe ansatz and the exact S-matrix conjectured in [58]. The resulting integral equation was written down explicitly in [19], where it was used to obtain the exact mass gap of the model. In [24] the ground-state energy was computed as a power series in the coupling α, defined in (2.21) (although ξ = 0 in this case). At leading order in the 1/N expansion, one obtains for E (0) (α) the same result that we presented in (2.34). At next-toleading order in 1/N one finds which is available up to order α 42 in [24]. Interestingly, this expansion is almost identical to the bosonic result (2.35), but it only keeps its transcendental part. Our goal in this section is to test the result (5.4) against a perturbative calculation. Like before, it is more convenient to use the linearized version of the model, which is obtained by introducing three auxiliary fields: one scalar field X to impose the constraint (2.1), a Majorana fermion λ to impose the constraint (5.1), and a Hubbard-Stratonovich scalar field τ to integrate out the quartic fermionic term in the Lagrangian (5.2). The resulting Euclidean Lagrangian, which includes the coupling to the chemical potential h for both S and Ψ fields, is given by 5 5 In Euclidean signature, the gamma matrices satisfy the anti-commutator relation {γ μ , γ ν } = 2η μν E I, where η μν E is the Euclidean metric. After a Wick rotation of the action, the Euclidean gamma matrices γE that enter in the Euclidean Lagrangian We now expand the Lagrangian around the following classical vacuum, (5.6) Up to linear terms in the fields, one finds that the Lagrangian can be written as The tree-level Lagrangian is given by The quadratic terms are given by where ψ = (ψ 1 , . . . , ψ N −2 ). Finally, the interaction terms are given by From (5.9) we can compute the propagators in momentum space. The propagators of the boson fields η,σ 1 ,σ 2 andχ were already obtained in (3.7). There is, however, an additional bosonτ with propagator The fermion propagators are (5.12) The propagators of the bosonic and the fermionic fields are represented diagrammatically as in Figs. 5 and 6. The interaction terms are represented by the vertices in Fig. 7.
We can now calculate the effective potential in an expansion in powers of Δ, as in (3.8). The leading order term comes from the tree-level Lagrangian together with the one-loop contributions of the η bosons and ψ fermions, for which there are 1/Δ = N − 2 of each: (5.13) Using dimensional regularization, we obtain (5.14) The extremization procedure for χ and σ is identical to the bosonic case. As for τ , much like χ, there is the non-perturbative choice τ = 0 and the trivial one τ = 0. We choose the latter to connect our result with perturbation theory. After renormalizing the coupling like in the bosonic case, we find

Ring diagrams
There are three types of ring diagrams that contribute to the effective potential at next-to-leading order in Δ: ring diagrams with bosonic η bubbles, ring diagrams with fermionic ψ bubbles, and ring diagrams with mixed ψ-η bubbles. The first type of ring diagrams are the same ones of the bosonic sigma model, so we do not have to compute them again. Since the propagators for ψ andτ are both h independent, fermionic ring diagrams do not contribute to δF(h) after subtraction at h = 0. Thus we only need to consider mixed ring diagrams.
The ψ-η bubbles are connected by λ lines. Then the contribution of mixed ring diagrams to the effective potential is A minus sign arises from the single fermion loop that runs across the entire diagram. This sign gets canceled by an additional sign that appears in the computation of the effective potential. The polarization loop is in this case This integral can be computed with standard Feynman techniques and, after evaluation at the vacuum χ (0) = τ (0) = 0, we obtain Again, as in the bosonic case, we have expressed the polarization loop in terms of a hypergeometric function and the variable x = k 2 /h 2 . The free energy at nextto-leading order in Δ from bosonic plus mixed ring diagrams can now be written as The integrals I m , corresponding to bosonic diagrams, were already defined in (3.24), and we have a set of new integrals from the mixed diagrams given by The front factor in (5.19) has been extracted for better comparison with the non-supersymmetric result of (3.23). As we already mentioned, one important difference in the present model is that we do not need a renormalization constant to cancel the divergences of the ring diagrams. Instead, there is a total cancellation of divergences between bosonic and mixed diagrams.
We will see this explicitly in Sect. 5.2, thus proving that the beta function at subleading order is β (1) (λ) = 0. At this stage, it is convenient to compute the trace in (5.20). Expanding with the binomial theorem, we obtain

c) Mixed
We can calculate the trace recursively in arbitrary dimension d: As is the standard procedure, we take the dimension of the spinor representation to be a fixed integer 2 d 2 = 2. In d-dimensional spherical coordinates we can always pick k 0 = k cos θ 1 . Then we obtain Tr (γ 0 / k) 2 = 2k 2 cos(2 θ 1 ), (5.23) and the momentum angular integral can now be computed term by term in (5.21), yielding where we singled out the integrals For ≥ 2 the integrals are finite at = 0. Since the factor 1/Γ (d/2 − ) vanishes in the sum of (5.24) for = 0, none of the terms with ≥ 2 will contribute to the free energy.

Cancellation of divergences
The goal in this section is to check that the subleading free energy (5.19) in the supersymmetric model is already finite without the need of renormalization. For that, we follow similar techniques to those in Sects. 3.2 and 3.3. That is, we want to isolate the singular part of the integrals I m, , as we did in (3.33)-(3.34). We start by expressing the hypergeometric function in (5.25) as a sum of two terms, using fractional linear transformations: We now go back to (5.25), plug in (5.26) and change the variable of integration from z to 1 − z. We get (5.28) By expanding the square bracket with the binomial theorem and integrating term by term with the Euler beta function, we obtain (5.30) The coefficients d (s) k are defined by the Taylor expansion The Laurent expansion of B m, up to order m is obtained by summing only the terms k = 0, 1 for = 0 and k = 0 for = 1. This follows from a computation similar to the one in appendix A. We then find This result is slightly different to the bosonic case (compare (5.32) to (A.1)), and it greatly simplifies the computation of mixed ring diagrams, because the terms in (5.32) already incorporate the singular part plus the finite part of I m, . Equation (5.30) gives a divergent result for m = 1, = 1 (even for arbitrary ). In this special case we have to compute the integral (5.25) explicitly. Instead of using the hypergeometric representation, we use the integral representation of (5.18) and commute the x and y integrals. Then we can compute both integrals analytically in terms of the Euler beta function B(x, y). We obtain, We are now ready to compute the divergent part arising from mixed ring diagrams in (5.19). We start by considering the sum which can be written in the form (5.36) We can now perform the sum in s, as we did in (3.42), and we obtain the result (5.37) The first term in the r.h.s. is equal to S m up to order m−1 , so that term is sufficient to compute the singular part of the integrals I m . The second term, contributing to order m , will be needed for the computation of the finite part of I m . At this point, it becomes a simple exercise to verify that the singular parts of I m and I m cancel each other in (5.19). We first observe that If we now rewrite the sum of the terms = 0, 1 appearing in (3.27) and (5.24) using S m and S m , respectively, one can verify from the observation above that and mixed ring diagrams. In addition, this shows that β (1) (λ) = 0 from a direct computation in perturbation theory.

Finite part of the free energy and comparison with the Bethe ansatz
We can also organize the finite parts of the mixed ring diagrams contributing to the free energy in a manner similar to what we did in (3.61) for bosonic ring diagrams. We note, however, that there are a few simplifications in this case. The analogue of Z, corresponding to the sum of the integrals I m, for ≥ 2, is equal to 0. The analogue of X, corresponding to the contribution from the terms we missed in the sum (5.34), is also 0, as is easily inferred from the result (5.32). Moreover, the analogue of W is present in the mixed ring diagrams, but it exactly cancels with that of the bosonic diagrams (this can be deduced from the derivation that led to (5.38)).
In short, we only need to consider the analogue of Y , which arises from the second term in the r.h.s. of (5.37). After simplifying f m (y, 0) we find (5.40) and the computation of the (m + 1)th derivative follows naturally from the geometric series. Then, the contribution to the free energy is simply (5.41) Putting all our results together we obtain Thus, the ground-state energy is E susy (1) (α) = −2πα 2 X (α; 1) +Y (α; 1)+Z (α; 1)+Y (α; 1) . (5.43) Here, α is defined by Eq. (2.21) with β 0 given in (5.3) and ξ = 0. We find perfect agreement between (5.43) and the coefficients obtained from the Bethe ansatz, which were calculated up to order α 42 in [24]. As we previously observed in (5.4), the perturbative expansion of E susy (1) (α) does not contain any rational term (at least to the available order). In our perturbative computation, we verified that Y(α; 1) cancels the rational part of X (α; 1) + Y (α; 1), so that the coefficients are linear combinations of Riemann zeta functions evaluated at odd arguments.
As we did in the bosonic case, we can obtain from these results the location of the renormalon singularities and their associated trans-series. The mixed diagrams (5.41) contribute only to the IR renormalon singularity at ζ = 2, and we find where only the sign of the IR term differs from the bosonic case. The asymptotic behavior extracted from this discontinuity matches the coefficients with the expected precision. Let us note that, although divergences cancel between bosonic and mixed ring diagrams, the IR renormalon at ζ = 2 does not cancel. This is in contrast to the cancelation of leading IR renormalons that occurs in some supersymmetric theories according to [44,45].

Conclusions
The Bethe ansatz calculation of the free energy in twodimensional integrable models is one of the most interesting exact results in quantum field theory. It makes it possible to understand quantitatively many important aspects of asymptotically free theories, like dynamical mass generation and the presence of renormalons and condensates. For this reason, it is important to test the predictions of the Bethe ansatz against more conventional methods in quantum field theory. This has been done in the past, in particular in the case of the Gross-Neveu model [17,18]. However, a direct test against perturbation theory was limited until very recently by the difficulty of extracting perturbative series from the Bethe ansatz. The results of [22][23][24] have opened the possibility of such a test, and this has been the first goal of this paper. We have performed a calculation of the free energy at next-to-leading order in the 1/N expansion and at all loops, and we have verified the predictions of the Bethe ansatz up to very high order in the coupling constant, both in the non-linear sigma model and its supersymmetric extension. The second goal of this paper has been to use this analytic, all-orders result to obtain detailed information about renormalon singularities and their associated trans-series, at leading order in the 1/N expansion.
There are many problems open by our investigation and by closely related efforts. First of all, one could extend the tests presented here to other integrable field theories. In [24], explicit results for the perturbative series of the free energy have been obtained for the Gross-Neveu model and the principal chiral field with different choices of chemical potentials. There are also results for integrable, non-relativistic models, like the Gaudin-Yang model [25]. All these perturbative expansions, obtained from the Bethe ansatz, could be tested by conventional techniques, and these tests would provide additional insights. It should be mentioned that, in the case of the principal chiral field, which is a matrix model, a direct perturbative calculation at higher loops is more challenging than for the vector models that we studied in this paper. Another interesting problem is to extend our results to higher orders in the 1/N expansion, even for the models considered here.
In our calculation of the free energy we did an expansion around the trivial large N vacuum, to make contact with perturbation theory. However, it is well known that both the non-linear sigma model and its supersymmetric version have a non-trivial large N saddle point, which leads to a non-perturbative mass gap. In fact, previous analysis of renormalons in the non-linear sigma model have been based on expansions around the non-trivial large N saddle point [59][60][61]. A preliminary calculation of the free energy considered in this paper around this non-trivial vacuum seems to lead to a purely non-perturbative result, at next-to-leading order in the 1/N expansion. It would be very interesting to understand more precisely the relationship between the expansions around these two very different vacua, and the rôle of renormalons in each of them. Note that similar issues are raised in the two-dimensional linear O(N ) model. There, the ground-state energy can be calculated in the trivial large N vacuum, similar to what we did in this paper [33], but it can be also calculated in the non-trivial large N vacuum [62].
The exact Bethe ansatz solution contains in principle all the information in the problem, both perturbative and non-perturbative. One of the goals of the resurgence program is to "semiclassically decode" this exact answer, by writing it as a Borel-Écalle resummation of a trans-series. Substantial evidence that this can be done was obtained recently in [28,29], in the case of the O(4) sigma model, by an impressive calculation. Additional evidence has been given in the context of integrable many-body systems, in [25,27]. The first step in the semiclassical decoding is to obtain the explicit form of the trans-series. So far, in all the problems solved by the Bethe ansatz, this has been done by looking at the large order behavior of the perturbative sector. A more challenging problem is to extract the transseries directly from the integral equation defining the free energy. This would require an extension of Volin's method by explicitly including the exponentially suppressed corrections.
Another important open question is to provide a physical interpretation of the trans-series that we have obtained. It was argued in [24] that the IR renormalon at ζ = 2 in the non-linear sigma model might be explained by the condensate of the operator O = ∂ μ S · ∂ μ S. Is it possible to devise some sort of perturbation theory in the background of the condensate which allows us to reproduce analytically the transseries? Such a generalized perturbation theory, akin to the one used in the calculation of OPEs [63], would provide one of the missing ingredients in our understanding of quantum field theory.
Funding Open Access funding provided by Universitè de Genève.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Integral representation of the series B m,
In this appendix, we will derive the integral representation of B m, that we presented in (3.58). Because the result has a Laurent expansion starting with m , this will also prove that the integrals (3.31) for = 0, 1 satisfy where the coefficients B m, ,k are given in (3.34). In particular, this means that the singular part of the integrals Im can be extracted from the combination in (3.39).
We start from (3.31) and expand the binomial. We then subtract the terms k = 0, 1 (for = 0) or k = 0 (for = 1) to obtain the following expression for the sums in (3.57) where g is given in (3.41), and Φ(x, z) is defined by We now Taylor expand g (y, x; z) in its first argument The binomial identity The O m+1 contains further terms of the Taylor expansion that we have ignored. Because the above integral is convergent at = 0, we can set = 0 in the integrand of (A.7) (further corrections will be of order m+1 ). This completes the proof of (A.1) and yields the result in (3.58), where g (y; z) ≡ g (y, 0; z).

B Analytic computation of the integrals
In this appendix we derive and present analytic results for the integrals I m, and B m, appearing in (3.61). Let us start with the integrals B m, , defined in (3.58). One way to compute them is to obtain a closed expression for the generating functional .

(B.2)
This expression is not suited for a Taylor expansion around y = 0 (because of the hypergeometric function), but it can be put in a more appropriate form by observing that

(B.3)
We then use Theorem 3 in [64] and expand the singular gamma functions around their poles, which leads to the expression where ψ (0) is the polygamma function. After some additional massaging, we obtain the following alternative expression for (B.2), in terms of Pochhammer symbols: Let us now consider the integrals I m, , defined in (3.28), with 2 ≤ ≤ m. We first focus on the particular case = m: Im,m. Knowing that the integrals converge without the need of regularization, we can set = 0, and the hypergeometric function in the integrand becomes

C Asymptotic expansions of the discontinuities
To obtain the discontinuity of X(λ; 1), as presented in (4.11), we first have to express the residues in (4.10) in terms of the position of the poles z1 and z2. This can be easily done with L'Hôpital's rule. Then it suffices to compute z1 and z2 perturbatively in powers of e 2/λ . These perturbative solutions can be conveniently obtained by defining variables v and ξ such that
For the discontinuity of Z(λ; 1), it is convenient to change variables as z = (1 − 2ξv) 2 . (C.7) From (4.13) we find that where v2 is the same as in (C.2) and (C.4). The asymptotic approximation to the residue of the pole in (4.7) now follows in the same way as for the pole of Y (λ; 1). Meanwhile, from (4.15) we obtain a much more simple equation for v4, v4 = e −2ξ log ξv 4 (1 − ξv4),

(C.12)
It is useful to introduce one further variable v = 1 + ξw, (C. 13) and vi = 1 + ξwi. The advantage of using A2, A4 and w is that we can expand the integrand of (C.12) as a power series in ξ, and we find the following structure: i≥0 P (i) l (w, w2, w4) log where P (i) l,r are polynomials. For example, up to order O ξ 2 , we find (C.15) The coefficients of the ξ series in (C.14) are functions of w which are integrable on the interval [w4, w2]. After integration, we replace w2 and w4 by their respective perturbative solutions through (C.4), (C.10) and (C.13). Summing up the result of this computation with the residue of the pole at z = z3 leads to (4.16).