Renormalons in integrable field theories

In integrable field theories in two dimensions, the Bethe ansatz can be used to compute exactly the ground state energy in the presence of an external field coupled to a conserved charge. We generalize previous results by Volin and we extract analytic results for the perturbative expansion of this observable, up to very high order, in various asymptotically free theories: the non-linear sigma model and its supersymmetric extension, the Gross-Neveu model, and the principal chiral field. We study the large order behavior of these perturbative series and we give strong evidence that, as expected, it is controlled by renormalons. Our analysis is sensitive to the next-to-leading correction to the asymptotics, which involves the first two coefficients of the beta function. We also show that, in the supersymmetric non-linear sigma model, there is no contribution from the first IR renormalon, in agreement with general arguments.


Introduction
Understanding the large order behavior of perturbative series in quantum theory is a possible route to unveiling non-perturbative effects. In quantum mechanics and in many super-renormalizable field theories, the coefficients of the perturbative series grow factorially, due to the growth of the total number of diagrams with the number of loops [1]. This behavior is controlled by instantons [2,3], and therefore it has a semiclassical description. The relation between instantons and large-order behavior has led to many beautiful results and it has evolved into the theory of resurgence, which provides a universal structure linking perturbative and non-perturbative sectors in quantum theories (see e.g. [4,5] for a presentation of instanton-induced large order behavior, and [6,7] for reviews of the theory of resurgence).
However, in the 1970s it was found that, in renormalizable field theories, the large order behavior of the perturbative series involves a different type of phenomenon [8][9][10][11][12]: one can find specific diagrams which grow factorially with the loop order after integration over the momenta. 1 These diagrams are usually called renormalon diagrams (see [17] for an extensive review and [18] for an invitation to the subject). They lead to singularities in the Borel plane of the coupling constant which, following [17], we will call renormalon singularities, or renormalons for short. Depending on the region in momenta which leads to JHEP04(2020)160 the factorial growth, one has UV or IR renormalons. The analysis of renormalon diagrams is mostly based on heuristic and plausibility arguments. It is relatively easy to find sequences of diagrams with the appropriate factorial growth, like the famous bubble chain diagrams in QED and QCD. However, the resulting behavior could be in principle corrected or even cancelled by some other set of diagrams. Usually, the identification of renormalon diagrams is combined with some type of large N limit (here N can be the number of components of a field in the theory, or the number of fermions, or any other convenient counting parameter), so that one can at least argue that renormalon effects appear unequivocally at large N . In cases where the operator product expansion (OPE) is available, IR renormalon effects can be shown to correspond to non-perturbative condensates in the OPE [11], and this is usually regarded as evidence for both IR renormalon physics and the validity of the OPE.
Renormalon effects are believed to control the large order behavior of perturbative series in many renormalizable theories, as they are typically more important than instanton effects. Therefore, the cleanest way of establishing the presence of renormalon effects is to show explicitly that the perturbative series has the asymptotic behavior dictated by them. However, it is in general difficult to produce explicit values for a large number of coefficients in the series, so this type of tests are difficult to make. A notable exception is the tour de force numerical computation in [19,20], which gives a beautiful and precise test of renormalon predictions in Yang-Mills theory.
Given the subtleties of renormalon physics, it is useful to look at simple field theories where one has more analytic control. For example, renormalons in the non-linear sigma model at large N were analyzed in some detail in [21][22][23][24][25], and very recently evidence for the dominance of the leading IR renormalon was obtained numerically in [26], in the case of the principal chiral field (PCF) [27]. Many asymptotically free theories in two dimensions turn out to be integrable, i.e. the S-matrix is known exactly. The Bethe ansatz can be then used to compute the free energy of these theories in the presence of an external field coupled to a conserved current [28]. This has made it possible to obtain the exact mass gap of these theories in various cases [29][30][31][32][33][34][35][36] (see [37] for a review). It was shown by Volin in [38,39] that one can extend the mass gap calculation and extract from the Bethe ansatz the full perturbative series for the vacuum energy, as a function of the running coupling constant. Volin worked out the example of the non-linear sigma model, where he addressed some aspects of the large order behavior of the perturbative series. Related work appeared before in [40,41], where the PCF was analyzed, but with a different choice of conserved charge than what is made in [31]. The leading large N contribution to the ground state energy was obtained at all orders in the coupling constant, and it was noted that the resulting factorial divergence is due to renormalon effects.
In this paper we generalize the results of [38,39] in two directions. First of all, we streamline the method of resolution of the integral equation by combining it with the Wiener-Hopf method, as we already did in [13]. This makes it possible to obtain the perturbative series for the vacuum energy up to very large order in various integrable models, namely, the supersymmetric O(N ) non-linear sigma model [42], the SU(N ) PCF [27], and the O(N ) Gross-Neveu (GN) model [8], in all cases for arbitrary finite N . Second, we do a precision analysis of the resulting perturbative series in order to test the predictions JHEP04(2020)160 of renormalon physics. Care is needed since the large order behavior mixes IR and UV renormalons, and one needs to disentangle their contributions to the asymptotics. As a result, we are able to test the predictions of renormalon physics to next-to-leading order in the asymptotics. This subleading correction involves the two coefficients of the beta function of the theory. Interestingly, we find that, in the supersymmetric non-linear sigma model, the contribution from the would-be leading IR renormalon is absent. This is in line with the results in [43,44], where it has been shown that the occurrence of IR renormalon singularities is restricted by supersymmetry.
The organization of this paper is as follows. In section 2 we review the theories that we will analyze, the observables that we want to compute, and we explain how to extract perturbative series from the Bethe ansatz. In section 3 we review the predictions or renormalon physics for the large order behaviour of these perturbative series. Then, we compare these expectations to our data. Finally, we conclude with some open problems raised by this investigation.

