Resurgence and $1/N$ Expansion in Integrable Field Theories

In theories with renormalons the perturbative series is factorially divergent even after restricting to a given order in $1/N$, making the $1/N$ expansion a natural testing ground for the theory of resurgence. We study in detail the interplay between resurgent properties and the $1/N$ expansion in various integrable field theories with renormalons. We focus on the free energy in the presence of a chemical potential coupled to a conserved charge, which can be computed exactly with the thermodynamic Bethe ansatz (TBA). In some examples, like the first $1/N$ correction to the free energy in the non-linear sigma model, the terms in the $1/N$ expansion can be fully decoded in terms of a resurgent trans-series in the coupling constant. In the principal chiral field we find a new, explicit solution for the large $N$ free energy which can be written as the median resummation of a trans-series with infinitely many, analytically computable IR renormalon corrections. However, in other examples, like the Gross-Neveu model, each term in the $1/N$ expansion includes non-perturbative corrections which can not be predicted by a resurgent analysis of the corresponding perturbative series. We also study the properties of the series in $1/N$. In the Gross-Neveu model, where this is convergent, we analytically continue the series beyond its radius of convergence and show how the continuation matches with known dualities with sine-Gordon theories.


Introduction
Since its discovery [1,2], the 1/N expansion has played an important rôle as a tool to study nonperturbative aspects of quantum field theory. Many phenomena which are invisible in ordinary perturbation theory, like spontaneous chiral symmetry breaking or the dependence on the theta angle, can be discovered already in the large N limit of four-dimensional gauge theories [3,4].
In two-dimensional models one can even obtain explicit quantitative results at large N for the chiral condensate [5] or the topological susceptibility [6].
The reason that is often given for these successes is that the 1/N expansion "resums" the perturbative series, and therefore it goes beyond what is available in perturbation theory. However, making this statement precise requires being more explicit about what we mean by resummation. It is known since the earlier work [7,8] that the factorial growth of diagrams in perturbation theory is tamed to just an exponential growth at the planar level in large N matrix model QFTs. A similar phenomenon applies to higher order in 1/N and QFTs based on vector models. The reduced number of diagrams at fixed order in 1/N leads in many theories to expressions that are analytic functions at the origin of the fixed large N coupling constant, order by order in 1/N . 1 Theories of this kind include zero-dimensional (0d) matrix models, N = 4 super Yang-Mills theory in 4d, 3d O(N ) models, or Chern-Simons-matter theories.
In many theories, however, the coefficients in the perturbative series still grow factorially after restricting oneself to a fixed order in 1/N . This happens when the growth is not dominated by the proliferation of diagrams, but by integration over momenta. The 1/N expansion can tame the first growth, but not the second one. This is the phenomenon of renormalons (see e.g. [10] for a review). Therefore, the way the 1/N expansion resums the perturbative series in a theory with renormalons must be very different from what happens in the theories without renormalons that have attracted more attention.
The most general framework to resum perturbative series is the theory of resurgence, which has been extensively studied in recent years (see e.g. [11,12] for a review and references). In this theory, conventional perturbative series have to be extended to more general objects called trans-series, which include exponentially small corrections. This trans-series can be obtained from the perturbative sector by a detailed study of the singularities in the Borel plane. There is growing evidence that many exact quantities in quantum theory can be obtained as Borel resummations of these trans-series. These include energy levels in quantum mechanics [13][14][15][16][17][18], 1/N expansions in matrix models [19,20], and perturbative expansions in some quantum field theories [21][22][23][24][25][26][27].
We can now ask the following question. In a theory which admits a 1/N expansion and has renormalons, each order in the 1/N series is a non-perturbative function which resums perturbation theory. Can we decode each of these functions in terms of the conventional perturbative series, plus its associated trans-series? Can we in principle recover each of these non-perturbative functions from the perturbative series? In other words, what is the interplay between the resurgent structure of perturbation theory and the 1/N expansion?
An additional reason to ask this question is the following. The resurgent structure of fullyfledged quantum field theories is quite intricate. We know however that quantum field theories 1 By "expressions" we refer here to physical, and hence renomalization scheme independent, quantities. At each order in 1/N , the analyticity properties in the coupling of unphysical quantities, such as beta-functions, depend on the scheme. For instance, in the limit of large number of flavours n f , the first orders in 1/n f of 4d QED β-functions are analytic in the MS scheme, while they are non-analytic in other schemes [9]. tend to become simpler and more tractable in the 1/N expansion. One can then hope that by looking at the large N limit one will find somewhat simpler resurgent structures which can be studied analytically.
Another, more difficult question concerns the nature of the 1/N expansion itself. It is well-known that, even in zero-dimensional models, this expansion grows factorially or doublyfactorially, and the resummation and resurgent properties of this expansion have been studied in detail in toy theories, like matrix integrals (see [11] for a review and references). There has been much less progress in quantum field theory, due among other things to the difficulty of going to large order in the 1/N expansion.
In order to address these questions as concretely as possible, it is convenient to look at models which have renormalons and at the same time can be studied in detail, both in perturbation theory and in the 1/N expansion. The ideal candidates for such a study are asymptotically free theories in two dimensions which are integrable, i.e. their S-matrices are known exactly. In this paper we will focus on three classical examples: the O(N ) non-linear sigma model (NLSM), the SU (N ) principal chiral field (PCF), and the O(N ) Gross-Neveu (GN) model. 2 It was noted long ago by Polyakov and Wiegmann that a Thermodynamic Bethe Ansatz (TBA) can be used to compute exactly the free energy of these theories in the presence of an external field coupled to a conserved current [31]. In addition, when the external field is large, one can use asymptotic freedom to calculate this observable in perturbation theory, and this was exploited in [32][33][34][35][36][37][38][39] to obtain the relation between the mass gap and the dynamically generated scale (see [40] for a review). In addition, a powerful method developed in [41,42] makes it possible to extract the perturbative series for the free energy at very high orders. This has led to many quantitative studies of renormalon physics and resurgence in relativistic [25,26,43] and nonrelativistic [22,23,44,45] integrable quantum field theories.
These quantum integrable models have been also studied in the 1/N expansion [36,[46][47][48][49][50]. Our aim in this paper is to elaborate on and extend this line of research in order to answer in detail the questions raised above. Before diving into all the details of the 2d QFT models, we start in section 2 by considering the 0d reduction of certain large N vector models. We will show explicitly that each order in 1/N is analytic in the 't Hooft coupling, the large N expansion is factorially divergent, while the reduction to 0d of the free energy F(h) defined in (3.3) turns out to be analytic at N = ∞. In section 3 we come back to field theory and introduce the key observable we will consider in this paper, the free energy F(h) as a function of a chemical potential h. We review how this can be computed by using the TBA in 2d integrable QFTs. The main results of the paper are reported in sections 4, 5, and 6.
In section 4 we consider the NLSM. We compute F(h) at the leading and next-to-leading order in the 1/N expansion, which resums an infinite number of renormalon diagrams appearing in ordinary perturbation theory, described in detail in [50]. We extract the exact answer both from a direct QFT calculation and from the TBA equations. At this order in 1/N there is a single IR renormalon singularity and the exact answer is obtained by the so-called median Borel resummation of the perturbative series. We have also studied the properties of the 1/N expansion of F(h) by exploiting the TBA to generate several terms. Our explicit results do not show factorial growth and are inconclusive. Either the asymptotic regime has not been reached yet or the 1/N series is actually convergent. A similar analysis has also been made in the PCF model, with the same inconclusive result.
In section 5 we consider the PCF model. We find a new, explicit solution for F(h) at leading order in 1/N from TBA, and for the choice of charges used in [34]. The exact answer can be understood as the median resummation of a non-trivial trans-series which can be obtained analytically and has an infinite number of IR renormalon singularities (in contrast to the solution of [46,47], which has a single IR renormalon singularity, and similar to the numerical results obtained in [25,26] for the O(4) NLSM). Therefore, in this case the large N limit provides an explicit, analytic, yet non-trivial example of resurgence and median resummation in a model with infinitely many IR renormalon corrections.
In section 6 we consider the GN model. In this case, a new phenomenon appears: at each order in the 1/N expansion, the exact answer includes an infinite number of non-perturbative corrections. While ambiguities in imaginary terms nicely cancel between one series and the next in the trans-series, as expected from resurgence, real non-perturbative corrections can not be obtained from the resurgent properties of the perturbative series. Therefore, in this case there is a tension between resurgence and the 1/N expansion. The 1/N series of F(h) in the GN model turns out to be convergent, with a finite radius of convergence. Using the TBA, we generate many terms in the 1/N series. The latter can be analytically continued beyond its radius of convergence. Interestingly, the analytic continuation of this series gives reliable results for small N , such as N = 4 and N = 2, which are in agreement with the well-known dualities between these models and sine-Gordon theories.
In section 7 we conclude with a detailed discussion of our findings in the more general context of the theory of resurgence, and we present various directions for future work. We report in appendix A a few technical details needed to reproduce some results of section 2. In appendix B we present the explicit form of the kernels entering the TBA equations (3.4) and discuss their analyticity properties in 1/N .
2 Ordinary integrals at large N Before analyzing the 2d QFT models it is useful to consider ordinary integrals, where we can get complete analytic results and show the generic divergent nature of the 1/N perturbative series.
A notable example is given by the 0d reduction of the large N quartic vector models, given by (2.1) with m ∈ R and g > 0. We can trivially rescale m, so we get three different cases: m = 1, m = −1, and m = 0. The integral in (2.1) can be computed analytically, but we won't need its exact expression. The large order behavior of the 1/N expansion can be obtained by using steepest descent methods, see appendix A for details. We have whereÎ c , ρ and θ are explicit functions of g and m reported in (A.8) and (A.9). As discussed in the appendix A, the 1/N expansion of (2.1) is divergent asymptotic and Borel resummable for any real value of m and g > 0. This result should be contrasted with what we would get by expanding in g at fixed N. In this case, taking N to be odd, we can use radial coordinates with radial variable r, so that in a g expansion the relevant saddle points of the integral (2.1) are those of the function f (r) = m 2 r 2 + 1 4 r 4 . (2.5) The qualitative and quantitative behaviors of such series are well known. In particular, for m = 1 we get one real critical point at r = 0 and a Borel resummable expression, for m = −1 three real critical points and Borel summability is lost, while for m = 0 the three critical points are degenerate and no expansion is possible. Given the relations (A.8) and (A.9), we can easily get the analyticity properties of the coefficient term c p as a function of the coupling g. In particular, we see that c p = c p (g) are analytic at g = 0 for any p and go like (2.6) Other useful examples are given by the 0d reductions of two of the three models considered in this paper, namely the NLSM and the GN models.
The 0d reduction of the non-linear sigma model is essentially the S N −1 sphere. We can define with Ω N = 2π N/2 /Γ(N/2) the volume of the S N −1 sphere. The large N expansion reduces essentially to the Stirling approximation of the Gamma function, which is well-known to be divergent asymptotic. In presence of a chemical potential h, the vacuum energy becomes where γ(a, z) = Γ(a) − Γ(a, z) is the incomplete Gamma function. The behavior of the 1/N expansion is now determined by the expansion of γ(a, z) for large a and z, at fixed ratio z/a. This can be found e.g. in [51], see eq.(8.11.6). After simple algebraic manipulations, we have where Q k (h 2 ) are polynomials of degree k in h 2 for k > 0 and Q 0 = 1. It can be shown that the above series is absolutely convergent for any real h for N > 3. Interestingly enough, while F NLSM (h) and F NLSM (0) are separately non-analytic at N = ∞, their difference F NLSM (h) − F NLSM (0) is a well-defined and analytic function. The 0d reduction of the Gross-Neveu model is given by the following Grassmann integral: where χ = (χ 1 , χ 2 , . . . , χ 2N ) is a set of 2N complex Grassmannian variables. Introducing an Hubbard-Stratonovich like parameter as in (A.1) we get (2.11) The large N expansion of this result is again manifestly divergent asymptotic. In presence of a chemical potential h, the vacuum energy becomes where the last line is readily computed again introducing an Hubbard-Stratonovich like parameter. We finally have Like in the NLSM case, F GN (h) and F GN (0) are separately non-analytic at N = ∞, but their difference F GN (h) − F GN (0) is a well-defined, simple and analytic function. Summarizing, we have found that the 1/N expansion in 0d reductions of large N vector models is asymptotic, and each coefficient in the 1/N expansion is analytic in the t' Hooft coupling. In agreement with what was anticipated in the introduction, the factorial growth of diagrams in perturbation theory is reduced to exponential growth, order by order in 1/N . In contrast, the coefficients in the 1/N expansion we will compute in the 2d models will generally be non-analytic in the t' Hooft coupling because (and only because) of the presence of renormalon singularities. We have also shown that the 1/N expansion of the relative free energy F (h) − F (0) is better behaved than F (0) and is convergent in the 0d reduction of both the NLSM and the GN models. This suggests that the relative free energy can have better convergent properties in 1/N also in the 2d models. It will be explicitly verified that the 1/N expansion of this quantity is indeed convergent in the 2d GN model, while we will not be able to draw firm conclusions on its nature in the NLSM and PCF models.

Free energy and integrability
We summarize in this section the general formulation which applies to the three integrable and asymptotically free models considered in the paper. More details for each model will be spelled out in subsequent sections.
The key observable we will study in this paper is the free energy F (h) as a function of an external field h coupled to a conserved charge. Let H be the Hamiltonian of the theory and Q the charge associated to a global conserved current. The external field h can be regarded as a chemical potential, and we can consider the ensemble defined by the operator H − hQ. (3.1) The corresponding free energy per unit volume is defined by where V is the volume of space and β is the total length of Euclidean time. More precisely, the observable of interest will be the relative free energy From now on, for simplicity, we will refer to F(h) just as the free energy. It was pointed out in [31] that in integrable quantum field theories one can calculate F(h) by using the exact S-matrix and the TBA ansatz, in terms of a linear integral equation. The basic physical intuition behind is the following. Let m be the mass gap of the integrable theory. If the lightest particles in the theory are charged under Q, for h > m the ground state of the theory will no longer be the vacuum, but a state with non-vanishing number density ρ. The latter can be determined in terms of Bethe roots χ(θ) by the TBA equation where θ is the particle rapidity and χ(θ) is supported on the interval (to be determined) [−B, B].
The integral kernel appearing in the Bethe ansatz equation is given by where S(θ) is the S-matrix element of the particles populating the ground state. In all the cases considered in this paper only one species of particles with definite charges are present, so S is a scalar quantity. The energy per unit length e and the density ρ are given by The value of B is fixed by the density ρ and one can eventually obtain an equation of state relating e to ρ. The free energy F(h) is finally obtained by a Legendre transform of e(ρ): In an alternative formulation of the TBA equations, the basic quantity is a function (θ), with support on an interval [−B, B], which describes physically the excitation of holes. This function satisfies the integral equation where now the value of B is determined by the condition (±B) = 0, (3.9) and depends on the external field h. The free energy is then given by We refer the reader to appendix B for the explicit form of the kernel K(θ) in the three models, for a detailed discussion of the existence and uniqueness of the solutions of (3.4) and (3.8), as well as for the analyticity properties in 1/N of K(θ) in each case. Given the free energy F(h), we denote by F k (h) its coefficients in a 1/N expansion (see (3.17)- (3.19) below for the precise definition for each model). The coefficients F k (h) are nonperturbative functions of the external field h and the mass gap m in each model. In order to recast the results in terms of asymptotic expansions of ordinary perturbation theory, we have to define a running coupling constant of some kind. Let us denote by g the coupling constant (before the large N limit) appearing in the Lagrangian description of our models, with beta function given by with β 0 > 0. Since in all the three models considered β 0 ∝ N we can conveniently define a 't Hooft-like coupling as α ≡ 2β 0 g 2 , (3.12) so that the β-function up to two loops reads (3.14) As it is well-known, the first two terms are renormalization scheme-independent, while all the others are not. A useful definition of running coupling is 3 Applying µ∂ µ to (3.15) gives So we see that α is a plausible coupling of the integrable model, in a renormalization scheme where the β-function has exactly the form (3.16). We will refer to this renormalization scheme as the TBA scheme in the following. In the three models we expand the free energy as Several clarifications are in order. In (3.17)-(3.19) the numerical factors have been chosen for convenience and α is the coupling (3.15) evaluated at a convenient scale µ ∼ h that will be spelled out in detail for each model in the next sections. In order to clearly distinguish exact quantities from their asymptotic expansions, we have denoted by F k (h) and Φ k (α, C ± ) the exact and the asymptotic expansion in α of the 1/N coefficients of F(h). The factors Φ k are trans-series of the form where ϕ ( ) k are in general asymptotic divergent series. ϕ (0) k coincides with the perturbative series, while ϕ ( ) k with > 0 are the non-perturbative contributions associated with the trans-series. The latter can present a two-fold ambiguity, encoded in the trans-series parameters C ± . This ambiguity is balanced by the one due to the non-Borel summability (because of the IR renormalons) of the asymptotic series ϕ ( ) k in a way that will be detailed for each model in the next sections. The symbol ∼ stands for "asymptotically equivalent to". 4 Note finally that in the NLSM and in the GN models the Φ k 's do not precisely correspond to the expansions of the F k 's because the relation between α(µ ∼ h) and h given by (3.15) depends on ∆ by means of the factor ξ.

The non-linear sigma model
In this section we present our results for the NLSM. After reviewing general aspects of the theory, we will present two approaches to obtain the exact F(h) up to next-to-leading order, one based on the diagrammatic calculation in the 1/N expansion, and another one based on the expansion of the integral equation from the Bethe ansatz. We will then compare this non-perturbative result with the resummation of the perturbative calculation.

General aspects
The NLSM is described by the Lagrangian density where φ = φ 1 , . . . , φ N is an N -uple of real scalar fields satisfying the constraint In our conventions (3.11) we have The NLSM has a global O(N ) symmetry. The conserved currents are given by where I, J = 1, . . . , N , and we denote by Q IJ the corresponding charges. As in [32,33], we can add a chemical potential h associated to the charge Q 12 . In the NLSM it is convenient to define the 't Hooft coupling α ≡ α(µ = h) , (4.5) where α(µ) is the TBA coupling defined in (3.15). The free energy F(h) was computed in perturbation theory in the coupling constant up to two-loops in [32,33,52]. In terms of the asymptotic expansions defined in (3.17) and (3.21), we have at leading order in 1/N (k = 0) At next-to-leading order (k = 1) the series ϕ (0) 1 (α) has been calculated explicitly and at all orders in [50]. It is obtained by selecting Feynman diagrams with the appropriate power of N , which turn out to be ring diagrams. The resulting series has the factorial growth typical of renormalon behavior, due to integration over momenta, and one finds 5 We can then ask the question of how the 1/N expansion resums this series. We will now present two different ways of computing F 0 (h) and F 1 (h) as exact functions of h.

The 1/N expansion from QFT
The standard way to perform the 1/N expansion is to introduce an auxiliary field σ which implements the constraint (4.2) (see e.g. [53,54]). Once this is done, we obtain the action (4.8) The subscript B denotes bare quantities and the second line is a rewriting in terms of renormalized fields and counterterms. The additional couplings for h = 0 are We do not need a vertex renormalization for h because it couples to a conserved current. Moreover turning on h does not introduce any new UV divergence, so all the counterterms in (4.8) can be taken independent of h. Fixing also the finite part of the counterterms to be h-independent means that we choose the same scheme for h = 0 and h = 0. The only exception is the counterterm C, for which we find that a finite h-dependent shift C h is needed in order to satisfy a certain "renormalization condition" on the observable that we explain below. In order to compute the vacuum energy up to NLO, we need to plug in the action the VEVs including their 1/N corrections where capital letters denote the VEVs and hatted fields denote the fluctuations. Similarly we plug the expansion of the renormalization constants, out of which only Z g and C are non-trivial already at leading order (4.11) It will be convenient to collect some combinations of counterterms δm 2 is the total counterterm for the mass-squared coupling, and δZ 3 the counterterm for the cubic coupling. As we will show, h = 0 induces a non-zero VEV for Φ in the 12 plane, which in turn induces a quadratic mixing betweenσ,φ 1 andφ 2 . In the following we will draw diagrams with the following conventions φ propagator: ; σ propagator for h = 0: ;σ-φ 1 -φ 2 mixed propagator for h = 0: ; σ tadpole vertex: ; VEV insertion: ; Counterterm: .

Leading order
The LO (leading order) vacuum energy density for h = 0 is given by the following diagrams subject to the vanishing tadpole condition (this is equivalent to minimizing the effective potential) (4.14) We did not write explicitly the tadpole condition forφ I but it is easily checked that it is solved by Φ = 0. Evaluating the diagrams in the tadpole condition we find The solution to this equation gives us the physical mass-squared of the scalars at leading order, that we denote as m 2 0 , i.e. Σ = m 2 0 2 . Plugging in the diagrams for the vacuum energy we obtain The LO vacuum energy density for h = 0 is given by subject to the tadpole conditions (now we cannot avoid considering also the tadpole forφ, without loss of generality we assume the VEV to be in the direction 1) Assuming Φ 1 = 0 the tadpole condition forφ has a unique solution for Σ, namely Σ = h 2 2 . Plugging this in the tadpole condition forσ we obtain Using the determination of Z (0) g g 2 in terms of the mass-squared for the theory with h = 0 in eq. (4.15), we get Plugging everything back to the vacuum energy density we obtain

(4.22)
Note that we did not really need the result for (Φ 1 ) 2 here because the dependence on this VEV canceled when plugging the value of Σ.
Taking the difference we get the following LO result for the observable 6 Note that we are slightly abusing notation: the symbols F 0 (h) and F 0 (0) do not denote the same function evaluated at two different arguments, as is clear from the fact that the difference does not vanish for h = 0. This is because we are considering a different stationary point of the effective potential for h = 0, so we are actually taking the difference between the free energies of two different states, whose energies cross when h = m 0 where indeed the observable vanishes (recall that m 0 denotes the physical mass of the bosons at LO). For h > m 0 the vacuum with the condensate Φ 1 = 0 is energetically favored.

Next-to-leading order vacuum diagrams
The NLO (next-to-leading order) vacuum energy density for h = 0 is given by the following diagrams It is a general fact that the various diagrams linear in the NLO correction to the VEVs δΣ, δΦ drop from the NLO energy density thanks to the LO tadpole condition, so we avoided drawing such diagrams. We will postpone the NLO tadpole condition forσ to the subsection 4.2.3 and leave δZ g as undetermined for the time being. The last contribution comes from the NLO correction to the physical mass-squared The leading free energy F0 was computed in [48], see eq.(2.7), where µ there = h here . However, [48] missed the non-perturbative term proportional to m 2 0 , crucial to establish that F0(h = m0) = 0.
when F 0 is re-expressed in terms of the physical mass-squared m 2 . Note that in all the NLO diagrams we can just use m 2 as the mass-squared of the bosons. The first diagram in (4.24) is a closed loop of theσ field. The propagator coming from the resummation of theφ bubbles is where B(m 2 , p) is the "bubble function" Summing up with the other diagrams, we find The NLO vacuum energy density for h = 0 is given by the following diagrams We again used the LO tadpole condition to avoid drawing all diagrams with insertions of δΣ and δΦ. Note that due to the quadratic mixing at this order we receive a contribution from the bubble of theσ-φ 1 -φ 2 propagator (first diagram) and to avoid overcounting we need to subtract the contribution ofφ 1 andφ 2 to the LO closed loop of the bosons (second diagram).
Using a matrix notation, we can express the quadratic action involvingσ,φ 1 andφ 2 as Here p τ is the momentum in the Euclidean time direction. Therefore we have (4.32) We substituted (4.21) inside the determinant. Summing up with the counterterm diagrams, we find The subscript h in the mass-squared counterterm reminds us that the definition (4.12) of δm 2 contains a factor of the VEV Σ, which depends on whether we are at h = 0 or h = 0. Taking the difference F 1 (h)−F 1 (0) we note that all the contributions involving the propagator counterterms cancel and we are left with (4.34) Here we used (4.23) to evaluate the terms involving m 2 0 ∂ ∂m 2 0 . We still have a UV divergent integral and some counterterms that did not cancel in the difference, so we need to fix those to get the finite result for the observable.

Tadpole condition and propagator correction
The NLO Tadpole condition for h = 0 is The two-loop diagram gives Going to the second line we performed the integral in q, which gives Plugging this result and evaluating the one-loop diagrams the NLO tadpole condition can be rewritten as The left-hand side is precisely the combination that appear in the observable in eq. (4.34) multiplied by − h 2 −m 2 2 , so substituting this relation we find (4.39) However, we are still left with a UV divergent integral and some combination of counterterms. To fix this remaining combination we need to consider the correction to the physical mass-squared. The 1/N correction to the inverse propagator of φ for h = 0 is given by Summing up this diagrams, and imposing that the physical mass-squared at NLO is given by eq. (4.25) we obtain Plugging (4.41) in (4.39), the residual dependence on the counterterms and on r cancels and we obtain The integral appearing in this final expression is UV finite, confirming that δC h is a finite shift. In order to fix the finite shift δC h we impose as a further renormalization condition that also at NLO the observable vanishes for h = m. This ensures that after the inclusion of 1/N corrections it remains true that the state with a condensate Φ 1 = 0 becomes favored precisely in the range h ≥ m (recall the comments below eq. (4.23)). By performing the integral for h = m, one can check that this condition is satisfied by fixing δC h = h 2 4π . To write more explicitly the k-dependent part, it is convenient to choose k in the Euclidean time direction and simply plug k = (im, 0), so that We can then drop the imaginary part because it is odd in p τ . Therefore we obtain the final result (4.44)

The 1/N expansion from the Bethe ansatz
As discussed in section 3, since the NLSM is integrable, the free energy F(h) can be calculated exactly by solving the integral equation (3.8). One could think that the 1/N expansion of the free energy can be obtained by simply expanding the kernel in a power series in 1/N . However, this turns out to be subtle, since the kernel is singular in the limit ∆ → 0, see appendix B. The 1/N expansion of the integral equation (3.8) for the NLSM was studied in [48] at the first nontrivial order. The calculation of higher order corrections in similar models introduces additional subtleties, as shown in [49], but in our analysis we will restrict to the next-to-leading order term. It turns out that the 1/N expansion of the kernel requires the introduction of distributions. The expansion of the kernel reads, where δ(θ) is Dirac's delta function. The first two terms in the expansion read (4. 46) In these expressions, and where ψ (1) is the first derivative of the digamma function. The expansion (4.45) has to be understood in the following sense: when acting on a test function f (θ) which is differentiable and vanishes at the boundaries, one has where P means the principal value.
We can now solve the integral equation by using a large N ansatz for the distribution (4.50) and for the endpoint The leading term in the integral equation (3.8) is then given by By splitting the kernel L 1 (θ) as is regular at θ = 0, we can write (4.52) as After integrating by parts and an integration w.r.t. θ, one can also write (4.52) as which is the form used in [48]. We conclude that the large N limit of the distribution appearing in the TBA equation (3.8), 0 (θ), can be obtained as a solution of the singular integral equations (4.52) or (4.56). (4.56) is very similar to the integral equation that one would find for the density of states of a large N matrix integral. The singular part of L 1 (θ) is what one would find for a conventional, Hermitian matrix model. The additional term L r 1 (θ) would correspond to a non-conventional eigenvalue interaction. We note that the integral equation (4.56) does not determine by itself the value of B 0 as a function of h. It was argued in [48] that this value is fixed by requiring the following behavior near the edges of the distribution: When L r 1 (θ) = 0, as it happens in the PCF model studied in the next section, one can show that the condition (4.57) follows from (4.55), by requiring 0 (θ) to be regular at ±B 0 .
The solution to the singular integral equation (4.56) is not known explicitly. However, there is both analytic and numerical evidence that the leading order free energy obtained from this solution, agrees exactly with (4.23). By assuming this analytic value for F 0 (h) one can deduce the value of B 0 [48]: Let us now write down the equation for the next-to-leading correction 1 (θ). By using (4.57), it can be seen that B 0 does not get corrected at that order, and one finds the equation This equation can be solved numerically to obtain the next-to-leading correction to the free energy, given by We checked, for some values of h, that F 1 (h) computed as above agrees with (4.44). The numerical resolution of the singular integral equation is quite time consuming: with few hours of computation we reached an agreement with a relative error of order 10 −7 . Proceeding to higher orders in ∆ in the expansion of the kernel (4.45) and solving singular integral equations to compute F k for k > 1 becomes very challenging. We have been able to extract more coefficients of the series in ∆ by following an indirect procedure: we first solve numerically with high precision the TBA for fixed h 0 and different values of ∆, and from this sequence of F(h 0 ) we compute the value of the asymptotic expansion at each order, using Richardson transforms to accelerate the series. This allowed us to compute, for a given h, the functions F k (h) up to k = 6. These first coefficients do not present a factorial growth as it would be expected from the properties of the kernel discussed in appendix B. However, they are not enough to allow us to make claims on the nature of the 1/N series.

Trans-series expansion and comparison with perturbation theory
We can now compare the exact results for F 0 and F 1 found in the previous subsections with perturbation theory. Up to order ∆, (3.15) gives Given (4.23) and the definitions (3.17) and (3.21), it is straightforward to get The perturbative series ϕ (0) 0 agrees with (4.6), as it should, but in addition we see a nonperturbative single trans-series term, which cannot be captured in perturbation theory. The mismatch between the exact free energy and its perturbative calculation is due to the h-independent term proportional to m 2 . It is given by the non-perturbative free energy at h = 0, evaluated at the non-trivial large N point, which has been calculated in [55] and has the value This suggests that the difference between perturbation theory and the 1/N expansion is only due to the difference between a perturbative and a non-perturbative evaluation at h = 0. We will give further evidence of this when we consider the next-to-leading term in the 1/N expansion. A direct expansion in α of the NLO term F 1 in (4.44) is challenging. Luckily enough, we will verify that it has the same structure of F 0 , namely its trans-series expansion is composed of only two terms ϕ (0) 1 is the perturbative result and its first terms can be read from (4.7). This series is factorially divergent, and its Borel transform has a Borel singularity in the positive real axis (i.e. an IR renormalon), which was analyzed in detail in [50]. We can then use lateral Borel resummations s ± (ϕ (0) 1 )(α), i.e. integrating the Borel transform slightly above or below the real axis avoiding in this way the singularities, indicated respectively by s + and s − . They lead to an imaginary piece. Let us denote by s ± (ϕ (0) 1 ) the function obtained from ϕ (0) 1 by using the lateral Borel resummations of the series (4.7). The results of [50] imply that Im s ± ϕ which is the imaginary contribution of the IR renormalon unveiled in [50]. A detailed numerical comparison indicates that the exact 1/N result (4.44) agrees precisely with the real part of the lateral Borel resummations. Defining where in (4.66) α = 1/ log(h/m). We can use the so-called median resummation, which in this case is simply to write the above result as These results indicate that the trans-series expansion of F 1 is rather trivial, namely we have where C ± = ±i. In fig. 1 we illustrate the agreement between the median resummation and the numerical evaluation of the exact result. In order to show the high precision of the resummation and the agreement with the exact value, in table 1 we report the comparison between the  This is similar to the result found in [46,47] for the PCF model at large N , in which the exact answer is the real part of a Borel-resummed series. One additional insight of the result above is that we have an explicit diagrammatic interpretation of the underlying perturbative series, in terms of the ring diagrams considered in [50].

The principal chiral field model
In this section we present our results for the PCF model. In the first subsection we review some general aspects of the model. Then, we present the large N exact solution for F 0 (h). In the final subsection we "decode" this exact solution in terms of an explicit, resurgent trans-series: we find a trans-series extension of the perturbative expansion, and we show that the exact solution can be obtained from this trans-series by Borel resummation.

General aspects
The PCF model is a quantum field theory for maps Σ : R 2 → SU (N ), with Lagrangian density In our convention (3.11) we have The PCF has a global symmetry SU (N ) L × SU (N ) R , therefore there are conserved charges that can be used to construct the free energy F (h). We will consider a vectorial symmetry SU (N ) V ⊂ SU (N ) L × SU (N ) R , and the corresponding charge will be denoted by Q. Its eigenvalues in the fundamental representation will be denoted by q = (q 1 , · · · , q N ). Two choices of charges have been studied in the literature. In [46,47] one takes This leads to an explicitly solvable large N limit. There is a different setting, considered originally in [34] to determine the mass gap of the theory, in which one chooses the charge , · · · , − 1 2(N − 1) .
We will consider in this paper the setting (5.4). As we will see, the large N limit is also solvable in this case. In addition, the resulting large N free energy leads to an infinite series of IR renormalon contributions which can be calculated explicitly. The relevant kernel for the choice of charges (5.4) is reported in (B.26). In the PCF model it is convenient to define the 't Hooft coupling where α(µ) is the TBA coupling defined in (3.15). The free energy F(h) can be computed in perturbation theory. The one-loop result was presented in [34] for an arbitrary choice of charges, and it is given by where g 2 is the MS coupling defined by It follows from (5.6) that in terms of the 't Hooft coupling α, we have where F 0 is the leading 1/N term of F appearing in (3.18). The free energy F(h) can also be computed from the TBA equations (3.8), (3.10). The TBA solution can be used to extract the perturbative expansion of F 0 (h) up to very high orders in the coupling constant and at finite N , by using the methods in [41,42]. This was done in [43] for the normalized energy density, but the results in that paper can be easily translated into an expansion in α for F 0 (h), and one finds

Exact solution at large N
As in the case of the NLSM, the kernel (3.5) given by (B.26) has a 1/N expression which involves distributions: where the first term is simply The integral equation at leading order in 1/N , (4.56), reads in this case Here, we have denoted B 0 by B, since we will not consider subleading corrections to the value of the endpoint. Due to the simplicity of the leading kernel, (5.12) is the equation for the density of eigenvalues of a Hermitian matrix model with a potential Since the support of 0 (θ) is the full interval [−B, B], we are considering a so-called one-cut solution. This solution can be obtained by using standard matrix model techniques, see e.g. [54]. The density is given by where the function M (θ) can be written as a contour integral around z = 0 To obtain the explicit solution, it turns out to be useful to write V (x) in (5.13) as a power series around x = 0, where, in our case, Then, by using the expansion we find the expression The value of B is determined by the condition (4.57), as in [48], and one finds By expanding the cosh, we find are the moments of 0 (θ). As it is well-known from the matrix model literature, they can be calculated by using the expansion at infinity of the resolvent By doing this, one finds (5.26) Note that, from the point of view of quantum field theory, this is a strong coupling expansion, since small B corresponds to small h/m. Fortunately, this expression can be summed up in closed form in terms of a generalized hypergeometric function, and we obtain in the end From these expressions it is also possible to obtain exact formulae for the large N limits of the density of particles and of the energy density, defined by

Trans-series expansion
We can now address the question of what is the relation between the exact large N result (5.27), and the trans-series expansion of the free energy, including exponentially small corrections. It turns out that all the special functions appearing in (5.27) and (5.21) have simple trans-series expansions which can be used to obtain a trans-series expansion of F 0 (h). Let us start with (5.21). By using the trans-series asymptotics of the Bessel functions, we have the following equality, valid for B > 0: Here, are Gevrey-1 series, s with no subscripts ± denotes the standard Borel resummation, available when there are no singularities in the positive real axis of the Borel plane, and C ± = ∓i. (5.32) As usual in Borel-Écalle resummation, the value of C ± is correlated with the choice of lateral resummation. We also have where Finally, we have the following formula for the generalized hypergeometric function, It follows from (5.27) that F 0 (h) has a trans-series expansion in terms of the small parameters 1/B, e −2B . We are however interested in obtaining the trans-series expansion in terms of the 't Hooft coupling α introduced in (3.15), which makes contact with conventional perturbation theory. To do this, we need the relation between α and B, which is obtained by combining (3.15) and (5.21). This relation can be written as a trans-series equation, and it has a trans-series solution of the form The leading term in this trans-series is given by We can also compute the first exponential corrections. This is better done in terms of b. We find, for the very first terms, (5.40) We obtain in this way the following trans-series structure for F 0 (h): where C 2 ± = −1. To make this completely explicit, we can write (5.41) with the notations introduced in (3.21) as a trans-series in α, e −2/α . We have where we have factorized C ± and dropped the unnecessary subscript 0: The series ϕ ( ) (α) can be read from (5.41), (5.44) Then, we have the following equality: There are many aspects of the above result which are worth commenting in detail, both at the physical and the mathematical level. From the physics point of view, let us note that the trans-series (5.42) has an infinite number of exponentially small corrections, corresponding to an infinite number of IR renormalon singularities. This is in contrast to the NLO term in the 1/N expansion of the NLSM, studied above, and also to the planar solution of the PCF model [46,47] with the choice of charges (5.3). However, all corrections are built up of a finite number of trans-series in the variable B, appearing in (5.27). A similar phenomenon appears in the trans-series solution of certain Riccati ordinary differential equations [56,57], in which all the exponential corrections in the trans-series are obtained from a finite number of building blocks.
From a more formal point of view, one can ask whether the exponentially small corrections are determined by the perturbative series. It turns out that, in this case, Φ(α; C) satisfies the same resurgent equations that the trans-series solution to Painlevé II described in detail in [19]. Namely, we conjecture the following equality of laterally resummed trans-series where S = 2i (5.47) and C is now an arbitrary complex parameter. As shown in [19], (5.46) leads to the following relationships We have explicitly verified some of these equations numerically, by including up to the fourth term in the trans-series. (5.48) gives the so-called Stokes automorphism of Φ ( ) across the positive real axis, and it expresses it as an infinite linear combination of the higher order terms in the trans-series, Φ ( +k) (α). The coefficients in the r.h.s. of (5.48) are called Stokes constants and are explicitly known. Equivalently, (5.48) says that the Borel transform of Φ ( ) (α), Φ ( ) (ζ), has singularities at ζ = 2k, k ∈ Z >0 . By expanding Φ ( ) (ζ) around the k-th singularity, one can obtain Φ (k+ ) (α). In particular, when applied to = 0, (5.48) implies that all the higher order terms in the trans-series can be obtained from the perturbative series. It follows from (5.48) that the Borel resummed trans-series is real for any real C [19]. This is sometimes called a median resummation (see also [58]). Since C ± = ∓i = ∓S/2, (5.45) is a median resummation corresponding to C = 0. For illustration we compare in fig. 2 the exact value of the leading rescaled free energy ( 1 5 ) n = 0 2.68541385336. . . · 10 −9 ±0.0020200401313. . . n = 1 −2.685413850(1) · 10 −9 ∓9.162556(5) · 10 −14 n = 2 −3.9(1) · 10 −18 ±4.581282(5) · 10 −14 n = 3 1.4(1) · 10 −18 ∓3(5) · 10 −20 Table 2: First orders of the truncated trans-series. In the first column the difference between the exact value and its real part and in the second column its imaginary part. The uncertainty correspond to the error associated to the Borel resummation, as explained in the main text.
Φ ≡ s ± (Φ) (C ± ) with the real part of its approximations given by the truncated trans-series The term Φ Other subdominant contributions -such as the one obtained by introducing one or more dummy variables (e.g. a Borel-Le Roy parameter) and minimizing with respect to them-have been neglected (see e.g. section 4 in [59] for an overview of possible numerical recipes to estimate the error when using Borel resummation techniques). We see how the exact result is approached when considering more and more terms in the trans-series expansion. In table 2 we show more in detail the cancellations happening between the different orders of the trans-series. We report, for fixed α = 1 5 , the difference between the exact value and the real part of the truncated trans-series, and its imaginary part. It is interesting to notice how both values approach zero with a change of magnitude happening every two orders, a consequence of the fact that the Φ ( ) 's satisfy the relationships (5.48).
In order to better appreciate how the exact result is approached when taking more and more terms in the trans-series, we show in fig. 3 the resummation of Φ for a fixed α = 1, and as a function of the number of trans-series terms n included in the resummation. In this plot the error interval is given by the left-over ambiguity in the trans-series, namely by the imaginary part of s ± (Φ (n) ) at each order. The numerical error associated to the Borel resummation is sub-leading and has been neglected. Note how quick is the convergence to the exact result and how the choice of the uncertainty gives a reliable estimate of the error. All the statements above were made for the observable F 0 (h), but a similar trans-series expansion can be made for the normalized energy density using (5.29). In this expansion it is convenient to introduce the coupling a defined by where the numerical factor inside the log has been chosen to match this definition with that in [43]. Note that the coupling a here was denoted by α there. The resulting trans-series is The first line reproduces the result obtained in [43]. The trans-series appearing here, in terms of a, has the same formal properties of its close cousin (5.41), like for example (5.46).
Our main conclusion is that, in this example, the exact large N free energy of the PCF model can be obtained by a median Borel resummation of a non-trivial resurgent trans-series, and therefore provides a beautiful success for the program of resurgence in an asymptotically free quantum field theory. Mathematically, this works because the building blocks of the exact solution (5.27) are special functions with known trans-series representations. From the physical point of view it is however gratifying to have a non-trivial example of resurgence at work with infinitely many non-trivial IR renormalons, yet analytically tractable.
In a recent tour de force, Abbott and collaborators [25,26] were able to obtain detailed information on the trans-series expansion for the normalized energy density in the O(4) sigma model, which is nothing but the PCF model we are studying at N = 2. By extrapolating numerical results to analytic results, they obtained an expression very similar to (5.53), and they gave evidence that the exact answer can be recovered by median resummation of the transseries. Our results are an analytic version, at large N , of their result for N = 2.

The Gross-Neveu model
In this section we present our results for the GN model. In subsection 6.1 we review some general aspects of the model. Then, in subsection 6.2 we present the trans-series expansion for F 0 (h) and F 1 (h) and show that these trans-series cannot be reconstructed using resurgence. Finally in subsection 6.3 we numerically study the 1/N series expansion of F up to high order, we analytically continue the series beyond its radius of convergence, and see how this continuation matches with known dualities between GN models at low N and sine-Gordon theories.

General aspects
The Lagrangian density describing the Gross-Neveu (GN) model [5] is where χ = (χ 1 , . . . , χ N ) is a set of N Majorana fermions in two dimensions. As is well-known, a non-perturbatively generated mass gap and spontaneous breaking of a Z 2 chiral symmetry occur in this theory. For N > 4 the lightest particle in the spectrum is the fundamental fermion in the Lagrangian (6.1). In our conventions (3.11) we have 2) The GN model has a O(N ) global symmetry, under which the N fermions transform as vectors, with conserved currents given by J IJ µ =χ I γ µ χ J . where α(µ) is the TBA coupling defined in (3.15). The first two terms F 0 and F 1 in the 1/N expansion (3.19) have been analytically computed in [35,36] using both QFT and TBA techniques. They read and Shi(x) is the hyperbolic sin integral function (6.7)

