Resurgence in the O(4) sigma model

We analyze the free energy of the integrable two dimensional O(4) sigma model in a magnetic field. We use Volin's method to extract high number (2000) of perturbative coefficients with very high precision. The factorial growth of these coefficients are regulated by switching to the Borel transform, where we perform several asymptotic analysis. High precision data allowed to identify Stokes constants and alien derivatives with exact expressions. These reveal a nice resurgence structure which enables to formulate the first few terms of the ambiguity free trans-series. We check these results against the direct numerical solution of the exact integral equation and find complete agreement.


Introduction
Perturbation theory has proved to be a useful tool in calculating physical processes for the electromagnetic and weak interactions, but has had only a limited success for their strong counterpart. Important phenomena such as confinement and dynamical mass generation are inherently non-perturbative, and cannot be accessed from the few known perturbative coefficients. Perturbation theory in QCD is expected to be asymptotic, with coefficients growing factorially, see e.g [1,2]. This factorial growth can be traced to the proliferation of Feynman diagrams [3,4], or to integrals in specific renormalon diagrams, when loop momenta lie in various IR and UV domains [5]. It is a signal of non-perturbative contributions, which usually originate from non-trivial saddle points in the path integral. Because these are exponentially suppressed, they do not appear directly in perturbation theory, but they can be extracted from the large-order behaviour of perturbative coefficients. The tools for doing this are known as resurgence theory [6,7,8].
While it would be ideal to apply this theory to non-perturbative phenomena in QCD, too few perturbative coefficients are known for this to make progress. Thus we turn to toy models, which share important features with QCD, but are more tractable. The two dimensional O(N ) symmetric sigma models are exactly of this type: They exhibit a dynamically generated mass gap, and are asymptotically free in perturbation theory [9]. At the same time they are integrable, allowing physical quantities such as the mass gap, scattering matrices, and ground state energy to be calculated exactly [10,9]. Our aim in the present paper is to use the O(4) sigma model to reveal the relation between the perturbative and non-perturbative effects, as the first steps in the full resurgence program.
The simplest cure for the factorial growth of perturbative coefficients is the Borel transform, which is obtained by dividing out this factorial in order to ensure constant asymptotics. The resulting function has a finite radius of convergence, and exhibits pole and cut singularities, often on the real line. Since the inverse Borel (Laplace) transformation involves an integration of the analytically continued function along the positive real line, singularities there lead to ambiguities in the result, as the contour must be shifted. Such cases are often called non Borel summable, but in fact these ambiguities, which are imaginary and exponentially small, contain useful information. While isolated poles give single exponential terms, each cut gives an exponential multiplied by a power series, which is itself asymptotic. The coefficients of the original series are said to resurge into those of these further series, which describe the contribution of non-trivial non-perturbative saddle points. The maps between these different sectors are called alien derivatives, and their algebraic properties provide useful constraints. The full physical expression is what is called a trans-series, a sum over different exponential factors each multiplying a real-valued power series, which is the final result.
The resurgence program has been pushed forward for several quantum field theories including supersymmetric theories [11,12,13,14,15,16], various quantities in the maximally supersymmetric 4D gauge theory in the large N limit [17,18,19], the simple φ 4 theory in 2D, [20,21]. Factorial growth can be seen in lattice simulations [22]. Much is known in the large N , [23,24] and also semiclassically [25,26,27]. They also appear in the hydrodynamics of the Yang Mills plasma [28,29]. Recently in many body systems the authors of [30,31] could manage to extract the leading exponentially small corrections exactly in the Gaudin-Yang model. In the Hubbard model the knowledge of the exact perturbative coefficients enabled them to calculate all the alien derivatives [32]. We are not aware of any asymptotically free QFT, however, where the full resurgence theory had been rigorously established.
The O(4) sigma model could be the first example, as it allows an exact treatment. This model was one of the first where the scattering matrix was exactly determined [10]. It was also the first model where the dynamically generated scale Λ M S was analytically related to the mass of the particles m [9]. This seminal calculation was done by introducing a magnetic field coupled to one of the conserved O(N ) charges, and determining the free energy in two different ways. For large magnetic fields one can establish a standard, renormalization group improved perturbative expansion in the parameter h/Λ M S . On the other hand, the magnetic field forces the positively charged particles to condense into the vacuum, and this vacuum condensate consists of particles with rapidities in (−B, B), with density χ(θ). Based on the scattering matrix the thermodynamic limit of the Bethe ansatz leads to a linear integral (TBA) equation for χ(θ), which determines the density ρ and the groundstate energy , whose Legendre transform is the sought for free energy. By neglecting exponentially small contributions the TBA equation can be expanded as a function of h/m. Comparison to the perturbative expansion then led to the relation between the UV parameter Λ M S and the IR parameter m. We emphasize that the TBA equation is exact, in that it contains also all the non-perturbative exponentially small corrections. Thus its analytical expansion could lead directly to the exact trans-series of the groundstate energy and free energy, which would provide a veritable gold mine for the resurgence literature. Unfortunately, the exact calculation of the exponentially small terms is beyond the scope of the present day research. Even the calculation of not just the first few perturbative coefficients resisted an analytical treatment for decades.
A breakthrough was obtained by Volin, who invented a way to calculate the perturbative coefficients systematically [33,34]. His idea was to expand the resolvent of the spectral density χ(θ) both in the middle of the interval θ ∼ 0 and in the edge region θ ∼ B, and then to match the two expansions. In the middle region the TBA equation determines the analytical structure of the resolvent, while in the edge region the Wiener-Hopf technique can be used to parameterize its Laplace transform. Surprisingly, matching the two representations fixes all the perturbative coefficients in terms of algebraic equations, which can be solved iteratively. Using this method Volin calculated the first 26 perturbative coefficients for generic O(N ) models, which was later extended to 44 coefficients in [35]. These authors also extended the method for other relativistic and non-relativistic theories [30,31]. The perturbative coefficients in the O(N ) models are linear combinations of products of odd zeta functions degrees not larger than the perturbative order. Although this information allowed the authors to gain qualitative information about the analytical behavior of the Borel transform, and establish the location of the leading singularities [33,35], it is not sufficient to investigate resurgence properties. Focusing on the O(4) model, we managed to solve the algebraic equations recursively in a closed form. This enabled us to calculate the first 50 perturbative coefficients exactly, and switching to work numerically, the first 2000 terms with 12000 digits precision. Numerical data is equally useful for studying resurgence, and with this data, we have been able to fix the first few terms in the trans-series. On the way we observed a very interesting resurgence pattern between the physical observables, which deserves to be investigated further, and we hope that our results will spark new research in this field. We now summarize our method and results.