Integrable asymptotically free theories in two dimensions
In this paper we will consider integrable quantum field theories in two dimensions which are also asymptotically free. Our convention for the beta function is and we will denote With the convention above, asymptotically free theories have β 0 > 0. We will consider two types of theories: the "bosonic" theories include the non-linear O(N ) sigma model, its N = 1 supersymmetrix extension, and the SU(N ) PCF. The "fermionic" theory will be the O(N ) Gross-Neveu model. Let H the Hamiltonian of any of these theories, Q the charge associated to a global conserved current, and h an external field coupled to Q. The external field h can be regarded as a chemical potential, and as usual in statistical mechanics we can consider the ensemble defined by the operator 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.
When h is large, and since the theories we are considering are asymptotically free, one can calculate F (h) in perturbation theory. We can use the renormalization group (RG) JHEP04(2020)160 to re-express the perturbative series in terms of the RG-invariant coupling g 2 (µ/h, g), defined by g can be expressed in terms of Λ/h, where Λ is the dynamically generated scale, which we define as At leading order we have, The theories we will consider are also integrable, and their S matrix is known exactly. It was shown in [28] that this makes it possible to calculate the free energy (2.4) by using the Bethe ansatz. After turning on the chemical potential h beyond an appropriate threshold, there will be a density ρ of particles charged under the conserved charge Q, 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 where S(θ) is the S-matrix appropriate for the scattering of the charged particles. The energy per unit volume and the density are then given by Finally, the free energy F (h) can be obtained as a Legendre transform of e(ρ): Note that the first equation defines ρ as a function of h. Integrable asymptotically free theories in two dimensions have been a useful laboratory to test general expectations from QFT. For example, in asymptotically free massless theories, the masses of the particles in the spectrum are expected to be proportional to the dynamically generated scale Λ, but the calculation of the proportionality constant is a difficult non-perturbative problem. It was noted in [29,30] that, in integrable models, this constant can be calculated exactly. The reason is as follows: in the calculation of F (h) JHEP04(2020)160 from the Bethe ansatz, the answer is naturally expressed in terms of m/h, where m is the mass of the charged particles under Q. On the other hand, the perturbative calculation gives the answer in terms of Λ/h. By matching these two expressions, one can find an exact expression for the physical mass as a function of Λ. This typically requires just a one-loop calculation in the "bosonic" theories and a two-loop calculation in the GN model. The original calculation of [29,30] was done for the non-linear sigma model, but it was quickly generalized to the PCF [31,34], the Gross-Neveu model [32,33], and to supersymmetric models [35,36] (see [37] for a review). In order to perform these computations, one has to solve the integral equation (2.8) for large h, which corresponds to large B. This is technically challenging and it was done in [29][30][31][32][33][34][35][36] by using the Wiener-Hopf method. In the next section, following [38], we will present a more powerful method to solve the integral equation, which can be used to generate the full perturbative series for F (h), or, equivalently, for the ground state energy e(ρ).
Before doing this, let us list some of the basic ingredients in the four theories that we will study.
(i) Non-linear O(N ) sigma model. The basic field of the non-linear sigma model is a scalar field S : R 2 → R N satisfying the constraint S 2 = 1. (2.12) The Lagrangian density is where g 0 is the bare coupling constant. The first two coefficients of the beta function are [45] β 0 = 1 4π∆ , 15) and the coefficient ξ defined in (2.2) is given by This theory has a global O(N ) symmetry. The conserved currents are given by and we will denote by Q ij the corresponding charges. As in [29,30], we take as our charge in (2.3) the quantum version of Q 12 . The exact S matrix of the O(N ) non-linear sigma model was found in [46]. For particles charged under Q 12 , it is given by and ∆ is given in (2.15).