Trans-series expansion
Let us now work out the explicit form of the first terms of the trans-series Φ k (α) defined in (3.19), for k = 0, 1, using (6.5) and (3.15) relating α and h. This map depends on ∆ through the ξ factor in (6.2). Up to order ∆ we get Note that the perturbative expansion ϕ (0) 0 is trivial and all trans-series terms ϕ ( ) 0 with ≥ 1 are truncated, yet non-vanishing. We see that here is no way to reconstruct the non-perturbative terms ϕ ( ) 0 (α) with ≥ 1 from ϕ (0) 0 . Since the leading order is somewhat trivial, it is useful to go through the next-to-leading order Φ 1 , where each term ϕ ( ) 1 (α, C ± ) has a non-trivial asymptotic expansion in α. After some algebra, the first terms read This series is non-Borel resummable, but can be studied analytically. Its Borel transform is The simple pole at t = 2 hinders Borel summability. We can deform the contour to avoid the pole, passing either above (C + ) or below (C − ) it. The Borel resummation of this series gives then Nicely enough, the ambiguity in the imaginary part in (6.13) is exactly the same appearing in ϕ (1) 1 (α, C ± ). The two contributions cancel each other if we choose the contour C ± for C ± , respectively. The terms in parenthesis in the series (6.10) for ϕ (1) 1 (α, C ± ) form the same asymptotic series (6.11). Their Borel resummation gives where for simplicity we have not reported the first four terms of ϕ 1 appearing in (6.10) (that's why the ⊃ sign instead of the equality sign). Again, if we pick up the contour C ± , the imaginary part in (6.14) cancels respectively the imaginary terms proportional to −C ± appearing in the second trans-series ϕ (2) 1 . In the spirit of resurgence, imaginary parts nicely match between one series and the next, but we see a plethora of real non-perturbative terms which cannot be detected. The asymptotic series ϕ (1) 1 and those with ≥ 2 cannot be reconstructed from the knowledge of ϕ (0) 1 only. In order to quantify and illustrate the phenomenon, in fig. 4 we compare the exact free energies F 0,1 , rescaled by h 2 , with the median Borel resummation of their perturbative series expressed in terms of h and m, i.e s med (ϕ (0) 0,1 )(1/ log(2h/m)). 7 We restrict to the perturbative series because higher trans-series terms can not be obtained from a resurgence analysis. As expected, for both k = 0, 1, the perturbative and full results are in good agreement at large h, i.e. weak coupling, but they significantly differ at strong coupling, when the terms with ≥ 1 are no longer negligible.
All the above analysis could be repeated for the energy density e(ρ), Legendre transform of F(h), expanded in ∆ and written as a trans-series as The asymptotic expansion in this case is more conveniently written in terms of the coupling a ≡ α(µ = 2πρ) .  while for ϕ 1 we have We see that the expansions of Φ 0 and Φ 1 in terms of a are very similar to those of Φ 0 and Φ 1 in terms of α. The perturbative series ϕ (0) 1 in (6.19) agrees with the perturbative expansion found in eq.(A.12) of [43] using the techniques of [41,42], while ϕ ( ) 1 , with > 0, are non-perturbative terms that could not be captured in that analysis. All the considerations made above about imaginary part cancellations and impossibility of recovering the non-perturbative terms from the perturbative series of F(h) apply also for e(ρ) and will not be repeated.

Higher orders in the 1/N expansion
In contrast to the NLSM and PCF models, the kernel in the GN model is analytic at N = ∞. This implies that the TBA solution can easily be expanded in powers of 1/N , with each term a regular function of θ, making it possible to solve the integral equations (3.8) in a systematic 1/N expansion: The kernel coefficients K k (θ) are trivially derived from the GN kernel reported in (B.32). The B k with k ≥ 1 can be expressed in terms of the values of the m (B 0 ) with m ≤ k and their derivatives by solving recursively the condition (3.9) at each order in ∆. For instance, for the first few orders we have where we used the fact that K n (θ) and n (θ) are even functions. Solving the equation at each order in ∆, we can compute iteratively all the k (θ) knowing the values of m (θ) for m < k and the B q with q < k − 1. Finally, in order to compute F k (h) it is enough to expand (3.10) in ∆: where we used once again the fact that n (θ) is an even function.
In order to make more explicit the procedure and show how the iterative process starts, let's rederive F 0 (h) and F 1 (h) given in (6.5). At order ∆ 0 we simply have The leading order free energy reads as already reported in the first equation of (6.5). At order k = 1 we have and Performing the integral we get and from it (6.29) where Chi(x) is the hyperbolic cos integral function Given that 0 (B 0 ) = 0, the next-to-leading term of the the free energy reads Plugging (6.28) in (6.31) and computing the integral gives the second equation of (6.5). At higher order in ∆ the computation becomes analytically prohibitive. On the other hand, it is straightforward to proceed numerically and, for a given B 0 , automatize the iteration procedure to compute higher order terms in ∆. We have been able to compute with high precision, for different values of h, F k up to k = 28. This allowed us to study the large order behavior of the series. We get 8 F k ∝ ρ −k sin(k ϑ) , (6.32) where ρ and θ are two parameters that we can numerically evaluate. For example, for h = 3 we get ρ = 0.50 ± 0.02 and ϑ = 0.35 ± 0.07 . It is now natural to ask if the analytic continuation of the series in ∆ beyond |∆| ≥ 1/2 contains any physical information. This analytic continuation can be obtained by considering Padé approximants PF [m/n] (∆, h) of the series of F k we computed. It is known that for convergent series, parametrically diagonal Padé approximants converge (in capacity) to the exact function and the location of their poles and zeros define an appropriate locus of branch-cuts connecting branch-point singularities [61] (see e.g. app. D of [62] for a brief overview and [63] for a comprehensive introduction). Moreover, as we will see below, F(∆, h), at fixed h, is analytic at ∆ = ∞ and non-vanishing, hence diagonal approximants are the optimal choice to reconstruct the function.
We calculated PF [14/14] (∆, h) for a given set of values of h, and compared the result with F(∆, h) computed by directly solving (numerically) (3.8) at given h and ∆. We find full agreement for all values of h sampled and for 0 < ∆ < 1 2 , and consider it a sanity check of the correctness of the coefficients F k . The location of the poles and zeros of PF [14/14] (∆, h) shows that the point ∆ = 1/2 is a branch-point of F(∆). On the other hand, no singularities appear for ∆ < 0 (as expected from the form of the GN kernel) and hence we can reliably continue F(∆) for ∆ < 0 using its approximant PF [14/14] (∆, h). The two interesting points to discuss are ∆ = 1/2 (N = 4) and ∆ = −∞ (N = 2).
For N = 4 the stable particles in the model are the kinks and their mass equals Since the fermions are exactly at threshold and are marginally unstable, we can use (3.8) to compute the free energy by choosing either kinks or fermions as particles populating the vacuum, but some care is needed. It is useful to briefly review how the analysis goes [35].
Kinks are in the (±1/2, ±1/2) spinorial representation of O(4), so their charges are half those of the fermions. When h/2 > m/2 (or h > m) the vacuum is populated by kinks with O(4) components (1/2, 1/2) and (1/2, −1/2), which do not interact with each other. The Smatrix is identical for the two chiralities and the associated kernel is If we consider kinks, the associated TBA equation is where k describes the excitation of the kink holes. Given k (θ), the free energy is computed as where the factor 2 counts the two-fold degeneracy of the kink and exactly compensates for the 1/2 factor in the mass. The kink kernel (6.35) is naively 1/2 of the fermion kernel in (B.33) for ∆ = 1/2. If we take the limit carefully, however, we also get a δ function because lim y→0 1 π y y 2 + x 2 = δ(x) . (6.38) Hence lim where we denote by K f the fermion kernel. So, for ∆ → 1/2, the TBA equation (3.8) for the fermion excitation holes f becomes which is identical to (6.36) with Note that we would not get the correct result by setting N = 4 directly in the fermion case, because in this way we would not detect the δ(θ) term in (6.39). On the other hand, the analytic continuation of PF [14/14] (∆, h), computed using fermion states for ∆ < 1/2, gives the correct result at ∆ = 1/2. The N = 4 model is also equivalent to a pair of decoupled sine-Gordon models: where b is the inverse radius of the compact scalars. For b 2 > 4π the only asymptotic states in the sine-Gordon models are given by kinks and anti-kinks. In presence of a chemical potential h, the vacuum gets populated by kinks (and no anti-kinks) with a kernel whose Fourier transform is given by [64] K SG (ω, p) = sinh π(p+1)ω where the parameter p is defined in terms of b as follows: On the other hand, the Fourier transform of the kink kernel (6.35) reads (6.45) Interestingly enough, lim p→∞ K SG (ω, p) = K k (ω) , (6.46) so there is an equivalence of the free energy F (h) in the two models provided b 2 = 8π (6.47) in the two sine-Gordon models. Note that when b 2 → 8π, the coupling of the sine-Gordon interaction at the same time becomes marginal and vanishes. From (6.42) it might seem that we get in this limit a pair of decoupled free field scalars, but in fact this is an artefact. The only asymptotic states are kinks and these are still interacting, as evident from the non-triviality of the kernel (6.46). 9 In the correspondence the kink mass of the GN model is mapped to the mass of the sine-Gordon kink: m k = m N=4 SG . For N = 2, i.e. |∆| = ∞, the TBA equations (3.8) are physically meaningless, since fundamental particles are no longer asymptotic states. Yet, we can wonder if its analytic continuation is of physical interest and is related in some way to the actual N = 2 Gross-Neveu model. The latter is nothing else than the Thirring model, famously dual to the sine-Gordon theory [66]. If we approach the infinite limit from negative values of ∆, we see from direct inspection of either (B.32) or (B.33) that the GN kernel is analytic at infinity (recall that the digamma function ψ(z) is meromorphic with simple poles at z = −n, with integer n ≥ 0). By direct inspection we can also check that this kernel is a contraction: All iterated kernels are hence bounded and the solution of the TBA equation is analytic in ∆. By taking the analytic continuation of (B.34) for ∆ < 0 and then the limit ∆ → −∞ we immediately see that the Fourier transform of the kernel equals the kernel for the kink in the sine-Gordon theory as p → ∞ (and the kernel for kinks in the N = 4 GN model). Quite interestingly, the in principle meaningless |∆| = ∞ point of the GN TBA equation (3.8) is related to the free energy F(h) of a sine-Gordon model at b 2 = 8π! In the correspondence the fermion mass appearing in (3.8) is identified with the mass of the sine-Gordon kink: m f = m N=2 SG . We can match the values of h in the N = 2 and N = 4 theories, since h multiplies a conserved quantity and does not renormalize. In particular, we can compute the ratio h N=2 /h N=4 . Recall that h is the chemical potential for a single U (1) current of the formχ 1 γ µ χ 2 . Bosonizing we haveχ -0.7028. . . -1.4056352. . . PF [14/14] (∆, h 0 = 3) -0.7029(2) -1.4056353(10) Table 3: Comparison between the numerical values of the free energy at h 0 = 3 obtained by a direct numerical evaluation of the TBA equation and by using the analytic continuation given by PF [14/14] (∆, h 0 ) for ∆ = 1/2 and ∆ = −∞.
When N = 2 the scalar is free and its Lagrangian reads 50) while for N = 4 the Lagrangian of the two scalars at b 2 = 8π reads 10 j=1,2 Rescaling the fields to have a canonically normalized kinetic term in (6.49) we get When h < m the vacuum is empty and F(h = m) = 0 for both N = 2 and N = 4. This implies that the kink mass ratio is the same as in (6.53), m N =2 SG = 2m N =4 SG and hence that the free energies of the two models are related, namely F(∆ = −∞, h) = 2F(∆ = 1/2, h) . (6.54) In figure 5 we plot PF [14/14] (∆, h 0 ) as a function of ∆ at fixed h 0 . At ∆ = 1/2 the black dot corresponds to the (correct) free energy numerically computed using (6.37), while the red dot is the (wrong) value one would get by naively setting ∆ = 1/2 in the fermion kernel appearing in (3.8). We see that the analytic continuation given by PF [14/14] (∆, h 0 ) gives the correct value. The dashed black line corresponds to the asymptotic value for ∆ = −∞ which should equal to 2F(∆ = 1/2, h), according to (6.54). In table 3 we compare the numerical values of F(h 0 ) obtained by a direct numerical evaluation of the TBA equation and by using the analytic continuation given by PF [14/14] (∆, h 0 ) for ∆ = 1/2 and ∆ = −∞. The two results are in total agreement. 10 Note that the relation between the inverse radius of the sine-Gordon and the GN coupling is different in the two cases, namely (6.51) Both for N = 2 and for N = 4 the kernel on the GN side (for fermions and for kinks, respectively) coincides with the kernel of the sine-Gordon kink, when their corresponding b 2 parameters are set to 8π. Curiously, this corresponds to g 2 = 0 for N = 4 and g 2 = −π 2 for N = 2.