Summary
We started to investigate the perturbative expansion of the groundstate energy as the function of the running renormalized coupling. Having observed factorial growth we switched to the Borel transform and got insight into the analytical structure by its Padé approximant. We observed a pole singularity at 1, a cut starting at 2 and another cut starting at −1, signaling both UV and IR renormalons. In order to get a more precise analytical continuation of the Borel transform we applied the conformal mapping method. We then used the inverse of the Borel transformation by integrating a bit above and below the real line, so avoiding the singularities. The results had unwanted imaginary parts, being the complex conjugate of each other. Since we wanted to understand how precise the real part was we analyzed the TBA equations directly. As the TBA equation provides an exact answer, we solved it numerically with very high (30-50 digits) precision. Comparison with the real part of the inverse Borel transformed revealed exponentially suppressed corrections. Our aim was to understand these non-perturbative imaginary and real deviations directly from the perturbative corrections, thus to establish the first steps into resurgence.
In doing so we used the fact that the leading singularity on the Borel plane is encoded in the large n behavior of the perturbative coefficients χ n . Constant asymptotics determines the residue of the pole singularity, called Stokes constant, while consecutive 1/n, 1/n(n − 1), corrections provides the perturbative expansion of the function multiplying the logarithmic cut, called the alien derivative of the original functions, ∆ ω χ, where ω is the position of the cut. The available large number of precise perturbative coefficients enabled us to determine these coefficients with high (more than 100 digits) precision. By using simple assumptions on the Stokes constants, (ratios of powers of e and π) and the structure of the alien derivatives (have the same transcendental structure as the original perturbative coefficients, involving products of odd zeta functions with increasing transcendentality) enabled us to guess the first ≈ 9 coefficient exactly and determine the next 60 with reasonable but decreasing precisions. The explicit knowledge of the pole term at 1 and the leading perturbative coefficient of the first cut starting at 2 completely agreed with the imaginary ambiguity of the inverse Borel transform. In order to cancel the ambiguity we had to add these exponentials to the perturbative series. Since these new non-perturbative terms are asymptotic by themselves we had to analyze the Borel transform of the function multiplying the logarithmic cut starting at 2. Similar asymptotic analysis revealed a cut starting at −1 and another one starting at 2. The ambiguity coming from the former one in the inverse Borel transform was real and its canceling exponential left a real contribution, which seemed to agree numerically with the deviation we observed in the comparison with the TBA result. This was a very reassuring sign, but in order to establish a more precise matching and to reveal the full trans-series we needed to relate the various alien derivatives to the original perturbative series. Since we could not recognize any resurgence structure of the free energy as the function of the renormalized coupling at this point we started to analyze the groundstate energy, and the density ρ as functions of the TBA variable B.
Our asymptotic analysis revealed that ρ has a cut starting at −1 and another one starting at 2, while has the same structure with an additional pole at 1. For our big surprise we managed to relate their alien derivatives at −1 with themselves, i.e. the alien derivative of any of these functions is proportional to itself multiplied by . This enabled us to calculate all consecutive actions of ∆ ±1 , which can be expressed in terms of and ρ, that is they resurge to themselves. We called these functions first generation. We could not relate, however, the new functions appearing by their alien derivatives ∆ 2 in any way so we called those functions second generation. These second generation functions resurge also themselves, once we calculated their ∆ ±1 alien derivatives, but their alien derivatives at ∆ 2 are again independent. We called these new functions the third generation. Our original numerical precision allowed as to see the structure only up to this point, but already there a beautiful structure appeared, which definitely deserves a better understanding and a quantitative description, which might be obtained by analytically calculating the trans-series expansion directly from the TBA equation.
By using the formula for the alien derivatives of composite functions we managed to calculate the resurgence properties of the free energy as a function of the renormalized coupling. These formulas allowed us to express the D ±1 alien derivatives (in the running coupling) of the free energy in terms of first generation functions, while the D 2 derivatives with first and second generation functions. These allowed us to calculate the expression D 2 D 2 f , which contributes to the leading real deviation from the TBA results. Using the median resummation prescription based on the alien derivatives D 1 f , D 2 f and D 2 D 2 f we could reproduce both the imaginary and real deviations from the TBA result, i.e. we reproduced the real physical value including non-perturbative exponentially suppressed terms purely from the perturbative coefficients. These completed the first steps in the resurgence program and provided the leading terms in the trans-series ansatz, which we also formulated. Unfortunately we cannot see yet, how any bridge equation could be derived, which would relate the functions from different generations to each other.
In summarizing, by determining a large number of high precision perturbative coefficients in the O(4) model we managed to extract non-perturbative information and construct the first few terms in the trans-series exactly. These results are in complete agreement with numerical data obtained from the conformal mapping method and the direct numerical solution of the TBA equation.
Our results provide the location of the first few non-trivial saddle points together with the exact perturbative coefficients coming from fluctuations around them. It would be very fascinating to confirm these numbers by direct field theoretic calculations based on the uniton or other non-perturbative solutions [25,36,7,27,37]. The first non-perturbative saddle is particularly interesting as it does not have any fluctuation part.
We observed important resurgence properties of the various functions but these relations showed only the tip of the iceberg. A more systematic extensive analysis should reveal the full trans-series and their resurgence structure, i.e. the formulation of the bridge equations. We believe that our research will spark new activities in this field. As the O(4) nonlinear sigma model is equivalent to the SU (2) principal chiral model generalizations for other N are possible in two directions. In this respect a double scaling limit [38] could simplify the analysis. Generalizations to other models including those in [35,30] would be also very interesting.