JHEP04(2020)160
(ii) N = 1 non-linear O(N ) sigma model. The N = 1 supersymmetric version of the non-linear sigma model has two fields: the field S of the purely bosonic version, satisfying also S 2 = 1, and an N -uple of Majorana fermions ψ satisfying the constraint The Lagrangian density is where we follow the conventions of [42]. The first two coefficients of the beta function are (see [35] and references therein) where ∆ is given in (2.15). It follows that The exact S matrix of the supersymmetric O(N ) sigma model was obtained in [47]. The precise choice of conserved charge in (2.3) is discussed in detail in [35]. One eventually obtains an integral equation of the form (2.8) in which the kernel is given by and ∆ is again as in (2.15).
(iii) SU(N ) principal chiral field. Here, the field is a map Σ : R 2 → SU(N ), with Lagrangian density The first two coefficients of the beta function are [48] β 0 = 1 16π∆ , In this case, The S matrix of the principal chiral field was obtained in [49,50]. The choice of conserved charge in (2.3) is as in [31], and the relevant S matrix element is

JHEP04(2020)160
(iv) SU(N ) principal chiral field with FKW charges. The SU(N ) principal chiral field can also be explored by using a different set of conserved charges, discussed in [40,41], which despite exciting multiple particles can be more convenient for certain large-N analysis. The kernel in (2.8) is presented in [41]: , (2.30) where ∆ is given by (2.27).
and we follow the conventions of [32,33]. The first two coefficients of the beta function are (see e.g. [51]) where ∆ is given in (2.15). Therefore, The full S matrix of the GN model was found in [46]. As in [32,33], we take as our charge in (2.3) the quantum version of Q 12 , associated to the global O(N ) symmetry. The relevant S matrix is then where x and ∆ are again given by (2.19) and (2.15), respectively.

General solution of the integral equation
As it is clear from the discussion above, the perturbative regime of the integrable field theory corresponds to large h, which means large B in the integral equation. This is a singular limit which is difficult to study analytically. This problem appears already in much simpler models solved by the Bethe ansatz, like the Lieb-Liniger [52] and the Gaudin-Yang [53,54] models. In a tour-de-force paper [38,39], Volin reformulated in a powerful way the matching method which was used to study the integral equation of the Lieb-Liniger model [55,56]. He applied this method to the non-linear sigma model and he was able to compute analytically the perturbative series to large order. We have recently generalized Volin's method to solve the long-standing problem of deriving the perturbative series in the Lieb-Liniger and the Gaudin-Yang models from the Bethe ansatz [13,14]. In this paper we will generalize it to the integrable quantum field theories listed above. We will obtain in this way analytic results for the solution of the integral equation (2.8) when B is large, and in particular we will find explicit expansions for JHEP04(2020)160 ρ, e(ρ) in power series of 1/B and log B. In particular, we will streamline Volin's method and derive one of its key ingredients directly from the Wiener-Hopf decomposition of the kernel.
A crucial ingredient in [38,39] (see also [57]) is the resolvent of the density of Bethe roots, This function is analytic in the complex θ-plane but it has a discontinuity in the interval [−B, B], given by From its definition we deduce that It follows from (2.10) that we can compute the density ρ as a function of B from the residue at infinity of the resolvent. The weak coupling regime corresponds to large B, so we should study the resolvent in a systematic expansion in 1/B. To do this, we consider the resolvent in two different regimes. The first one is the so-called bulk regime, in which we take the limit in such a way that is fixed. This is therefore appropriate to study χ(θ) near θ = 0. The second regime is the so-called edge regime, in which we also have (2.38) but we keep fixed the variable This is therefore appropriate to study χ(θ) near the edge of the distribution θ = B.
In order to study the bulk regime, we need an appropriate ansatz for R(θ). We will propose concrete formulae for the relevant models in the next section, following previous studies in [13,14,38,55,56,58]. To study the edge regime, we write R(z) = R(θ(z)) as a Laplace transform, A general property ofR(s) is that, as explained in [13,38,39], it has an expansion as s → ∞ in integer, negative powers of s: where χ n , n ≥ 0, are the coefficients of the Taylor expansion of the density of eigenvalues near θ = B,  Figure 1. The contour C in the integral (2.46) can be deformed to (minus) the Bromwich contour B appearing in the inverse Laplace transform.