Conclusions
The 1/N expansion provides a non-perturbative resummation of the conventional perturbation theory. In this paper we have discussed the interplay between resurgence and the 1/N expansion in three integrable theories with a continuous global symmetry: the non-linear sigma model, the principal chiral field and the Gross-Neveu models. All these theories have marginal interactions, they are UV free and are affected by IR renormalon singularities. A notable observable is the relative free energy F(h) defined in (3.3). Its special role comes from the fact that it is possible to compute it exactly using TBA techniques [31] and has a non-trivial structure (unlike, e.g., S-matrix elements in integrable theories). Standard large N QFT and/or TBA techniques also allow to analytically determine the first 1/N coefficients F k (h) defined in (3.17)- (3.19).
We have computed F 0 and F 1 in the NLSM, given respectively in (4.23) and (4.44), and determined F 0 analytically in the PCF model (for the choice of charges in [34]), given in (5.27). Crucial for the latter computation has been the observation that the NLSM and the PCF kernels can be expanded in 1/N if the non-analytic term is treated separately. In this way we have also been able to check F 1 in the NLSM by using TBA techniques. These expressions, as well as the previously known coefficients F 0,1 in the GN model [35,36], have then been compared to the asymptotic series expansion, one gets in terms of the coupling constant defined in (3.15). While the perturbative asymptotic expansion agrees with the ones previously determined [43], we get a plethora of non-perturbative trans-series terms which are associated to the non-Borel summability of the series due to the presence of IR renormalons.
The final results turned out to be different in the three models. In the NLSM F 0 contains a non-perturbative term which cannot be captured from the (trivial) perturbative expansion. On the other hand, the median Borel resummation of the perturbative series reconstructs the full next-to-leading coefficient F 1 , see fig.1. The series for F 1 is non-Borel resummable because of the presence of an IR renormalon, yet somehow unexpectedly no non-perturbative terms are missed. In the PCF model the expansion of F 0 gives rise to the non-trivial trans-series (5.41), with a perturbative series affected by an infinite number of IR renormalon singularities. In this case, resurgence techniques work nicely and allow us to reconstruct the full answer from the perturbative series. In the GN model the expansion of both F 0 and F 1 give rise to transseries (6.9) and (6.10) which can not be reconstructed from the perturbative series only, using resurgence.
We also studied the behavior of the 1/N series for F(h). The non-analyticity of the kernel at N = ∞ for the NLSM and the PCF models suggest that the 1/N expansion of (θ) and χ(θ) should be divergent asymptotic. This points towards a divergent 1/N expansion of F(h), as well as of its Legendre transform e(ρ). In contrast, the 1/N expansion of F(h) (and e(ρ)) in the GN model is expected to be convergent. We have numerically computed higher values of F k in each model (for some values of h) in order to verify these expectations. Our results for the NLSM and PCF models are inconclusive. The number of coefficients F k we computed does not allow us to establish whether the series are convergent or divergent asymptotic. The first possibility is not in contradiction with the 1/N non-analyticities of the kernel, because F(h) is obtained by integrating (θ) over rapidities and we cannot exclude that these non-analyticities are smoothed out by the integration procedure. It would be nice to settle this issue in future studies. In the Gross-Neveu model the expected convergence of the 1/N series is numerically confirmed. We analytically continued the series beyond its radius of convergence (see fig.5), where the TBA equations (3.8) and (3.4) no longer make sense, and showed how this continuation gives values of F(h) in complete agreement with those obtained for the sine-Gordon theories dual to the GN models with N = 2 and N = 4.
There are several directions worth exploring in future studies. From the point of view of the general theory of resurgence, our most important finding is its breakdown in certain models, when combined with the 1/N expansion. By a breakdown of resurgence we mean that the structure of non-perturbative corrections at each order in the 1/N expansion can not be predicted from the study of the perturbative series only. A better understanding of this breakdown is perhaps the most important problem open by our investigation. There are two possibilities here. One possibility is that this is a feature of the 1/N expansion and does not apply at finite N . It might happen that, when fixing the order in the 1/N expansion, the resulting perturbative series in the coupling constant is not sufficiently generic and cannot be used to predict non-perturbative corrections. This would be somewhat similar to the "Cheshire cat resurgence" in supersymmetric theories [67]. The other possibility is that the phenomenon we have found is generic in theories with renormalons. The detailed study performed for the O(4) sigma model in [25,26] validates standard resurgence expectations and seems to favor the first possibility. Clearly, additional studies are necessary in order to clarify this fundamental issue. A detailed resurgent analysis of the free energy of the GN model at finite N , along the lines of [25,26], would be very useful. It would be also important to understand why resurgence seems to be so successful in the PCF model, at least at leading order in 1/N , but fails in the GN model, though both models belong to the same universality class of integrable, gapped, and UV-free theories. In particular, it would be interesting to clarify if this is related to the different analyticities properties in 1/N of the kernels in the two theories. More generally, it would be important to study observables other than F(h) which can be computed exactly in the 1/N expansion and can be analytically decoded in terms of trans-series.
Perhaps it is useful to distinguish three different levels of validity of the theory of resurgence, in order to understand what is at stake. The first, more general level of validity is the statement that observables in quantum theory are given by ambiguity-free Borel-Écalle resummations of trans-series. This statement is probably true and it is implicit in many of the early studies of renormalons, like e.g. [28,29]. All of our results in this paper, including the example of the GN model, vindicate this first level of validity. The second level of validity is the stronger statement that the trans-series is fully determined by its perturbative part, up to the numerical values of the trans-series parameters. It is this second level of validity that breaks down in some of the examples that we have studied, when restricting to a fixed order in the 1/N expansion. Finally, a third and largely independent issue is whether renormalon singularities and the associated trans-series have a semi-classical interpretation, in terms of expansions around saddle-points of a classical action. It has been argued that, after a twisted compactification, the renormalon sectors of the PCF model and the NLSM can be interpreted semiclassically [68,69]. Note however that the successful examples of resurgence that we have considered (like e.g. the PCF model at large N ) are independent of such a semiclassical interpretation. It is perfectly conceivable (and, in our opinion, quite likely) that, for theories with renormalons at infinite volume, one does not have a semiclassical interpretation of the trans-series, but some version of resurgence will be still valid. In that scenario, a crucial open question will be to find a generalization of perturbation theory which makes it possible to calculate the trans-series from first principles, and without relying on resurgence properties or integrability. The use of the OPE, as in QCD sum rules, goes along this direction, but it is clear that a more general procedure has to be devised in order to compute general observables for which OPE techniques are in principle not applicable, as it is the case of the free energy studied in this paper.
The results that we have obtained apply to the TBA renormalization scheme defined in (3.15) and might not be valid in others. It would be useful to better understand if and to what extent resurgence methods depend on the renormalization scheme (e.g., do they apply in MS?), given also the impact that the choice of scheme can have when resumming perturbative series even in absence of renormalons, see e.g. [70,71].
Finally, it would be very interesting to extend these considerations to non-integrable theories. One possibility would be to compute F(h) in the quartic linear O(N ) model when m 2 < 0 at some order in 1/N and check if the result can be reconstructed from a perturbative expansion around the naive vacuum, where IR renormalons have been shown to appear [72].
like parameter σ to rewrite the quartic term (x · x) 2 as Inserting in (2.1) and integrating over x gives and we rescaled σ → √ N σ. The function K is complex and the original contour of integration is not a downward flow, so that we have to decompose it in terms of Lefschetz thimbles. We get two critical points The function K has a branch-cut singularity at For p 1 we then get (2.4), which agrees with the earlier result (2.16) of [74] for m = 1. Note that the large order behavior for m = 1 and m = −1 are related, since θ → π − θ andÎ c → 1/Î c when m → −m. We see that the large order coefficients oscillate with a period given by θ.
When θ = π/2, c p = 0 for even p and one has to look at the sub-leading O(1/N ) corrections. This situation is realized when m = 0. In this case, the integral (2.1) simply equals to , (A. 10) and (2.4) simplifies to The sub-leading contributions are captured by the full resurgent (asymptotic) formula [73] whereĉ η q are the first terms of the expansion around the adjacent saddleẑ c , normalized so that c η 0 = 1. Interestingly, c η q = c q , because these saddles are all equivalent for m = 0. We easily get This series alternates every two terms and has odd terms parametrically larger than the even ones by a factor of p, for large p.
where L(θ) is a given continuous function, K(θ) is an integral kernel, and f (θ) is the function to be determined. Equations of this kind are known as non-homogeneous Fredholm linear integral equations. In this appendix we would like to show the existence and uniqueness of the solutions of (3.4) and (3.8) in the three models considered in the paper, as well as some analyticity properties in 1/N of the kernel K(θ). Let us start by briefly reviewing basic mathematical facts. The above equation (B.1) can be compactly written as a fixed point equation Importantly, the operator T is a contraction if Indeed, for arbitrary functions f 1 and f 2 we have If T is a contraction, by the contraction theorem (also known as Banach fixed point) the solution T f = f exists and is unique, and it can be written schematically as The key quantity to study is the kernel K, which will be discussed separately for each model in the next subsections. Our goal will be two-fold: we will first show that K is a contraction, hence proving the existence and uniqueness of the solution, in each case. Secondly, we will study the analyticity in 1/N of the solution. Thanks to (B.5), the analyticity of the solution in a certain region can be established by proving that all the iterated kernels K p are analytic in that region (in particular it is not enough to establish this property for K alone). We will see that there is analyticity in a disk close to the origin of the 1/N plane, but not including it, for the NLSM and PCF, while in the case of GN the solution is analytic in the origin within a certain radius. The analysis will be quite detailed for the NLSM and more sketchy for the PCF and GN models.