Outline
The paper is organized as follows: In the next section we summarize our setup for the O(4) model in a magnetic field and explain the perturbative calculation of the free energy in the renormalization group improved running coupling. We also explain how the same perturbative series can be obtained from the TBA equation. In section 3 we use the numerically determined high order perturbative coefficients and analyze the analytical structure of the Borel transform. We also calculate the inverse Borel transform and compare it to the TBA result, which we compute numerically. In Section 4 we perform the asymptotic analysis of the perturbative coefficients and quantitatively understand their leading singularity structure, including their nearest alien derivatives, which we investigate in a similar fashion. This enables us to fix the first few terms in the trans-series, but not enough to see any resurgence. We then perform a similar analysis for the density and the energy density as the functions of B, starting in section 4.5. Here we observe a very nice resurgence structure, which we translate to the original variables in section 5, and summarise in section 6. We use the various alien derivatives in section 7 to formulate the median resummation, which agrees with the TBA results and provides the first few terms in the trans-series. Some technical details are relegated to two Appendices.

The O(N ) model TBA and perturbative coefficients
The 2D O(N ) sigma models are exactly soluble, and provide useful testing grounds for phenomena appearing in QCD. Particles transform with respect to the fundamental representation of O(N ), and scatter on each other with an integrable elastic scattering matrix, which is exactly known. In a magnetic field, positively charged particles condense in the vacuum. These particles scatter diagonally on each other, and from the thermodynamic limit of the Bethe ansatz an integral (TBA) equation can be derived for their spectral densities. The systematic expansion of this integral equation provides a tool to calculate the perturbative coefficients at very high orders, as we explain in this section.