q 6 l i k b c + N n 8 3 i k 5 s 8 q A h L G 2 p Z D M 1 d 8 T G Y 2 M m U S B 7 Y w o j s y y N x P / 8 7 o p h j d + J l S S I l d s s S h M J c G Y z J 4 n A 6 E 5 Q z m x h D I t 7 K 2 E j a i m D G 1 E J R u C t / z y K m l d V L 3 L q n t / V a n V 8 z i K c A K n c A 4 e X E M N 7 q A B T W A g 4 R l e 4 c 1 5 d F 6 c d + d j 0 V p w 8 p l j + A P n 8 w e j C o + y < / l a t e x i t >
The connection between (2.42) and (2.43) makes it possible to test the solution forR(s) against a numerical study of the distribution χ(θ) at the edge. One consequence of (2.42) is that, sinceR(s) decreases as 1/s at infinity, one can use the Bromwich inversion formula to writeR The functionR(s) can be also used to calculate the energy as a perturbative series in 1/B and log B, as pointed out in [38]. To see this, let us neglect exponentially small contributions of the form e −B in the first equation of (2.10). Then, we can write where the contour C is shown in figure 1. It can be deformed to (minus) the Bromwich contour, and we obtain at the end of the day up to exponentially small corrections at large B. It follows from the above considerations that, if one knows the resolvent R(θ) and its inverse Laplace transformR(s), it is possible to calculate the functions ρ, e as a function of B. To determine R(θ), there are two possible routes. In the first one, followed in [38], one writes the kernel K(θ) in the form where O is a difference operator. The integral equation can then be written as a difference and discontinuity equation for the resolvent, This can be used in the edge regime to determine the analytic structure ofR(s), which makes it possible to derive its functional form. However, in [13] it was shown, based on observations in [32,38], that one can findR(s) directly from the Wiener-Hopf decomposition of the kernel K(θ). Let us present the argument in [13], adapted to the setting of this paper. We first use the variable (2.40) to inspect the edge limit of (2.8). It is useful to define χ(z) = χ(θ(z)) and By extending the above equation to all real z, and ignoring exponentially small terms when B is large, we obtain where ξ(z) is an unknown function. Let us then define one last Fourier transform so that we can transform the full equation: The subscripts ± denote that something is analytic in the upper/lower half complex plane (including the real axis but possibly excluding the origin). We also introduce the Wiener-Hopf decomposition of the kernel is even (as it happens in all the cases we will consider). This decomposition can almost always be done provided 1 − K(ω) is well defined along the real axis (see e.g. [32]), though some care might be necessary at ω = 0. With the above definitions, we take the Fourier transform of (2.53) and obtain

JHEP04(2020)160
We can rewrite (2.56) as From (2.57) it follows that C(ω) must be analytic in both the upper and lower half complex planes, reducing it to an entire function. However, in their respective half planes (including the real axis), F − (∞) = 0, X ± (∞) = 0 and G ± (∞) = constant (where the first two come from their definitions and (2.53), while the latter can be checked explicitly). Both sides are thus bound at infinity, and by using Liouville's Theorem we find that C(ω) = 0. We conclude that, at leading order in the edge limit, We can now relate this function to the inverse Laplace transform of the resolvent (2.44). We considerR(is) for Im(s) < 0 and bend the Bromwich contour around the negative real axis, without crossing it, as in figure 1. We obtain This is the leading solution in the strict large B limit, but one has additional corrections as a power series in 1/B. Taking into account constraints on the allowed poles and behaviour at infinity ofR(s), we can writê and Q(s) is of the form The result (2.59) is of the form obtained by Volin in [38] in the case of the non-linear sigma model. In our derivation, the function Φ(s) is determined simply by the Wiener-Hopf decomposition of the kernel, therefore it can be written immediately in a large number of cases. The coefficients Q n,m appearing in this expansion are not fixed by the Wiener-Hopf decomposition, but as shown in [38,55,56], they can be calculated recursively by comparing the edge solution (2.59) to the ansatz for R(θ) in the bulk regime, see [13,38] for details on the matching procedure. JHEP04(2020)160

Solving the bosonic models
As we mentioned before, we will refer to the non-linear sigma model, its supersymmetric extension, and the PCF, as "bosonic" models. In the bosonic models the perturbative expansion of the free energy has the structure where κ i , i = 0, 1, 2 are calculable constants andḡ is the RG-invariant coupling defined in (2.5). In order to write down the perturbative series, it is useful to follow [59] and introduce an intermediate coupling α defined as 2 where ξ is given in (2.2). Similar schemes in QCD have been considered in [60,61]. We have, at leading order, This suggests to introduce yet another coupling constant α through the relation where c is a model-dependent constant. As we will see, the quotient e/ρ 2 can be expressed as a formal power series in α, without logarithms. A wise choice of c simplifies the resulting expressions (removing e.g. terms involving log (2)).
Let us now consider the solution of the integral equation for the bosonic models. We have, as noted in [31] in the case of the PCF, where k, η, a are constants. φ(s) is typically a rational function of Γ-functions, and there are no log(s) on the r.h.s. before expanding the last exponential. An interesting observation is that η does not affect the final result for the energy expansion. In [31] it is noted that we must have a = 1 − 2ξ. It is convenient to normalize We re-define A to keep the form ofR(s) unchanged:

JHEP04(2020)160
From this formula we extract directly, by using (2.47), In order to fix the remaining coefficients, we need an ansatz for the resolvent in the bulk.
(iv) SU(N ) principal chiral field with FKW charges. From (2.30) we retrieve (2.86) Due to the existence of multiple particles we must change the overall factors in the definitions of the observables, though we keep the definitions ofẽ andρ from (2.69) and (2.71) respectively, 2ẽ .

(2.87)
As in the standard PCF case discussed above, only the c n,m,0 coefficients are non-zero. This is similar to the analysis done in [62]. From now on, we will remove the subscript FKW in e, ρ and work exclusively with the quantities defined in (2.87).
We introduce the coupling constant defined in (2.65) and we choose c to be where γ E is the Euler-Mascheroni constant and ψ (m) is the polygamma function. With

JHEP04(2020)160
this choice of α we find e Bρ 2 = α + For compactness we have introduced the auxiliary function We have calculated the first 50 terms of the expansion (2.89).
It is instructive to test the above expansion against one-loop perturbation theory (in the other models compared in this paper, this comparison has been done in [29][30][31][32][33][34][35] to derive the mass gap). In order to do this, we have to compute the free energy as a function of the external field h. We find . (2.91) We have to express first α in terms of h/m: By using all these results, one eventually finds, (2.94)

JHEP04(2020)160
Let us now compare this result to the one-loop expression for the free energy. After using the relationship found in [31] between m and the dynamically generated scale Λ MS , where the charges q i are given by [40] the h-dependent term of (2.94) matches the one in (2.96). In order for the h-independent terms to match, the following equality must hold (2.99) We do not have an analytic proof of this identity, but we have checked it for many values of N .
In [40] and more recently in [62], the PCF with FKW charges has been studied in the large N expansion. Although we work at finite N , it is straightforward to expand our results in series in 1/N , and we find for example which agrees with the results in [40,62].

Solving the Gross-Neveu model
In the case of the GN model, the free energy F (h) was calculated in perturbation theory in [32], and it reads, and g is the RG invariant coupling constant defined by (2.5). One can introduce an intermediate coupling exactly as in (2.63), but in order to simplify the perturbative expansion it is better to use the analogue of (2.65), which in the GN case it is given by (2.108) We also need an ansatz for the resolvent in the bulk. In the GN case, we have where e(k) = 0 if k is odd and 1 if k is even. This has a different structure from the bosonic case, but not surprisingly it is similar to the bulk ansatz used for the Gaudin-Yang model in [13,14]. From the bulk and the edge ansatz, we obtain the following expressions for ρ and e:

JHEP04(2020)160
As in the bosonic models, the unknown coefficients can be obtained by matching the bulk and edge answers. One obtains in this way the perturbative expansion 4ẽ We have calculated analytically the first 45 terms of this expansion. It is also possible to analyze the integral equation (2.8) in a 1/N expansion, as already noted in [32,33]. In this way one can obtain all-order results in α for the contribution of order ∆ to the above perturbative expansion. This is explained in appendix A. The result (A.12) provides a further test of (2.112).