B.2 Principal Chiral Field
In the PCF model the kernel equals We can similarly study the analyticity in 1/N =∆ of K. We will be very brief since the analysis is similar to the one performed in the NLSM. Let∆ = z/2 be complex and focus on a small disc D around the origin, defined as in the NLSM case. The function F (θ, z) is analytic around z = 0 and can be expanded for small z as |F (θ, z)| = |z 2 F (θ) + O(z 3 )| < (λ + δ) 2 |F (θ)| . (B.31) For sufficiently small z, the analyticity properties of the PCF kernel coincide with those of both the NLSM model and the Lieb-Liniger model. The iterated kernels K (p+1) are bounded as in (B.18), with the replacement (λ+δ)|F (δθ)| → (λ+δ) 2 |F (δθ)| in the last term of the second row of (B.18). The integral involving |F (δθ)| is easily bounded by a finite constant, so we conclude that the kernel in the PCF model is analytic for 0 < |∆| < 1, but not at∆ = 0 and at∆ = ±1, as evident from (B.28).

B.3 Gross-Neveu model
In the Gross-Neveu model the kernel equals where ∆ = 1/(N − 2), which can also be written as (−1) n n − 2∆ (2∆ − n) 2 + (θ/π) 2 − n n 2 + (θ/π) 2 . (B.33) Its Fourier transform for ∆ < 1/2 reads where M is finite. So, for sufficiently small ∆, K is a contraction. Numerically we see that K is a contraction for arbitrary ∆, not necessarily small. So the unique solution (B.5) is guaranteed to exist. The analyticity in ∆ is trivial in the GN model (see also Appendix A of [35]). Let again be ∆ = z/2 and z = λ + α + iβ, where α and β span a disc of radius δ around the point λ, with λ > 0. Suppose that for a certain p ≥ 1 The relation (B.36), if valid for p, is then also valid for p + 1. Since it applies for p = 1, it follows that it is valid for any p ≥ 1. We see that K and the associated solution f (θ) are analytic for small ∆ > 0, including ∆ = 0. By looking at (B.33) we can determine the radius of convergence of the small ∆ expansion. The above expression is analytic up to the point 2∆ = z < 1. Replacing z = 1 + w in (B.33) gives π 2 K(θ) = w w 2 + (θ/π) 2 − F (θ, w) .
(B. 38) with the function F as in (B.11). We see that around z = 1 the kernel is locally identically to the one of the NLSM close to the origin. Hence the point z = 1 is non-analytic. We conclude that the large N expansion should be convergent with a radius of convergence around ∆ = 0 ρ = 1 2 . (B.39)