Perturbative definition of the model
We are interested in the Euclidean theory in the case when one of the conserved O(N ) charge, say Q 12 = (S 1 ∂ 0 S 2 − S 2 ∂ 0 S 1 )dx, is coupled to a magnetic field: Here λ is the bare coupling and the h-dependent terms are chosen such that the Hamiltonian is simply H(h) = H(0) − hQ 12 . In the perturbative calculations one introduces an infrared regulator −2ω 2 S 1 (which is put to zero at the end) in the Lagrangian to fix the ground-state to be S 1 = 1 and S i>1 = 0. By expressing S 1 = 1 − λ 2 (s 2 2 + · · · + s 2 N ) with λs i = S i for i > 1 we can perturbatively expand the free energy which is regulated in dimensional regularization. The first few terms in the perturbative expansion read as [39] where γ = Γ (1) + ln(4π). UV divergences can be get rid off by introducing the renormalized couplingg, writing λ 2 = (µe  . Thus the renormalized free energy can be expressed as The running of the coupling guaranties that the result is independent of the renormalization scheme. It is thus natural to introduce a renormalization group invariant scale Λ = By taking this scale in the M S scheme we can introduce the running coupling as The free energy density then has the expansion It is quite complicated to proceed with the perturbative calculations at higher orders. Fortunately the theory is integrable and we will be able to calculate the higher orders systematically. The quantities, however, on the integrable side are the density ρ and the ground-state energy density (ρ). They are related to the free energy by Legendre transformation: In order to express the ground-state energy density in terms of the density one can introduce a new running coupling such that In the following subsections we calculate the higher order terms from the TBA equation.

TBA calculation of the ground-state energy density
In the infrared description we start with the particle spectrum and their scattering matrices.
The O(N ) model has a particle multiplet which transform in the fundamental (vector) representation of the O(N ) group. This is a relativistic theory and the dispersion relation can be parametrized in terms of the rapidity as E(θ) = m cosh θ, p(θ) = m sinh θ. The scattering matrix S kl ij (θ), which depends on the difference of the rapidities, is non-diagonal and can be calculated exactly [10]. We are interested in the ground-state energy in a magnetic field, when the Hamiltonian is modified as H(h) = H(0) − hQ 12 . If h > m particles of type + (corresponding to the field S + = S 1 + iS 2 ) condense into the vacuum. In order to describe this condensate one can introduce a finite volume L and analyse momentum quantization via the Bethe Ansatz equation For h > m the groundstate is filled with particles of rapidities {θ j } ∈ I in an interval. Taking the L → ∞ thermodynamic limit we also have M → ∞ such that the density ρ = M L is finite. By introducing the rapidity density of states, χ(θ) such that Lχ(θ) dθ 2π is the number of states in the interval (θ, θ + dθ) the derivative of the BA equation has a thermodynamic limit: where B is a function of h and the kernel is where Ψ is the digamma function: Ψ(θ) = ∂ θ log Γ(θ). The density and energy are obtained simply as These equations depend on B, which can be related to the magnetic field via: h = ∂ ρ (ρ) = ∂ ∂B / ∂ρ ∂B , following from minimizing F = − hρ as a function of ρ. Thus either the magnetic field h, or the density ρ, or the parameter B can be used as the control parameter. Large magnetic fields correspond to large densities and large Bs, thus the perturbative expansion goes in 1/B. In the following we explain our understanding of Volin's perturbative solution of the TBA. Figure 1: Analytic continuation of the resolvent R(✓) with a short cut is displayed in the middle. Analytic continuation through the cut from the upper half plane results in R + (✓) on the left, while continuation from the lower half plane results in R (✓) displayed on the right. R ± (✓) agrees with R(✓) on the upper/lower half planes and they both have long cuts.

Perturbative expansion of the TBA
The basic idea is to solve the TBA equation in the bulk ✓ ⇠ 0 and in the edge regions ✓ ⇠ B perturbatively and to match the two expansions. The calculations are simpler for the resolvent which is analytic on the whole complex plane except on the interval ( B, B) where it has the jump R(✓ + i0) R(✓ i0) = 2⇡i (✓). The density is obtained from the residue of the resolvent at infinity, while the energy density through the Laplace transform This is related to the Fourier transform of (✓): Then, using the (✓) = ( ✓) symmetry In proceeding with the bulk solution we can introduce two functions R ± (✓) obtained by analytical continuations of R(✓) through the cuts, see Figure 1. R + (✓) is obtained by diving into the cut from above. Thus R + (✓) = R(✓) on the upper half plane and has two long cuts ( 1, B) and (B, 1). R (✓) is defined analogously by continuing from the lower half plane and it also has long cuts. The difference of the two functions is R + (✓) R (✓) = 2⇡i (✓). By adding the TBA equation at ✓ + i⇡/2 to that of at ✓ i⇡/2 the right hand sides cancel. By taking further into account the (z + 1) (z) = 1/z property of the digamma function the sum of the TBA equations for the resolvent takes the form In particular it means that f (✓) = R + (✓+ i⇡ 2 (2 1))+R (✓ i⇡ 2 (2 1)) satisfies f (✓) = f (✓ 2i⇡ ). By making an assumption for the asymptotic large B behaviour of the resolvent R(✓) = 8 Figure 1: Analytic continuation of the resolvent R(✓) with a short cut is displayed in the middle. Analytic continuation through the cut from the upper half plane results in R + (✓) on the left, while continuation from the lower half plane results in R (✓) displayed on the right. R ± (✓) agrees with R(✓) on the upper/lower half planes and they both have long cuts.

Perturbative expansion of the TBA
The basic idea is to solve the TBA equation in the bulk ✓ ⇠ 0 and in the edge regions ✓ ⇠ B perturbatively and to match the two expansions. The calculations are simpler for the resolvent which is analytic on the whole complex plane except on the interval ( B, B) where it has the jump R(✓ + i0) R(✓ i0) = 2⇡i (✓). The density is obtained from the residue of the resolvent at infinity, while the energy density through the Laplace transform This is related to the Fourier transform of (✓): Then, using the (✓) = ( ✓) symmetry In proceeding with the bulk solution we can introduce two functions R ± (✓) obtained by analytical continuations of R(✓) through the cuts, see Figure 1. R + (✓) is obtained by diving into the cut from above. Thus R + (✓) = R(✓) on the upper half plane and has two long cuts ( 1, B) and (B, 1). R (✓) is defined analogously by continuing from the lower half plane and it also has long cuts. The difference of the two functions is R + (✓) R (✓) = 2⇡i (✓). By adding the TBA equation at ✓ + i⇡/2 to that of at ✓ i⇡/2 the right hand sides cancel. By taking further into account the (z + 1) (z) = 1/z property of the digamma function the sum of the TBA equations for the resolvent takes the form In particular it means that f (✓) = R + (✓+ i⇡ 2 (2 1))+R (✓ i⇡ 2 (2 1)) satisfies f (✓) = f (✓ 2i⇡ ). By making an assumption for the asymptotic large B behaviour of the resolvent R(✓) = 8 Figure 1: Analytic continuation of the resolvent R(✓) with a short cut is displayed in the middle. Analytic continuation through the cut from the upper half plane results in R + (✓) on the left, while continuation from the lower half plane results in R (✓) displayed on the right. R ± (✓) agrees with R(✓) on the upper/lower half planes and they both have long cuts.

Perturbative expansion of the TBA
The basic idea is to solve the TBA equation in the bulk ✓ ⇠ 0 and in the edge regions ✓ ⇠ B perturbatively and to match the two expansions. The calculations are simpler for the resolvent which is analytic on the whole complex plane except on the interval ( B, B) where it has the jump R(✓ + i0) R(✓ i0) = 2⇡i (✓). The density is obtained from the residue of the resolvent at infinity, while the energy density through the Laplace transform This is related to the Fourier transform of (✓): Then, using the (✓) = ( ✓) symmetry In proceeding with the bulk solution we can introduce two functions R ± (✓) obtained by analytical continuations of R(✓) through the cuts, see Figure 1. R + (✓) is obtained by diving into the cut from above. Thus R + (✓) = R(✓) on the upper half plane and has two long cuts ( 1, B) and (B, 1). R (✓) is defined analogously by continuing from the lower half plane and it also has long cuts. The difference of the two functions is R + (✓) R (✓) = 2⇡i (✓). By adding the TBA equation at ✓ + i⇡/2 to that of at ✓ i⇡/2 the right hand sides cancel. By taking further into account the (z + 1) (z) = 1/z property of the digamma function the sum of the TBA equations for the resolvent takes the form In particular it means that f (✓) = R + (✓+ i⇡ 2 (2 1))+R (✓ i⇡ 2 (2 1)) satisfies f (✓) = f (✓ 2i⇡ ). By making an assumption for the asymptotic large B behaviour of the resolvent R(✓) = 8 Figure 1: Analytic continuation of the resolvent R(θ) with a short cut is displayed in the middle. Analytic continuation through the cut from the upper half plane results in R + (θ) on the left, while continuation from the lower half plane results in R − (θ) displayed on the right. R ± (θ) agrees with R(θ) on the upper/lower half planes and they both have long cuts.

Perturbative expansion of the TBA
The basic idea is to solve the TBA equation in the bulk θ ∼ 0 and in the edge regions θ ∼ B perturbatively and to match the two expansions. The calculations are simpler for the resolvent which is analytic on the whole complex plane except on the interval The density is obtained from the residue of the resolvent at infinity, while the energy density through the Laplace transform This is related to the Fourier transform of χ(θ): Then, using the In proceeding with the bulk solution we can introduce two functions R ± (θ) obtained by analytical continuations of R(θ) through the cuts, see Figure 1. R + (θ) is obtained by diving into the cut from above. Thus R + (θ) = R(θ) on the upper half plane and has two long cuts (−∞, −B) and (B, ∞). R − (θ) is defined analogously by continuing from the lower half plane and it also has long cuts. The difference of the two functions is By adding the TBA equation at θ + iπ/2 to that of at θ − iπ/2 the right hand sides cancel. By taking further into account the Ψ(z + 1) − Ψ(z) = 1/z property of the digamma function the sum of the TBA equations for the resolvent takes the form In particular it means that f (θ) = R + (θ+ iπ 2 (2∆−1))+R − (θ− iπ 2 (2∆−1)) satisfies f (θ) = f (θ− 2iπ∆). By making an assumption for the asymptotic large B behaviour of the resolvent R(θ) = one can show that f (θ) must be a constant. Its antisymmetry then implies that f (θ) = 0.
Although it is possible to proceed for generic O(N ) models we restrict our attention to the O(4) model where ∆ = 1/2, since formulas are much simpler there. In particular f (θ) = 0 implies that R − (θ) = −R + (θ), thus they both are analytic except for long cuts on the real line. By using the conformal mapping w = log x−1 x+1 the function R + (w) is analytical everywhere, thus must have the form k(x) = sinh w 2 F(cosh 2 w 2 ). In the original variable it implies the following expansion where the c n,m coefficients are numerical constants, but the overall constantÃ may depend on B.The density is obtained from the residue of the resolvent at infinity Let us focus on the edge region, where the Wiener-Hopf technique can be used. Once exponentially small corrections of the form e −B are neglected the Fourier transform of the TBA equation can be separated for an equation analytic on the upper and another analytic on the lower half plane [35]. This determines the Fourier transform of χ, which is related toR(s) by (17). Using this asymptotic behaviour leads to the following ansatẑ with constant coefficients Q n,m . By re-expanding R(θ) obtained in the bulk in the edge region, and performing the Laplace transform, the result can be compared toR(s). This givesÃ = A, and by matching the two parametrizations, all the unknown coefficients Q n,m and c n,m can be determined. This was used in [33,34] to calculate the perturbative coefficients for generic O(N ) models. The calculations in the O(4) model take a simpler form, which we present in appendix A. Having calculated the Q n,m and c n,m coefficients the energy density and the density can be written as and Eventually we are interested in the expansion of /ρ 2 in terms of the running coupling α, which is defined by the following relation Here we used the exact m/Λ M S value [9]. (25) can be solved by expanding B in powers of α.
The first few terms of this expansion are where z 3 = ζ(3) is the zeta function at 3. Using this B(α) relation we can express /ρ 2 in terms of α. As a result we could calculate this expansion analytically at 50th order, which goes beyond the result of [35].
We then switched to a high precision numerical implementation. This resulted in 2000 coefficients with 12000 digit precision (for ρ and ). The calculation ran on a PC for 5 days.
3 Perturbation theory and non-perturbative effects Above we described how to expand the TBA to produce results equivalent to those of ordinary perturbation theory, although allowing thousands of terms (as powers of α or 1/B), not just three.
But we can also solve the TBA equation (12) directly at finite coupling, without expanding, and find ρ and numerically. These direct results can be compared to those from summing up the perturbative series, for which we employ Borel resummation methods. In this section we set up and compare these two techniques, and observe that the most straightforward resummation omits exponentially small instanton-like contributions, of order e −8/α . These come multiplied by another power series in α, and, thanks to our very precise numerical results, we are able to fit a few terms.

Solving the integral equation
The TBA equation (12) is a linear integral equation, which can be solved analytically for large values of the B by the Wiener-Hopf method [9,35].
We are interested in the energy and the density ρ as the function of the running coupling α. Rather than work at fixed α, it is much easier to solve the equation at fixed values of B, and then with the help of the formula (25) recover points on the functions ρ(α) and (α) numerically.
The numerical method we applied was as follows. The unknown χ(θ) in the TBA equation (12) is expanded in even Tschebyshev-polynomials on [−B, B], up to order n c − 1 for some odd n c : T n (x) = cos(n arccos x).
Inserting this into (12) and evaluating at the zeros of the next polynomial T nc , namely θ k = B cos[(k− 1 2 )π/n c ] for k = 1, ..., n c , leads to a set of linear algebraic equations for the coefficients s j . Then formulas in (14) give and ρ.
We did this for B ∈ {1, 2, ..., 20}, which covers the range from non-perturbative to highly perturbative. The results can be found in table 1. These used n c = 171, which is sufficient for 37 digits of precision at B = 20, and more at lower values.  Table 1: Some of our results from the numerical solution of the TBA equation (12). At fixed B, and ρ are calculated from the numerical χ(θ), and˜ , α from them. We also display the numerical uncertainty of the result together with its deviation from the inverse Borel transformation based on the conformal mapping method. Note that the entries are truncated to 6 digits, while our numerical precision is much higher.

Borel-Padé resummation
Let us consider the renormalization group improved perturbative series for the normalized energy density as a function of the running coupling (8): The first few coefficients can be read off from (9) and a general method [33,34] to determine in principle all coefficients was explained in section 2. The application of this method allowed us to get 2000 coefficients of (28) with 12000 digits of precision, which makes it possible to implement the conformal mapping improved Borel resummation technique with very high precision and compare to the exact TBA data. This perturbative series has zero radius of convergence, because the coefficients χ n grow factorially at large n. The Borel transform is defined by dividing this out, and we choose the following conventions (removing also some powers of 2): The inverse of this is a Laplace transform: Applied term-by-term, this would trivially recover (28). But since the Borel series (29) converges inside the unit circle, it defines an analytic function. It is the integral of this function, analytically continued to infinity, which gives a resummation of˜ (α).
One commonly used representation of the analytic B(t) is given by Padé approximants. Knowing N coefficients c n , we can uniquely fix the coefficients of this rational function: From this, we observe that B(t) has singularities on both the positive and negative real axes, at |t| ≥ 1. They are drawn in figure 2, which shows that there is an isolated pole singularity at t = 1 and two dense sets of other poles which appear to condense into two cuts (starting at t = −1 and t = 2) as N is increased. The integration contour must then run either slightly above or below the positive real axis to avoid the singularities. We indicated this in (30) by endpoint ∞ ± i0, but in practice rotate the contour off the axis by some small angle. We write˜ (±) (α) for the results with either choice.
Both˜ (+) and˜ (−) have an imaginary part, thus they cannot give the correct physical result. This non-perturbative ambiguity is defined as |u|=0.5 |u|=0.7 |u|=0.9 Figure 2: Left, positions of the poles of a 100th-order Padé approximant of the Borel transform ofˆ , in the complex t plane. These accumulate along cuts t ≤ −1 and t ≥ 2, plus an isolated pole at t = 1. Right, lines of constant |u| for the conformal mapping (34), under which the disk |u| < 1 covers the whole t plane, minus these cuts.
From our Padé analysis we conclude that the singularities (along the positive real axis) are: and a logarithmic cut starting at t = 2. The leading ambiguity is coming from the pole term: (+) and˜ (−) share a common real part, in terms of which we try to approximate the physical result, and study its deviation from that as a function of the running coupling.

Conformal mapping
For the actual calculation of (30), instead of the Padé approximants, we used the conformal mapping method. This method is often applied even to real physical problems [2] to approximate physical quantities better from available perturbative data. In our case we can make a high precision comparison between the perturbative and the exact results getting a deeper insight into the mathematical structure of the deviation.
The Borel-transform (29) as it stands is not very useful, since it is still a Taylor-series, which is convergent inside the unit circle only and thus not applicable in the whole range of the integration contour. As explained above, the singularities of the analytically extended B(t) lay on the real axis for |t| > 1. This makes it possible to perform a conformal transformation of the form: which defines a convergent Taylor-series in the unit circle of the u-plane, such that the singularities of B(t) and the integration contours are mapped to the boundary and into the interior of the unit circle, respectively. For the actual computations we consider the inverse transformation: The coefficients b n of the Borel-transform on the u-plane can be determined from those of the original t-plane by matching the small u-expansion of the two different representations: In this way, we could determine the first 2000 coefficients b n . The improved (approximate) inverse Borel-transform we used in our calculations is

Fitting the non-perturbative corrections
We calculated˜ (+) (α) as explained above and compared it to the exact (TBA) results (which are of course real). We have numerically fitted the correction terms using the ansatz Since we know the exact value of the residue of the pole singularity (see next section) our fit provides the NLO term for the imaginary part of the ambiguity. As our ansatz for the real part shows we found that the leading non-perturbative ambiguity for the real part is extremely small, it is of fourth order in the non-perturbative expansion parameter e −2/α . We have studied the stability of the first few fit coefficients by taking higher and higher order fits to the correction functions. This way we obtained the following quite precise results for c and d coefficients:

Resurgence
Our data are precise enough to see that, in addition to the simple Borel resummation, extra non-perturbative corrections are clearly needed. So far we have only obtained these correction terms by curve fitting, but our goal is to recover them from the perturbative coefficients found in section 2. To do this we need resurgence theory, which this section sets up.
The following sections explore the pattern uncovered, translate back from an expansion in 1/B to an expansion in α, and ultimately use these results to recover the non-perturbative correction terms.

Asymptotic coefficients and cuts in the Borel plane
The non-perturbative corrections are intimately related to the singularity structure of the Borel transform B(t). To study these singularities we use the well-known relation between the cuts (and poles) of B(t) and the asymptotic behaviour of its expansion coefficients.
The series {c n } contains two regularly behaving sub-series: In order to see that we have devided by the correct factorial growth in (29) we investigated various Richardson transform of the a n coefficients. Demonstrative results are presented on Figure (3). By subtracting the leading order behaviour a similar analysis shows that for asymptotically large n the coefficients behave as Our method to calculate the asymptotic coefficientsÃ k ,B k is described in appendix B.
Having computed the asymptotic coefficients we can write down the singular part of the Borel transform, which is of the form where Around t = 1 we can expand the coefficient of the log term in powers of t − 1 and it becomes Similarly, around t = −1 we have the expansion The new expansion coefficients p m , q m are linear combinations of the asymptotic coefficients. The first few coefficients are and q 0 = 2B 1 , (52)

Results for the asymptotic coefficients
We observed from our numerical results that all p m coefficients m = 0, 1, . . . vanish. We established this fact to more than a hundred decimal digits for the first few p m . Of course the higher coefficients are less and less precise, but since this observation is in agreement with our previous finding, namely that apart from a pole at t = 1 there is no singularity up to the cut starting at t = 2, from now on we take it for granted that all p m exactly vanish. We will also assume that the precision of the non-vanishing coefficients q m is similar to that precision by which the corresponding p m vanish, that is about 147 digits for q 0 , 142 digits for q 1 , and gradually decreasing until it remains only 6 digits for q 90 .
Our next observation is that The numbers in parenthesis indicate that the above relations are satisfied to 153 (or 152) digits by our estimated results. Clearly we can safely assume that (53) are satisfied exactly.
With the help of the zeta function webpage EZ-Face-CECM * we were able to find exact expressions for the first few q m coefficients. We found q 0 = 0, (54) Here the notation z k = ζ(k) = Zeta[k] is used.

Alien derivatives
The notion of alien derivative is a concise and elegant way to characterize the logarithmic cut (and pole) structure of the Borel transform. We refer to [8] for details and definition here we merely summarize the connection between asymptotic coefficients and alien derivatives. For a formal asymptotic expansion we introduce the coefficients of its Borel transform and apply the asymptotic analysis explained in Eqs. (44)-(46) and Eqs.
(47)-(52). We calculate the asymptotic coefficientsÃ m ,B m and the expansion coefficients p m , q m . The alien derivative at t = 1 is then given by In this language the alien derivative of the energy density at 1 is a constant and the alien derivative at -1 is characterized by the coefficients (54), which look like perturbative expansion coefficients around some saddle point, but we were unable to find any obvious sign of resurgence in this structure.

Results for singularity around t = 2
We have established that there is only a pole singularity at t = 1 and no cut is starting there. Furthermore the residue of the pole is exactly known. If we remove this exactly known pole at t = 1 from the Borel transform there remains no singularity between −1 and 2. This allows us to re-expand it around t = 1/2. The convergence radius of this new series is 3/2. We have applied our asymptotic analysis described above † to the coefficients of this new series. We find ‡ for f = (2/πα)˜ and expansion parameter 1/x = α/2 where the first few coefficients arep  4,5,6 are also known analytically. The q m coefficients associated to the singularity at t = −1 are already known (more precisely) from our previous work and can be used to estimate the precision of the correspondingp (2) m coefficient. This method suggests that the precision of p (2) m between m = 7 and m = 51 gradually decreases from 88 digits to 6 digits. Again, the coefficients (60) suggest that they come from perturbative expansion around some saddle point but we could not recognize any resurgence here.
The cut characterized by D 2 f is responsible for the NLO corrections in (39). Having obtained the overall factor and the first few coefficients exactly, we can now give exact formulas † One has to rescale the variables by 3/2 appropriately to absorb the effect of an increased convergence radius.
‡ For the definition of the alien derivatives Dω see section 5.
for the expansion coefficients c i , i = 1, . . . , 5: The extremely good agreement between the exact coefficients and our previous fit results (42) (including the error estimates) makes us confident that our methods are consistent.

Resurgence properties of the basic functionsˆ andρ
In this subsection we reconsider the whole resurgence analysis going back to the more elementary building blocksρ(B) andˆ (B). Moreover, we will use the original B parameter to define the new expansion parameter 1/z, where z = 2B. Using Volin's method, we calculated the expansion coefficients in the asymptotic serieŝ where c n = 2 n+1 u n+1 n! .
The first two expansion coefficients are The symbol a is again shorthand for the numerical value a = ln 2. We have calculated the first 16 coefficients analytically and observed that the coefficients depend polynomially on a and that the u p coefficient has order p concerning its zeta-function dependence. The highest power of a occurring in u p is a p and more generally this coefficient is a linear combination of terms of the form f k a k , where f k is a zeta-function combination of order not larger than p − k. We calculated the first 2000 coefficients numerically (putting a = ln 2 numerically) with several thousand digits precision. The asymptotic analysis of the serieŝ is completely analogous. The first two coefficients are and by calculating the first 16 coefficients analytically we can establish that the zeta-function dependence and polynomial a-dependence of ξ p+1 is similar to that of u p . Remarkably, the ln 2-dependence completely cancels from the densityε, if expressed in terms of the running coupling α.
In this subsection we will use the language of alien derivatives [8] and will extensively use two of its basic properties: • The alien derivative ∆ ω is a proper derivative, it is linear and satisfies the Leibniz and chain rules.
We can go to the Borel plane and study the singularity structure of the Borel transform and similarly for B (t). We use again the method of asymptotic analysis described in subsection 4.1.
Having computed the asymptotic parameters for bothρ andˆ very precisely, we were able to recognize the functions appearing in the ∆ −1 alien derivatives. The ∆ ±1 alien derivatives of the two basic functions are given below.
Using these basic derivatives, we can derive further relations for the basic functions and for the combinations F is useful because it is related to the energy density and G has nice properties as we will see below. We use prime to denote derivative w.r.t. z. Using (68) it is easy to verify It is also easy to see that and for later purposes we note that

Further alien derivatives
Since the Borel transform of the function 1/ˆ has only a pole singularity at t = −1, after removing this exactly known pole no singularity remains between −2 and 1. Therefore we re-expanded the corresponding subtracted Borel transform around t = −1/2 and performed a (rescaled by 3/2) asymptotic analysis of the coefficients. We found that implying ∆ −2ˆ = 0. The function G was constructed so that hence its expansion around t = 0 has convergence radius 2. Applying the (rescaled by 2) asymptotic analysis we find This conclusion is based on the observation thatB 0 = 0, q 0 = 0 to 182 and 177 digits respectively, and q i = 0 (i = 1, . . . , 76) to (172, . . . , 6) digits. At t = 2 we find to 182, 177 digits respectively, and The first twop m coefficients arep 0 = 9 2 − 4a, (78) p 2,3,4 are also exactly known. The precision of the exactly knownp m , m = 0, . . . , 4 is (177, 171, 167, 162, 158) respectively and we calculated thep m coefficients up to m = 76 numerically. Between m = 1 and m = 76 the precision changes from 171 to 6 digits. The next observation is thatp m is bounded for m → ∞ implying ∆ ±1 H = 0 and that the Borel transform of H also has convergence radius 2. We note that ∆ −2 G = 0 and ∆ −2ˆ = 0 together imply that also We performed a (rescaled) asymptotic analysis of the coefficients of the Borel transform of H. We found that where the first two coefficients arẽ Using similar tricks we established that The first 5 coefficients r 1 , . . . , r 5 are known analytically: and similarly e 1 , . . . , e 5 , while the rest up to n = 50 numerically. The next trick is to consider the function We found that {h m } is bounded for m → ∞ implying ∆ ±1 η = 0. Thus Taking into account that also ∆ ±1 H = 0, we derive the relations There is also some evidence (a few digits) of the vanishing of ∆ −2 η. This, combined with We have also analysed the t = 2 singularity. Making the definitions the result for the leading expansion terms is 5 Translation back to the running coupling α It is natural to use the variable z = 2B both for the TBA calculation and for calculation of the coefficients in Volin's expansion. In perturbation theory on the other hand, the natural variable is the running coupling α. We thus introduce the expansion parameter The free energy density we were studying originally is given bỹ and in subsections 4.1 and 4.4 we studied the resurgence properties of Rewriting (25), the relation can be (perturbatively) calculated from with the result starting as In the following, we will denote (as before) by ∆ ω the alien derivative of functions expanded in 1/z, but use D ω for alien derivatives of functions expanded in 1/x. For a composite function C = γ • z, C(x) = γ(z(x)) the D ω alien derivative is given by the formula [40] D ω γ(z(x)) = e −ω(z(x)−x) (∆ ω γ)(z(x)) + γ (z(x))(D ω z)(x). (96) Applying (96) to γ =ρ and combining it with the alien derivative of (94) gives where dot denotes d/dx. We are interested in the alien derivatives of f (x). Taking we obtain For ω = ±1 we get as in subsection 4.1, and Combining the above two results we also get We have checked that D −1 f agrees with the result obtained in subsection 4.1 numerically and, using (100), we can verify that So far we have established that Using (98) D 1 F = D −2 F = 0 immediately follows and after a very long calculation we have verified that For the relevant exponantially small correction to TBA, we also calculated the expansion of  Table 2: Resurgence structure of first and second generation functions, from (108) to (115).
Here "res" means resurgence to functions encountered above, "const" means a number, and those marked "+" are also calculable and nonzero.

Resurgence patterns
Here we collect the results found in the last two sections, to show clearly the pattern of resurgence. Let us call the functionsρ,ˆ , G, and f elements of the first generation. Of course, only the first two of these are really fundamental, but we include G here because it has particularly simple resurgence structure and f because it is simply related to the free energy and we were originally interested in it. Similarly, the elements of the second generation are R, E, H, and F. Finally, the third generation consists ofR,Ẽ,H, andF. The resurgence pattern for the elements of the first generation is For the second generation it is This patten is illustrated in table 2. In this table the symbol "const" means a numerical constant, and the symbol "res" means resurgent, in the sense that the pertinent alien derivative can be written in terms of functions of the same or earlier generation. Note that the members of the families (ˆ , E,Ẽ), (G, H,H), and (f, F,F) resurge within the family. Finally, the meaning of the symbol "+" is that the pertinent alien derivative can be calculated by re-expanding the Borel transform around −1/2. But we do not know how to calculate ∆ 3 alien derivatives.
In addition to this a pattern of generations, there is another pattern, which is most easily seen in different variables. Consider the following representatives of the first generation: These are exchanged by ∆ ±1 as follows: Acting with ∆ 2 takes us to the second generation, and we observe that exactly the same pattern of ∆ ±1 holds: where the basis functions are This pattern is drawn in Figure 4, starting with the first generation at the top. Acting with ∆ 2 again will take us to the third generation, containing alsoR andẼ. Here we have not derived the actions of ∆ ±1 above, but the seemingly obvious conjecture is that the same pattern persists, which would allow us to solve for these: However, this pattern does not capture everything. We have expanded φ − on the Borel plane around t = −3/2 and performed an asymptotic analysis, which showed that ∆ −3 φ − = 0. This implies that the full resurgence pattern must be more complicated.

Median resummation and the cancellation of ambiguities
Having calculated the relevant alien derivatives of f = 2 απ˜ we are in the position to propose an ambiguity free resummation of the perturbative series. Clearly the lateral Borel resummations The ambiguity free median resummation involves the square root of the Stokes automorphism and takes the form [42,41] S med (f ) = S − (S Let us recall the alien derivatives we have calculated Since D 1 f is a constant all higher alien derivatives are vanishing D k D 1 f = 0. We also calculated the first few terms of the alien derivative of F : Using these results the median resummation takes the form: (126) § Observe that our definition of the alien derivative (57) is such that it is the logarithm of the inverse of the Stokes autormorphism. This is the opposite which is used in the literature [6,41].  Table 1) to median resummation (129) assuming D 1 D 3 f = 0. The lateral Borel resummation Re(S + (f )) has been subtracted from both, and the exponential prefactor e −8/α has been divided off. The points shown agree to 3 digits.
The expression in the bracket provides the ambiguity free trans-series of the free energy. We did not manage to calculate D 3 f and higher derivatives, but already these terms can be compared to the numerically obtained TBA results. First we can check that the result is real. Indeed, single alien derivatives are purely imaginary, and the S + (D 1 f ) and S + (D 2 f ) terms cancel the exponentially small imaginary c-terms (42) coming from S + (f ). This can also be seen from the definition of the Stokes automorphism and by noting that the imaginary part of S + (f ) is half of the difference. In calculating the leading real exponential contribution we point out that the term cancelling the imaginary contribution S + (D 2 f ) also has a real part, which can be read off from This combines with the direct D 2 2 f term giving We can compare these with the coefficients we determined previously − 2d 1 π = 0.58607 32 e 4 = 0.58610 and observe complete agreement within the available precision, which actually indicates that D 1 D 3 f = 0. This comparison is also shown in Figure 5.
The agreement between the TBA results and the median resummation gives a strong evidence of the first few terms of the trans-series. The form of the full trans-series is expected to be where χ (0) n = χ n are the perturbative coefficients, χ n = − 16i e δ n,1 are related to D 1 f , while χ (2) n to D 2 f .
We have also investigated the relation between the numerical solution of the TBA equation and the median resummation of the perturbative series directly for the basic building blockŝ ρ andˆ in terms of the original TBA variable B.
(137) In the Laplace transform we have to expand in s. The explicitly known terms contribute to m = 0 V Q (n, 0) = Γ( 1 2 ) Γ(n + 1) while the Q coefficients to m > 0: Q r,m−r−1 G n+r+1 ; G n = Γ( 1 2 ) Γ(n + 1) The coefficients, c n,m and Q n,m can be determined by demanding V Q (n, m) = V c (n, m). In solving these equations it is very natural to proceed in m and solve all the coefficients iteratively in terms of smaller m values. For m = 0 we can start the iteration as c n,0 = V Q (n, 0)/F n,0 . For m = 1 we can write a separate equation for n = −1 giving Q 0,0 = c 0,0 F −1,1 /G 0 and another for n > −1 resulting in c n,1 = (Q 0,0 G n+1 − c n+1,0 F n,1 )/F n,0 . Now let us assume that we have already used the equations up to m − 1 to determine Q i,j for i + j < m − 1 and c i,j for j < m. We then use the equations for m to determine Q i,j for i + j = m − 1 and c i,j for j = m. In doing so we start with the equation for n = −m and proceed one by one to n = −m + 1, −m + 2, . . . , −1 to obtain all the Qs: We implemented this iterative calculation in Mathematica. In order to speed up the calculation we use the Γ-function identity thus we have to expand the Γ-function only around 1. We then used the functional relation Γ(x)Γ(1 − x) = π sin(πx) together with the following recursive expression a k x k ; a n = na 1 a n − a 2 a n−1 + n k=2 (−1) k ζ(k)a n−k (143) to speed up the calculation.

B Asymptotic coefficients
In this appendix we summarize our asymptotic analysis.
Let us now define This can be calculated, starting from ξ (0) = ξ, using the recursive formula For the asymptotic coefficients we get with X (j) Most importantly, After j steps the j th coefficient of the original series is promoted to leading term.

B.2 Calculation of the constant term
Let us study the series of transformations After the first step, the 1/n term is eliminated and the asymptotic expansion becomes (1 − k)X k (n − 1)(· · · )(n − k) .
In the subsequent transformations, the O(1/n 2 ), O(1/n 3 ), etc. terms are eliminated step by step, and after s steps the asymptotic form is i.e.
The s th approximant of the leading (constant) term X 0 is thus where x [s] max is the last available member of the series (max is common to all s). Applying (159) for the transformed series ξ (m) we get an approximation for X m in the form

C Alien derivative of function composition
In this appendix, we prove the relation which we used for the alien derivative of composite resurgent function ψ[φ(x)], where both ψ and φ are resurgent functions of x. This relation was proposed in [40]. Here, we present its proof. We start by the following form of the composite function which is suggested by the theorem 0.3.2. of [40] and the following form is assumed for the resurgent function φ(x) φ(x) = x (1 + ε(x)) .
The alien derivative of the composite function (161), has two terms where, we used the Leibniz rule for ∆ ω . In the first term of (163), we should commute the alien derivative with the n-th ordinary derivative. We know that the ordinary and alien derivative do not commute, Some algebraic calculation shows that Using this result, we change the order of ∆ w and d n dx n in the first term of (163), where, we changed the order of the summation on n and in the second line. In the third line, we changed the variable from n to k = n − which runs from zero to infinity. In the last line, we used (161) to write the closed form of the second bracket of the third line, i.e. (∆ ω ψ) [φ(x)].
The second term of (163) is Combining the two parts of (163) gives the desired result for the alien derivative of the composite function