Large order behavior from renormalons
In many asymptotically free theories, renormalons are expected to dominate the large order behavior of conventional perturbation theory. However, as we mentioned in the Introduction, diagrammatic arguments are not fully conclusive and it is useful to have another type of reasoning which leads to a precise expectation for the large order behavior of perturbation theory. One such argument is provided by the connection between IR renormalons and the OPE, first pointed out by Parisi [11]. Let us briefly review this argument (see e.g. [5,17,60]). A generic observable in an asymptotically free QFT can be written as the sum of a perturbative and a non-perturbative contribution: is the perturbative series, and ϕ np (g) is typically exponentially small in the coupling constant g. In the terminology of the theory of resurgence (see [5][6][7] for reviews) we say that ϕ(g) is given by a trans-series with two different small parameters, namely g 2 and where A is an appropriate constant. In some cases, the observable ϕ(g) can be studied with an OPE, and this determines the form of ϕ np (g). 3 The OPE will involve contributions of JHEP04(2020)160 a series of operators O i , of dimension d i . Let us focus in the following on the contribution of a single operator O of dimension d. It is of the form, where Q is an external scale (it could be the external momentum in an Adler function, or the chemical potential h in the situation considered in this paper). In (3.4) we have indicated explicitly the dependence on the renormalization scale µ, and C(Q/µ, g) can be computed from perturbation theory. Since both ϕ(g) and ϕ p (g) are separately RG-invariant, the same must happen to ϕ np (g). Using standard RG arguments (see e.g. [17,60]), and evaluating the non-perturbative correction at µ = Q, we find is the anomalous dimension of O, and our convention for the β function is as in (2.1). In (3.5), g 0 is a reference coupling. At leading order in g we find, In (3.5) and (3.7) we have absorbed overall constants in C (1, g(Q)). As anticipated, ϕ np (g) involves an exponentially small parameter of the form (3.3), where Let us assume that where n 0 is a non-negative integer. By a standard argument [5][6][7], the exponentially small correction (3.7) due to the condensate of O gives the following contribution to the large order behavior of the coefficients a n appearing in (3.2):

JHEP04(2020)160
An equivalent statement of this relation can be obtained by considering the Borel transform of ϕ p (g), defined as ϕ p (ζ) = n≥0 a n n! ζ n . (3.13) Then, the large order behavior (3.11) means that each operator of dimension d appearing in the OPE gives a singularity in the Borel plane located at ζ = A. This is usually called an IR renormalon singularity. Note that the location of the IR renormalon singularities gives information on the one-loop coefficient of the beta function, and on the operators contributing to the OPE. In addition, the coefficient b + appearing in (3.12), which determines the next-to-leading correction to the leading asymptotics, has information on the two-loop coefficient β 1 of the beta function and on the anomalous dimension of the corresponding operator.
There is however another source of large order behavior in QFT: the UV renormalons. In the case of asympotically free theories, they lead to terms in the trans-series with positive exponents. These terms can be related to operators of dimension d + D, where D is the dimension of spacetime [10,17,63]. Their contribution to the large order behavior is of the form Here, A is as in (3.9), and where m 0 is a non-negative integer. In general, in an asymptotically free theory, we will have both IR and UV renormalons. Let us label the operators leading to IR renormalons by the indices i ∈ I IR , and the operators leading to UV renormalons by j ∈ J UV . Then, the perturbative series has to be extended to a general trans-series with exponential terms of the form i∈I IR (3.16) Here, b + i , b − j are given in (3.12) and (3.15), respectively, where we set d = d i,j and γ (1) = γ (1) i,j . As a consequence, we find the following large order behavior for the perturbative series: a n ∼ 1 2π where Here, we have restricted ourselves to the next-to-leading order for the asymptotic expansion. Further corrections can be also studied (see e.g. [17]), but we will not consider them in this paper.

Testing renormalon predictions
The leading large order behavior of the perturbative series will be determined by the IR and UV renormalons with the smallest possible value of d i . In two dimensions, there is an IR singularity in the Borel plane at which corresponds to condensates of dimension d = 2. There is also an UV singularity at ζ = −β −1 0 , due to operators of dimension 4 (see e.g. the discussion of [22] on the non-linear sigma model). Let us label the dimension d = 2 operators contributing to the IR singularity by i ∈ I (2) IR , and the dimension d = 4 operators contributing to the UV singularity by j ∈ J (2) UV . If we define we obtain the following renormalon prediction for the large order behavior: In this formula we have redefined the overall constants to absorb the factor (2π) −1 . Both the UV and the IR renormalon contribute to (3.21), and in order to test this prediction we have to disentangle their contribution. One obvious consequence of the presence of both singularities in the Borel plane is that odd and even terms of the c n series have different large order behavior. We have where k 1. Let us now introduce the auxiliary series: Then, we have the large order behavior, JHEP04(2020)160 we conclude that Therefore, the first sequence is sensitive to the first IR renormalon singularity, while the second sequence is sensitive to the first UV renormalon singularity.
We would like to test the above expectations from the theory of renormalons, against the large order behavior of the perturbative series that we have calculated above. Although the arguments we have presented are derived for observables in which the OPE can be used, we will assume that they also control the ground state energy in the presence of an external field studied in this paper. We first note that in our perturbative series the coupling is defined by (2.65) and (2.103) in the bosonic and GN models, respectively. In both cases we have α ∼ 2β 0 g 2 . (3.28) As discussed in [17,64], such redefinitions of the coupling constant change the location of the singularity in the Borel plane by the overall factor 2β 0 . This means the leading IR and UV singularities will be at ζ = ±2. However, they do not change the values of the coefficients b + i , b − j . In particular, the normalized version of the series (3.20) involves multiplying by 2 n the coefficients of the perturbative series obtained in section 2.
We will mostly focus on the IR renormalon singularity. This is because it involves dimension 2 operators, and there is only a small number of these. A detailed study of the UV renormalon singularity requires all the dimension 4 operators. For example, in the non-linear sigma model there is only one dimension 2 operator, but five dimension 4 operators [65,66]. To extract information about the leading IR singularity, we note that the leading large order behavior of S k is governed by the largest values of the b + i , which we will denote by b + * . This coefficient can be extracted from the auxiliary sequence which behaves as Finally, we note that the perturbative series for the bosonic models start with α, and not with α 0 . This simply amounts to a redefinition of the integer n 0 appearing in the asymptotic formulae listed above. We will now present a discussion of the four models we consider in this paper. (i) Non-linear O(N ) sigma model. A preliminary analysis of the large order behavior of the perturbative series (2.77) model was performed in [38]. A direct test of the presence of singularities in the Borel plane can be simply done by plotting the poles of the Borel-Padé transform of the perturbative series (2.77). As we see in figure 2 in the case of N = 4, these poles accumulate along branch cuts in the positive and the negative real axes. The location of the branch points is clearly at ζ = ±2, signaling the first IR and UV renormalons. The first IR renormalon corresponds to the operator of dimension 2 given by (see e.g. [22])

JHEP04(2020)160
This is the analogue of the gluon operator in Yang-Mills theory. The anomalous dimension of the operator given by the Lagrangian density is closely related to the beta function of the coupling constant (see e.g. [67,68]). One finds, so that γ (1) = 2β 0 . (3.33) In particular, in the asymptotics (3.26) there is a single term with In order to test this prediction, we use our data for the perturbative series to construct the sequence (3.29), and we study its asymptotic behavior for different values of N . To remove tails and to accelerate the convergence to b + * , we can use Richardson transforms (see e.g. [5]). Our results vindicate the result (3.30) with b + * given in (3.34), and n 0 = 0. Note that, by consistency, this provides a test not only of the value of b + * , but also of the assumptions leading to the prediction (3.29), namely the existence of an IR renormalon singularity at ζ = 2, and the factorial growth of the sequence. As an example, in figure 3 we show the sequence (3.29) and its second Richardson transform for N = 4 (left) and N = 6 JHEP04(2020)160 (right). After two Richardson transforms, the sequences approach the value b + * −1 = 2∆−2 with an error of 10 −5 and 10 −4 , respectively. We have tested this for many other values of N . The behavior persists even for rational values of N , although the tests become less precise as N becomes large.
It is also possible in this case to analyze the UV renormalon singularity at ζ = −2 in some detail, by looking at an auxiliary sequence similar to (3.30), but where we replace S k by D k . Our results indicate that the asymptotic behavior is controlled by an exponent The first two values are compatible with our empirical finding (and appropriate values of m 0 , namely m 0 = 0 and m 0 = 1, respectively). If this picture is correct, our numerical result verifies the fact that, in the UV renormalon, the next-to-leading asymptotics (3.15) involves the quotient −β 1 /(2β 2 0 ) (i.e. with the opposite sign than the IR renormalon).
(ii) N = 1 non-linear O(N ) sigma model. The behavior of this model is quite different from the other ones we are considering, since the IR renormalon singularity at ζ = 2 is absent. This can be seen by plotting the poles of the Borel-Padé transform of the perturbative series, as shown in figure 5 in the case of N = 5. Therefore, the only contribution to the leading asymptotics comes from the leading UV renormalon, and one finds a n ∼ (−2) −n Γ(n − 1). (3.37) The absence of the leading IR renormalon contribution to the asymptotics is expected from the arguments in [43,44], where it was shown that, in supersymmetric theories, many condensates vanish. In particular, in a composite superfield, only the lowest component can   (20,20). The poles accumulate in one branch cut starting at the singularity ζ = −2, corresponding to the first UV renormalon. There is no IR renormalon singularity at ζ = 2.
have a non-zero vev. This rules out for example a condensate of the operator of dimension 2 given by the Lagrangian density itself, (2.21).
(iii) SU(N ) principal chiral field. The PCF is very similar to the non-linear sigma model. The first IR renormalon corresponds to the operator of dimension 2 appearing in the Lagrangian density: Its anomalous dimension is given again by the formulae (3.32), (3.33). There is a single term in (3.26) with Note in particular that this is independent of N . Our numerical calculations, using for example the poles of the Borel-Padé transform, indicate clearly the presence of IR and UV renormalons at ζ = ±2. Moreover, they vindicate the value (3.39) with n 0 = 0. In figure 6 we show the sequence (3.29) and its second Richardson transform for N = 2 (left) and N = 3 (right). After two Richardson transforms, the sequences approach the value JHEP04(2020)160   (iv) SU(N ) principal chiral field with FKW charges. As one would expect, despite the different perturbative series, the asymptotics of the coefficients remains the same since the renormalon physics should be unaffected by the charge choice. Indeed, as shown in figure 7, an identical analysis finds b + * − 1 = −1 with an error of order 10 −4 for N = 2 and N = 7, for example.
The calculation of their anomalous dimensions can also be easily done, by using e.g. the results of [51] on the renormalization properties of this model. One finds, at one-loop, γ This result agrees with the calculation in [69]. Therefore, in the asymptotics (3.26) there will be two terms. In the first one, corresponding to the quartic fermion term, one has while in the second one, corresponding to the kinetic term, According to our numerical results, the leading term in the asymptotics has b + * = −2∆. This corresponds to (3.42) with n 0,1 = 1 or to (3.43) with n 0,2 = 0. In figure 8 we show the sequence (3.29), as well as its second Richardson transform, for N = 2 (left) and N = 3 (right). The numerical result agrees the expected result b + * − 1 = −1 − 2∆ with an error of 10 −4 and 10 −3 , respectively.
We conclude that the large order behavior of the perturbative series obtained in section 2 are indeed controlled by renormalons (provided some reasonable assumptions are made on the values of the exponents n 0 or m 0 ). In particular, we can isolate the contribution from IR and UV renormalons separately. This makes it possible to test the next-to-leading contribution of the IR renormalon to the asymptotics in the four models considered in this paper, which is given by (3.12). This value is sensitive to the first two coefficients of the beta function. In practice, this is clearly seen in our numerical analysis when we change N : in the O(N ) sigma model and in the GN model, this coefficient involves ±2/(N − 2), respectively, while in the PCF it is independent of N .

Conclusions
Since the pioneering work in [28], it is well known that in integrable field theories in two dimensions the ground state energy can be calculated exactly once a conserved charge is coupled to an external field. Extracting the perturbative series from the resulting Bethe ansatz equation turns out to be challenging. In this paper we have built upon [38,39] to JHEP04(2020)160 obtain the perturbative series expansion in a number of integrable field theories. In this way we have been able to provide a direct test of renormalon predictions for the large order behavior of these series. In particular, we have tested the next-to-leading correction to its asymptotics, involving the first two coefficients of the beta function, as well as the anomalous dimensions of the operators appearing in the OPE. We have also observed that, in the supersymmetric version of the non-linear sigma model, the first IR renormalon is absent, in accord with the observations of [44].
Although we believe that our tests of renormalon predictions are convincing, there are various issues that should be clarified. For example, we used predictions about the large order behavior based on OPE considerations. Our observable involves perturbing the Lagrangian with an integrated conserved current, so it seems reasonable that the OPE can be used to establish our working hypothesis, but one should address this point more carefully. Also, in our tests we made some additional assumptions on the exponents n 0 , m 0 appearing in (3.10) and (3.15). It would be interesting to derive these assumptions from first principles.
Some technical aspects of our analysis could also be improved. For example, although we have produced 40-50 coefficients in the different perturbative series, with the current implementation of the method there is a computational bottleneck beyond 50 terms. If we were able to generate many more terms, it is very likely that we could improve our tests of IR renormalon behavior, which become less precise as N becomes large. Clearly, our methods could be straightforwardly applied to other integrable field theories.
There are more ambitious research directions open by our results. In this paper we have extracted the perturbative series from the Bethe ansatz equation, but we have neglected exponentially small terms at large B. In principle one could incorporate these terms in a systematic way to obtain the full trans-series (3.16) associated to renormalons. In other words, one should be able to express the exact e, ρ as Borel-Écalle resummations of suitable trans-series, extracted from the integral equation (2.8). Such a trans-series solution would give an enormous wealth of information on the renormalon physics of these field theories, including exact values for the condensates (modulo Stokes jumps).
Another interesting direction is to understand the fate of renormalons after a (twisted) compactification on a circle. It has been suggested that, after such a compactification, renormalons disappear as such [70] and the corresponding singularities can be realized semiclassically [71][72][73][74][75][76][77][78] (see [79] for related work). In particular, [74,75,77] provide a concrete semiclassical picture in two of the models considered in this paper, namely the principal chiral field and the non-linear sigma model. It would be interesting to use integrability techniques to compute observables in the twisted compactification of these theories, as a function of the compactification radius. In this way one could study in detail the behavior of renormalon singularities and their eventual transmutation in semiclassical instanton singularities.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.