Treating Divergent Perturbation Theory: Lessons from Exactly Solvable 2D Models at Large $N$

We consider the operator product expansion (OPE) of correlation functions in the supersymmetric $O(N)$ non-linear sigma model at sub-leading order in the large $N$ limit in order to study the cancellation between ambiguities coming from infrared renormalons and those coming from various operators in the OPE. As has been discussed in the context of supersymmetric Yang-Mills theory in four dimensions, supersymmetry presents a challenge to this cancellation. In a bid to solve this problem we consider $O(N)$ as a toy model. A background field method inspired by Polyakov's treatment of the renormalization of the bosonic $O(N)$ model is used to identify explicit operators in the OPE of the two-point functions of bosonic and fermionic fields in the model. In order to identify the coefficient functions in the OPE, the exact two-point functions at sub-leading order in large $N$ are expanded in powers of the natural infrared length scale. The ambiguities arising from renormalons in the coefficient functions and vacuum expectation values of operators in the OPE are shown to cancel to all orders. The question of supersymmetric Yang-Mills theory without matter remains open.


Introduction
In this paper we study the supersymmetric (SUSY) two-dimensional O(N ) model in the 't Hooft limit [1,2,3,4]. This model is asymptotically free and, hence, the running coupling constant explodes at momenta of the order of the dynamical scale Λ. Perturbative expansions become meaningless. However, this model can be exactly solved in the leading and next-to leading orders in 1/N . The exact solution teaches us what happens to the divergent perturbation theory and how it conspires with non-perturbative terms in the operator product expansion (OPE). In the past some of the aspects of this problem were considered in non-supersymmetric version [5,6,7]. Supersymmetry introduces additional challenges due to vanishing of certain operators in OPE crucial in canceling perturbative divergences [8].
In Ref. [8] this issue was analyzed in the context of four-dimensional super-Yang-Mills theory (SYM) (see also the second paper in [9]). As is well known from multiple previous studies in (non-supersymmetric) QCD the leading infrared (IR) renormalon conspires with the gluon condensate operator G µν G µν (see e.g. [10]) in the operator product expansion (OPE). However, in SYM this mechanism of ambiguity cancellation fails: indeed, the vacuum expectation value (VEV) G µν G µν = 0 because G µν G µν is proportional to the trace of the energy-momentum tensor (up to an operator proportional to equation of motion). A way out proposed in [8] might be relevant in theories with matter. The puzzle apparent in the simplest supersymmetric Yang-Mills theory, N = 1 SYM without matter, remained unsolved. The correspondence between diagrammatic renormalon arguments and the OPE in this theory is not yet established.
In a bid to advance in this problem we turn to a simpler model in two dimensions which, however, preserves features of 4D SYM, namely the N = 1 O(N ) model at large N . Since the model is exactly solvable we can analyze the exact solution, represent it in the form of OPE, and explicitly verify that the renormalon ambiguities in the coefficient functions do conspire with those in the matrix elements of the OPE operators. In particular, we will consider the OPE of the two-point correlation functions of the bosonic and fermionic fields and demonstrate the cancellation of ambiguities to all orders using a form of the OPE similar to that considered in [7]. This form of the OPE somewhat obscures contributions of individual operator VEVs, but we will also verify the cancellation for the lowest few dimensions of operators in the OPE explicitly. Of course, we can be certain in these cancellations even before isolating renormalons since the exact solution is well defined and has no ambiguities.
The basic problem in 4D N = 1 SYM is that the trace of the energy momentum tensor is The VEV of the trace anomaly in a supersymmetric theory should generically vanish since it is a higher component of an anomaly supermultiplet. Furthermore, following the arguments in [9] it was argued that both the gluon and gluino terms in the trace individually vanish. As a result, there is an absence of low-dimension operators with non-vanishing VEVs which may appear in the OPE to cancel with low-dimension renormalon poles.
A certain parallel between 4D SYM and the 2D supersymmetric O(N ) model exists. In the non-supersymmetric O(N ) model the dimension-2 trace of the energy momentum tensor cancels with the lowest renormalon in the identity coefficient function. Using the equations of motion this trace may equivalently be expressed in terms of the Lagrange multiplier field enforcing the constraint of the O(N ) model. However, passing to the supersymmetric O(N ) model we observe that the corresponding fields in the supersymmetric O(N ) model are also higher components of superfields and so, as in 4D SYM, they must also have vanishing vacuum expectation values and can play no role in conspiracy with the renormalons in the OPE coefficients. So it might seem that there would be no lower-dimension operators to cancel potential renormalon poles.
At this point 4D SYM and 2D O(N ) diverge. The question of the existence of lowerdimension renormalon poles has been investigated in recent works considering the SUSY O(N ) model with a chemical potential [11,12]. The vacuum energy as a function of charge density was expanded as a power series in a parameter which is related to the coupling constant and the conclusion was that a lower-dimension renormalon pole does exist in this particular power series. So we wish to resolve the issue of presence of lower-dimension renormalons in the SUSY O(N ) model despite the hasty argument given against them in the paragraph above.
The question of infrared renormalons was also investigated in many other closely related 2D sigma models such as the principal chiral model [13] and the CP N model [14], where the power series considered was that of the VEV of "composite photon" condensate F µν F µν . Indeed, a renormalon ambiguity of the appropriate dimension in this VEV was seen. More interestingly when the model was considered on compactified space the location of the pole in the Borel plane was shifted compared to the uncompactified space, which raises some questions for the program of describing renormalons in terms semi-classical bions [15,16], which has been a major source of motivation for reconsidering the question of infrared renormalons in the last decade.
It is worth emphasizing that in non-solvable models such as QCD the OPE construction requires introduction of the normalization point µ such that µ Λ. Below µ the very notion of perturbation theory becomes senseless because at µ > ∼ Λ the effective Lagrangian cannot be formulated in terms of gluon and quark operators; the latter must be replaced by hadronic operators, see Section 7 for more details. There are no renormalons if µ is fixed. Renormalons in the sense usually discussed (and considered below) emerge only in the limit µ → 0 which can only be achieved either at weak coupling or in the exactly solvable theories. Thus, the bona fide verification of conspiracy between renormalons and VEVs of relevant operators in OPE can only be carried out in the abovementioned cases. This is our main goal in the present paper. The best one can do in the case µ Λ is to check that the artificial parameter µ drops out of OPE [6].

Outline of calculation
The paper is structured as follows. In Section 2, both the bosonic and supersymmetric O(N ) models are introduced and equations of motions which relate the various fields in the model are presented. It will be clear that indeed the direct analogues of the operators in the bosonic O(N ) model which cancel the renormalon ambiguities in the identity coefficient function can no longer fulfill that role due to the restrictions of supersymmetry. However there will be a host of other dimension-2 operators which may pick up the slack, so lower dimension renormalons can not be ruled out. The presence of these operators to some extent already resolves the difficulty stated above, and the bulk of this paper is about seeing in detail how the cancellation of ambiguities between operators and coefficient functions occurs in the operator product expansion.
To do this we will consider the correlation function of both the bosonic field n and its superpartner ψ in the SUSY O(N ) model. Schematically the correlation function of two general fields f in momentum space may be written as an operator product expansion where here j is an index referring to the engineering dimension of an operator and C j is its associated coefficient function which contains all dependence on the external momentum p.
In the scheme we are using here, the coefficient function will also contain all dependence on simple powers of the parameter g as well, and it can be calculated in ordinary perturbation theory in g. The vacuum expectation value of the operator O j will instead contain all dependence on the infrared parameter m obtained through dimensional transmutation as some power of exp − 1 g 2 . These coefficient functions C j expressed as power series in g will have some infrared renormalons and the goal of this paper is to show how they cancel in the full OPE.
The first step in showing this is to identify what the operators O j are, and to get a method that in principle could calculate the full OPE via perturbation theory. This is done in Sections 3 and 4. It is done via a background field method based on the scheme used by Polyakov to find the renormalization of the bosonic O(N ) model [17]. This method is used first on the bosonic O(N ) model in Section 3 in order to explain some of its features in a somewhat simpler context. The background field method applied the SUSY O(N ) model will be much more involved, and it is treated in Section 4. Along the way, we will provide a demonstration of the one-loop beta function of the SUSY O(N ) model in the Polyakov scheme, which was after all the purpose it was originally developed for.
Part of the reason that finding the OPE of the O(N ) model is tractable is because many simplifications occur at lowest order of the large N expansion. In particular, if we indicate the order in the large N expansion by a superscript, we will see in Sections 3 and 4 that the zeroth order coefficient functions C (0) j are particularly easy to find using this background field method, even for operators with vanishing VEV at lowest order. Such an operator would not contribute to lowest order correlation function itself, which looks like a free massive propagator in these models. However these operators will have nonvanishing VEVs O j (1) and will contribute to the correlation function expanded to O (N −1 ).
These leading order coefficient functions C (0) j and subleading VEVs O j (1) will be calculated in Section 5 for the first few terms in the OPEs of the n and ψ correlation functions.
Then the idea is that the renormalon ambiguities coming from the sub-eading coefficient functions C (1) j will in fact be canceled by ambiguities coming from the subleading VEVs O j (1) . This is an idea due to David [18,19], and it will also be briefly reviewed in Section 5 and the relevant operator ambiguities calculated.
The final step in our calculation is to find the renormalon ambiguities in the coefficient functions C (1) j themselves and verify that they cancel with the aforementioned operator ambiguities, and this will be discussed in Section 6 and its associated Appendix B. The coefficient functions are found by considering the exact correlation functions to subleading order in the large N expansion and calculating the asymptotic expansion in terms of the infrared scale m. This is similar to the method of considering the OPE for the self-energy in the bosonic O(N ) model in [7]. The various terms of this asymptotic expansion will have ambiguities, but the cancellation of ambiguities will be demonstrated to all orders in j. To a certain extent this asymptotic expansion will obscure the contributions of individual operators in the OPE, but certain terms in the asymptotic are clearly identifiable as coefficient functions due to dependence on the running coupling g(p), and the renormalon ambiguities from these terms will indeed cancel with the specific operator ambiguities calculated in Section 5.
The coefficient functions C (1) j and their renormalon ambiguities could of course be calculated in a way which is more analogous to the way they are calculated in 4d gauge theories. We could consider ordinary perturbation theory in g using the background field method of Section 3 and 4. At subleading order in large N , diagrams containing arbitrary numbers of "bubbles" of the perturbative field ϕ would lead to infrared renormalons. As a supplement, in Appendix A this arguably more direct method is considered for the identity coefficient function C (1) 0 in the bosonic O(N ) model. It is shown to lead to the same result as an asymptotic method [21] which is similar in principle to that considered in Section 6.
Finally section 7 presents a general discussion of our results, putting them in context, and explaining why exactly solvable models exhibit the conspiracy between renormalons and OPE operators.

Conventions and equations of motion
Here we will briefly summarize our conventions for the bosonic and supersymmetric forms of the O(N ) model, and in particular discuss the equations of motions relating the ordinary sigma model fields and the Lagrange multiplier fields, and their impact on VEVs at lowest order in the large N expansion. This section mostly follows the notation of [20] translated to Euclidean space.

Bosonic O(N ) model
The action for the bosonic O(N ) model involves the N component field n, and the Lagrange multiplier field λ. The constraint is taken to be n 2 = 1, with the coupling g −2 an overall factor multiplying the action. Then, We are free to add an arbitrary mass term for n since due to the presence of the constraint this is equivalent to an overall constant in the action. The mass m will be chosen in order to set the saddle point value of λ to zero after n is integrated out. This is equivalent to requiring the vanishing of all terms linear in λ, which implies Here M is a UV cutoff, and this condition gives the renormalization group flow of g 2 at lowest order in large N . Note that in the scheme used here, m 2 will always refer to this saddle point value given exactly by (2.2), and so it will start to differ from the physical mass of the excitations at subleading order in large N . The physical mass may be found directly from a pole of the correlation function (5.3) and is shown in Eq. (6.2). The emergence of the 1/N correction in (6.2) is due to the fact that the first coefficient in the β function for g 2 in the O(N ) model is N − 2 rather than N . It is more convenient to express all equations below in terms of m 2 defined in (2.2) rather than in terms of m 2 phys . It is certainly possible to re-express them in terms of m 2 phys , but then we will have to deal with bulky formulas. A propagator is generated for λ once n is integrated out, where J is the following frequently appearing one-loop integral, We will have occasion to make use of the following limiting behavior of J(p 2 , m 2 ), (2.6) Returning to the action (2.1), the equations of motion for n are −∂ 2 n + m 2 n − λn = 0. (2.7) The operator (∂n) 2 can be related to λ by dotting this equation by n and using the constraint, which leads to the zeroth order (saddle point) VEV However, note in passing that this manipulation is not an entirely innocent operation at higher orders in large N . The constraint is an equation of motion associated to the field λ itself, which appears at the same point as n in the combination λn. So dotting the equation of motion by n will involve something like a contact term. And in fact, at least in the one-dimensional O(N ) model [22], this λn term will also give a finite contribution to the VEV of (∂n) 2 at subleading order in large N . Analogous statements hold true for the SUSY O(N ) model as well, but since this will not affect the cancellation of renormalon ambiguities discussed in Section 6 we ignore this subtlety in what follows.

SUSY O(N ) model
In the SUSY O(N ) model the fields n and λ are extended to superfields Φ and Γ, where ψ, u are Majorana fermions, and λ F is labeled with a subscript since it will be redefined shortly.
The superderivatives are defined as and these can be used to form a kinetic term for the action, where a saddle point VEV from Γ has been explicitly separated out, as in the treatment of the bosonic case above. The action for F is quadratic, with equation of motion After integrating out F , the ordinary field Lagrangian (defined for later convenience without the conventional factor of 1/2) becomes 12) and the full action including the Lagrange multiplier fields becomes S = 1 2g 2 d 2 x (∂n) 2 + m 2 n 2 − iψ ∂ψ + mψψ + σ 2 n 2 + σ ψ ψ + 2m − λ n 2 − 1 + 2ūψ · n, (2.13) where λ F has been redefined as 14) The condition that λ and σ vanish at the saddle point, which relates g and m to the regularization, is given by the same condition (2.2) as in the bosonic case. The λ propagator is given by the same expression as the bosonic case (2.3) as well, which is repeated here along with the σ and u propagators. .

Equations of motion and VEVs in the SUSY O(N ) model
Now the equations of motion can be used to relate the fields of Φ and Γ. The action for σ is quadratic and leads to equation of motion ψψ = −2 (m + σ) . (2.16) The equation of motion for ψ is −i ∂ψ + (m + σ) ψ + un = 0.
Dotting by n, and using the constraints leads to Dotting byψ, and ignoring the potential problematic contact terms as discussed for the bosonic case, The similar equation of motion for n is So in total, the ordinary field Lagrangian (2.12) is just equal to λ F , Equation (2.19) in particular is a direct generalization of the equation in the bosonic case where of course −m 2 may be considered the saddle point value of a redefined λ if we wish. As discussed in Section 1, the SUSY O(N ) model would appear to immediately face a problem since both λ F and L, being F -terms of superfields, must have vanishing VEVs to all orders in N . But λ F and L are the analogues of the dimension-2 operators that cancel the renormalon in the non-supersymmetric (bosonic) case. So there has been some discussion whether a dimension-2 renormalon can appear at all in the SUSY O(N ) model (and the closely related CP N model) [12]. But of course it should be clear now that a number of other dimension-2 operators exist that can take the place of L and λ F in canceling any low-dimension renormalons, namely the individual terms of L in (2.12). Using the equations of motion (2.16), (2.18) (2.19), the vanishing of the zeroth order VEV of σ in conjunction with factorization at leading order of the large N expansion, we find zeroth order VEVs The last expression follows from (2.16).

Bosonic O(N ) model OPE
In this section we will consider the OPE of the ordinary bosonic O(N ) model, and in particular show how at lowest order in the large N expansion it reduces to the correct correlation function found easily from large N methods. The OPE of the 2-point correlation function of n in this model was first considered in [5] using an entirely different approach, and it was also considered in [6] using a background field approach similar in principle to that considered here. However the approach described here will differ in that we will use a background field method first considered by Polyakov [17] that maintains manifest O(N ) invariance in the background field at the cost of introducing O(N − 1) gauge redundancy. Polyakov's manifestly symmetric background field method can be used to easily find renormalization in the Wilsonian style for any value of N (there are some problems with the more straightforward approach of [6]) and it will allow us to identify operators in the OPE more easily, so it is worth spending some time reconsidering the bosonic case from this perspective before we move on to the more challenging supersymmetric case.

Polyakov-style background field method
The idea of a Wilsonian background field method is to separate the n field in the O(N ) model into a high frequency "quantum" field ϕ and a low frequency "classical" field n 0 . If we integrate out the high frequency field ϕ we find how the action for n 0 is renormalized. Perhaps less often understood is that if we take the expectation value of some correlation function of n and integrate out ϕ then we can find the OPE of that correlation function. We will split up n into the fields n 0 and ϕ in the manner of Polyakov [17], rather than the more straightforward linear manner n = n 0 + ϕ. Recall that the constrained action (with no λ field or m 2 term) is Now we will expand n i in terms of a spacetime dependent basis e i α (x), and refer to the components of n in this basis as ϕ, , where both the flat basis indexed by i and the space dependent basis indexed by α range from 1 to N .
The Lagrangian becomes Now we will partially fix our choice of e α by taking the N th component e N to be the low frequency coarse grained n field, which we call e N ≡ n 0 . We will use Latin indices for the remaining N −1 components e a . For convenience and consistency of notation with treatments of the O(N ) model elsewhere, we refer to the N th component of ϕ as not to be confused with the Lagrange multiplier field σ appearing in a different context in this paper. Now the Lagrangian may be written where D is the covariant derivative with respect to the O(N − 1) gauge symmetry of rotation of e a . If we further rescale ϕ to absorb the factor of 1/g 2 , and use the constraint to fix σ ϕ in terms of ϕ a , we have The first term is just the action for the background field. The second set of terms is quadratic in ϕ. The final set of terms involves corrections higher order in g 2 . Note that the equations of motion e · ∂ 2 n = 0 were used to throw away a term linear in ϕ. Of course since n 0 is integrated over in the path integral it will not always satisfy the equations of motion, but for our purposes this would only matter if there would be operators in the OPE involving factors of e · ∂ 2 n, and that is certainly not possible for low-dimension operators and at low orders of the large-N expansion.

Lowest order in large N
This expansion is dramatically simplified by considering the large-N limit, in which g 2 ∼ N −1 , so each factor of g 2 can only appear if it is compensated by a summation of indices producing a factor of N . In particular we see that corrections involving the higher order interaction terms will vanish at lowest order, since to compensate the factors of g we would need to contract the ϕ b ∂ϕ b factors together, but the single derivative inside a loop would cause this to vanish. Now let us consider the correlation function for the full field n. So as not to clutter our notation too much, the angled brackets will indicate integrating over the UV fields ϕ first, and then afterwards integration over the IR background field n 0 is implied.
The cross terms are those terms involving expectations of ϕ σ ϕ , and these vanish at lowest order for the reason just mentioned, since they require the use of the odd interaction terms. Now let us consider the first term The higher order O(N − 1) gauge dependent terms in the expansion of e a (x)e b (0) should in principle cancel the gauge dependence of ϕ a (x)ϕ b (0). But at lowest order in the large N expansion we will need a factor of N to compensate the overall factor of g 2 , and that may only be provided by the Kronecker delta term contracting with a corresponding Kronecker delta in the ϕ a ϕ b correlation function. So to produce this factor of N , we need only consider g 2 ϕ a (x)ϕ b (0) δ ab evaluated according to the reduced quadratic action where all ϕ are contracted with each other rather than factors of e, For consistency with the later treatment of the SUSY O(N ) model (where the analogous substitution will be more useful), the equations of motion (2.7) may be used to rewrite (∂n 0 ) 2 in terms of a background Lagrange multiplier field λ 0 . The term −m 2 in (2.7) may be considered the zeroth order VEV of λ 0 , which is more in keeping with the OPE philosophy.
So the reduced quadratic action is now In momentum space we have Here we are already taking into account the fact that λ 0 will be integrated over in the path integral, so the net external momenta introduced into the ϕ correlation function by operators constructed out of the background field is zero. A Feynman diagram notation that indicates integration over both ϕ and λ 0 may be helpful in clarifying this point, and is shown in . So at lowest order in large N we can completely ignore the effect of external momentum coming from the background field, (3.6) Figure 1: Diagrams involved in calculating the ϕ correlation function. A ϕ field is indicated in black, and the background λ 0 field as a red curly line. A single λ 0 insertion is indicated on the left. After integration over the background fields this will contract in a tadpole diagram (the trivial O N 0 VEV is shown as a stub) and there will be zero net external momentum introduced into the ϕ correlator. The middle and right diagrams show two λ 0 insertions, as in the third term on the RHS of (3.5). In the middle diagram each factor of λ 0 contracts with itself much like the left diagram. In the right diagram the different factors contract with each other, introducing momentum k into the ϕ correlator.
Now, since λ 0 (0) = −m 2 this is almost already our full OPE at lowest order in large N . But we also must consider the effect of the second term of Eq. (3.4), σ ϕ (x)σ ϕ (0) n 0 (x)·n 0 (0). So far we have not used the cutoff µ in any way, but the step function in (3.6) indicates that only this second term will matter for IR momentum p 2 < µ 2 .
Recall that the definition of σ ϕ after we impose the constraint is (3.7) At lowest order in large N each factor of ϕ a ϕ a must contract with itself, so we simply have This factor involving logs is just accounting for the field renormalization of n 0 compared to n, so that this entire term is just equal to the correlation function of the original field n at momenta p 2 < µ 2 as it should. It is easy to see that this is the field renormalization from Polyakov's original calculation that this background field method is based on [17]. An equivalent way to think of this is that after we take the expectation value of the n 0 correlator in momentum space, we get where g 2 µ is renormalized down to the new cutoff µ according to the renormalization group at lowest order in large N , So these log factors coming from σ ϕ just serve to renormalize g back up to the original cutoff M . Altogether, considering both Eq. (3.6) and Eq. (3.8), the OPE indeed produces the correct n correlation function N g 2 p 2 + m 2 , both above and below the scale µ. Note that we could formally take µ → 0, in which case the inverse powers of p 2 in the UV part of the OPE (3.6) would need to be regularized according to some scheme such as Hadamard's finite part (see e.g [28]). And in the IR part of the OPE (3.8), the factor n 0 (x) · n 0 (0) could be expanded in a Taylor series in x producing a series of derivatives of delta functions in momentum space. This is the form that the lowest order OPE takes in [5].

SUSY O(N ) model OPE
Now this entire framework will be extended to the SUSY case. While in the large N expansion it is useful to consider Lagrange multiplier fields, here we will leave the action as with constraints The variations of the fields need to be consistent with the constraints: And then given an arbitrary vector e such that e · n = 0, we have the following equations of motion,

Symmetries and Polyakov-style background fields
The constraints of the action (4.2) are obviously invariant under O(N ) transformations of the fields, but there are also invariant under two other not so obvious transformations. First we may leave n, ψ fixed and transform F by an arbitrary bosonic antisymmetric tensor We may also leave n fixed but transform both ψ and F by an antisymmetric fermionic matrix Using the previously defined superfield Φ in (2.9), all of these transformations may be written in terms of a single superfield E, So that the constraints (4.2) in superfield form, are preserved under the transformation, as long as the E superfield obeys the orthogonality constraint, which in component form is Note that if e is a general orthogonal matrix (not the identity), then d and f need not be antisymmetric. Now we have the machinery to do the equivalent of the trick Polyakov used in the bosonic O(N ) case. Expand Φ a in terms of a spacetime dependent orthogonal superfield basis Here to better distinguish the fields with a gauge dependent index α we will change our notation for the component fields (just as we called the components ϕ rather than n in the bosonic case), As in the bosonic case, this is just an arbitrary rewriting of our fields in a spacetime dependent basis, but it will be connected to the physics of renormalization by choosing the N th component of E to be the IR background field Φ 0 , Further, we will rewrite the N th component of Φ α in terms of σ fields which are solved in terms of constraints, So after rescaling the UV fields by g, the decomposition of Φ into IR and UV fields takes the component form, where the variational fields δn, δψ, δF are defined as and they are so called because they obey the variational constraints Eq. (4.3) together with the background field Φ 0 .
The constrained σ fields are solved in terms of these variations, (4.10)

Background-field Lagrangian
Now we will simply expand the Lagrangian (4.1) directly in terms of (4.8), using all the identities we've gathered so far, including the equations of motion for the background fields n 0 , ψ 0 , F 0 . As for the bosonic case, only the lowest order in g 2 will end up mattering at lowest order of large N , so we will only expand to this order. The quadratic terms coming from the bosonic term in the action is the same as for the bosonic O(N ) model in Eq. (3.2).
For the fermionic term, The final term is coming from the cross term involving both σ q and σ ϕ at lowest order. Now expanding δψ using the Eq. (4.9), To remove dependence on the fermionic gauge field d, a gauge transformation q − δn · d → q was used. This should also lead to a transformation of Q, but since we will be integrating out Q shortly, we will not show this explicitly. The Q fields appear in the Lagrangian through the term In Eq. (4.9), δF is decomposed into a term involving Q which is perpendicular to n 0 and a term parallel to n 0 which is simplified by the same gauge transformation above. Once Q is integrated out, only this term parallel to n 0 survives,

One-loop renormalization
The background field action given by substituting Eq. (4.11), (4.12), (4.13) in Eq. (4.1) is rather complicated, even at quadratic order, and in order to test it we shall make a brief detour and show that it gives the correct one-loop beta function, with no large N simplification involved. First note that the terms mixing ϕ and q in (4.12) will not be involved in the calculation since they can only appear in a loop leading to a dimension 3 operator at best, and this nonrenormalizable operator will be suppressed by an inverse power of the UV cutoff M . Likewise, the gauge fields in the covariant derivatives can only appear in a gauge-invariant combination at dimension 4, so for the same reason we may treat all covariant derivatives as ordinary derivatives.
The renormalizable contributions are as follows (use is made of the equation of motion e · F 0 = 0): • The terms proportional to ϕ 2 or mixed index ϕ a ϕ b can be integrated over straightforwardly, producing • A loop from two vertices of theqq term in (4.13) produces • A loop from two vertices of the ϕ a ∂ϕ b term in (4.12) produces • And finally, the loop from two vertices of the mixed index q a q b term in (4.13) vanishes, but the loop with one mixed q a q b and one contractedqq produces So in total, our renormalized action is leading to the correct beta function, (4.14)

Lowest order in large N
Once we consider correlation functions of n and ψ in the large N limit, much of the same discussion in bosonic case applies. Only favorably contracted indices will matter at lowest order, so we can further simplify the quadratic action described by (4.11), (4.12), and (4.13), This Lagrangian can be written in a more suggestive way if we introduce Lagrange multiplier fields σ 0 , u 0 , λ 0 in the background field action as in (2.13). Then, using (2.11),(2.17), and (2.20) we reduce the action to After a long journey, this looks rather similar to the action in terms of Lagrange multipliers we started with (2.13). Of course ϕ, q are unconstrained fields with N − 1 components, and there are all sorts of terms appearing in the full quadratic action found in Section 4.2 which are not shown here, since they will not contribute at lowest order in large N to the coefficient functions in the OPE. It may seem somewhat against the spirit of the OPE to explicitly indicate terms involving m, but this m may just be considered as the lowest order term of the σ 0 operator and in diagrams it will be indicated by a σ 0 line stub in the same manner as was done for λ 0 in the bosonic case in Fig 1. As a hint of what is to come, this form of the background field action will indicate the source of the terms in the OPEs for the n, ψ propagators which come from the lowest order in coefficient functions and the subleading order in operator VEVs (namely the quantities D (n) defined in Section 6). These terms will come from expanding the n and ψ propagators in the diagrams of large N perturbation theory in powers of m and loop momenta, while leaving the Lagrange multiplier field propagators unexpanded.

Operators in the OPE
Now we will return to the problem of finding the OPE associated to the correlation functions of n and ψ. This involves expanding the correlation functions in terms of ϕ and q using (4.8) and (4.9). As discussed in Section 3.2, if only the lowest order coefficient functions are needed then cross terms in n 0 and e 0 may be neglected, and e a 0 (x) · e b 0 (0) may be replaced with a Kronecker delta δ ab . Furthermore the σ 2 q ∼ O (g 2 ) term in the expansion of ψ may also be clearly neglected at this order. The n 0 ψ 0 · e a 0 ϕ a term in the expansion of δψ may also be neglected due to the extra factor of g 2 that would come from the ψ 0 propagator. So if we are only concerned with the zeroth order coefficient functions, the correlation functions may be expanded as As previously discussed for the non-SUSY (bosonic) case, the factor σ ϕ may be considered to be the field renormalization and those terms describe the IR behavior of the correlation functions for momentum less than µ. So we will instead focus on the terms involving the correlation functions of ϕ and q.

OPE at subleading order
The OPE for the correlation function of n at subleading order is schematically where the subleading coefficient functions C (1) j in the first term will be calculated in Section 6. This section will focus on the operator VEVs in the second term.
Considering the factored action (4.16), the zeroth order coefficient functions and operator VEVs up to engineering dimension 4 are n(−p) · n(p) (1) The diagrams associated with these terms are presented in Fig 3. Note that a possible operatorū 0 u 0 with dimension 3 has a vanishing coefficient function, and an m insertion in the fermion line is necessary. As discussed for the bosonic case, expanding the loop momenta in the ϕ, q lines in Fig 3 would lead to derivatives acting on the background field operators, but this is first seen at dimension 6.
In the same manner, the first few orders of the fermion correlation function OPE may be calculated, Here the fermion indices in u 0ū0 are contracted with the surrounding factors of ¡ p.

Operator ambiguities
The subleading coefficient functions C (1) j which were not shown explicitly in equations (5.3) and (5.4) above will have some IR renormalons which will lead to imaginary ambiguities which depend on the choice of integration contour in the Borel plane. However the full correlation function which may be calculated exactly in the large N limit is of course unambiguous. So the IR renormalon ambiguities must somehow cancel in the full OPE, and in fact they will cancel with ambiguities arising in the precise definition of the operator VEVs. This was first pointed out in the context of the bosonic O(N ) model by David [18,19] (see also the discussion in [10]).
In the presentation of the OPE so far, we have been using a hard cutoff µ > m separating UV fields from IR background fields as in [6]. However if this IR cutoff on the UV fields is taken seriously in calculating the coefficient functions no IR renormalons are seen, so the problem of ambiguities never arises. Instead it is the arbitrary cutoff µ that must cancel in the full OPE [6] (see also [9]). This lack of dependence on µ was seen explicitly at lowest order in the bosonic O(N ) model above in Section 3.2, although no IR renormalons would be seen in this simple example at lowest order in large N .
However if we formally take the limit µ, m → 0, or use dimensional regularization, the IR renormalons show up in the coefficient functions. And so the corresponding ambiguities in the operator VEVs must be taken seriously as well. The operator ambiguities will be calculated in a manner closer in spirit to the way they were originally calculated by David [18] in Section 6 below using the quantities D (n) defined there. Here we will instead discuss an auxiliary method to calculate these ambiguities using a cutoff regularization which was also first proposed by David in [18] in terms of lattice regularization. This will be illustrated by considering the VEV of σ 2 0 , which appears in the OPE of both correlation functions (5.3) and (5.4). Finding the ambiguity of the σ 2 0 operator will be very similar to manner that the λ 2 0 operator of the bosonic model is considered in [10] and [9]. Using the σ 0 propagator (2.15), the VEV of σ 2 0 to this order is , and the background field σ 0 has a UV cutoff µ. The main idea of the calculation is as follows.
If we calculate the leading power law UV divergence in µ we will in fact find a divergent series in g which has a Borel ambiguity. Since the original expression is unambiguous, the ambiguity of the divergent series must in fact be canceled by a corresponding term in the part of the VEV which is not power law divergent. Now if we instead consider dimensional regularization, then this power law divergent piece with a Borel ambiguity never appears, however the ambiguity in the part of the VEV which is not power law divergent remains. So the first step is to find the power law divergence in µ, which will come from the domain of integration where k 2 m 2 . Using the asymptotic expansion of J (2.5), Here any positive powers of m have been disregarded, and any factors of m inside the argument of a logarithm have been expressed in terms of g using the saddle point relation (2.2). For simplicity we may define a rescaled coupling constantĝ, and change the variable of integration,ĝ This series is clearly divergent, but it can be easily treated by taking the Borel transform, defined for a power series inĝ with arbitrary coefficients r j as and the corresponding inverse Borel transform of a function f (t) is defined as So then, which has a pole on the positive t axis so the inverse Borel transformation can not be applied without further consideration. The t integral may be regularized by deforming the contour of integration either above or below the pole. Using the definition ofĝ and the saddle point condition for g, the ambiguity in the inverse Borel transform is calculated as This is the ambiguity in the divergent series, but again the basic idea is that there is a corresponding ambiguity of opposite sign in the terms of O (m 2 ), and that is what remains when dimensional regularization is considered. To summarize, using the notation of curly braces to indicate the ambiguity at order 1/N , the ambiguity in the VEV of σ 2 0 is, The VEV of λ 2 0 also appears in the OPE (5.3), and can be calculated similarly. In fact since the propagator is the same, the calculation is exactly the same as the bosonic case considered in the references cited above [18,10], The λ 0 and σ 0 tadpoles also appear in (5.3) and (5.4) respectively. But after summing the contributions from the distinct tadpoles as in [20], This involves no power series in g 2 so there is no ambiguity, Finally there is the operator u 0ū0 (and its negative traceū 0 u 0 ), which has a VEV that can be easily related to σ 2 0 due to the similarity in propagators (2.15).
So now the total ambiguity in the OPE may be calculated for the first few operator dimensions.
{ n(−p) · n(p) } = {C where plus and minus are used to distinguish the coefficient functions of the bosonic and fermionic OPEs, respectively. The correlation functions on the left hand side of these equations are something unambiguous, and the ambiguities arising from IR renormalons in the coefficient functions on the right hand side should cancel with the ambiguities calculated here for the operators. Indeed in the following section this will be shown explicitly, and cancellation of ambiguities will be shown further to all orders in m.

Coefficient function asymptotic expansion and ambiguities
As shown in the last section for the ambiguities in operators, similar ambiguities are also revealed in the coefficients of an OPE. In this section, by looking at the large N limit of the supersymmetric O(N ) model, we demonstrate how the ambiguities of the coefficient functions emerge and their cancellation with those from condensates to all orders as well. In fact, a similar cancellation was illustrated in [7,23] for two cousin models, the non-supersymmetric (bosonic) O(N ) and Gross-Neveu model, in the large-N limit. First of all, let us look into the exact form of the two-point correlator of the n fields (on the Borel plane) and identify the contributions of the coefficient functions and the VEVs of operators. We will also discuss the same situation in the fermionic sector later on.

Effective mass and perturbation series
Within the scope of this paper, we focus on the large N limit of the theory and perturbative results here are obtained through the expansion around the saddle point. In Section 2, Eq. (2.2) gives the IR scale and the physical mass m 2 at the leading order. At this point, m 2 is written in terms of the UV scale M with the bare coupling constant g 2 , which introduces a natural parameter 1/ log(p 2 /m 2 ) for a perturbation series.
Yet, to study the ambiguity structure, it is not enough to stay only at the leading order as already mentioned in previous sections. We have to go one step beyond to look into the 1/N correction to the perturbation series. In this case, the two-point correlation functions, especially those of n a and ψ a , are further modified and different poles corresponding to the physical mass emerge. To be more concrete, in the subleading order, the propagator of n a s turns out to be n(p) · n(−p) = N g 2 p 2 + m 2 + n(p) · n(−p) (1) (6.1) where the latter term is shown in (6.4). The next-to-leading order term can be brought to the denominator of the propagator as a geometrical series and gives rise to the pole where m 2 phys is the corrected physical mass. [Comparison of (6.2) with the physical mass following from the β function is straightforward, the latter expression coincides with (6.2) up to a non-logarithmic factor.] Note that the arc terms in Fig. 4 can be safely ignored since they are proportional to p 2 + m 2 , which vanishes to leading order at p 2 = −m 2 phys , and so only the tadpole graphs contribute. The logarathmic term proportional to −2/N is precisely what is needed to correct the factor of N in (2.2) to the N − 2 that appears in the beta function (4.14).
From another perspective, the consistency condition with respect to the supersymmetric ground state λ 0 = 0 [20] suggests the same modification of the physical mass m 2 phys . This indeed coincides with the previous observation since the 1/N correction to λ 0 is nothing but considering the tadpole corrections. However, for our current purpose of studying the renormalon poles in SUSY O(N ) model, expanding in terms of m 2 phys makes these poles obscure so the perturbation series will still be expressed in terms of m 2 .

Bosonic propagator
To begin, the first order correction in the large N expansion [20] can be found by considering the Feymann diagrams in Fig. 4. This reads where J(k 2 , m 2 ) was defined in Eq. (2.4) and the n fields are convoluted for simplicity. The latter term is nothing but the tadpole contribution and does not contribute to the renormalon structure in the present case. For the time being let us focus on the first term in the square bracket and denote it as I 1 (p). To study the analytic structure of n(p) · n(−p) (1) , following [7], we can express the logarithm in the denominator in I 1 (p) as which can be put in a Mellin-Barnes representation, where the kernel function K(s, t) is defined as Under the above manipulations, n(p) · n(−p) (1) can be expressed in terms of a power series in m 2n /p 2n+2 . To see this is the case, we first evaluate the k-integral and compute the sintegral by closing the loop on the left and picking up the residues enclosed in this contour (see appendix B in detail). The final form of n(p) · n(−p) (1) has the structure (with tadpole terms ignored in the following discussion) where g 2 (p) is the effective coupling constant defined via the leading order saddle point equation .

The complete expression of A (n) [t], B (n) [t], D (n)
[t] functions are given in appendix B. Note that (6.6) can be interpreted as the Borel transform of the perturbative expansion of the correlation function n(p) · n(−p) (1) , especially A and B functions, due to the exponential prefactor of e −4πt/N g 2 . Then, what is the physical origin of these contributions and how can we connect the current presentation of the two-point correlator to the notion of the usual OPE? In the OPE expansion of n(p) · n(−p) , schematically it takes the form where C

Relation to coefficient functions and condensates
In this section we will be more concrete about the connection between the asymptotic expansion (6.6) given above and the usual OPE in terms of coefficient functions and condensates of operators.
First, let us discuss how the D (n) [t] terms in (6.6) are related to the condensates of operators in the OPE. Because D (0) [t] is mostly related to the renormalization, we can neglect it for the time being and focus on the meaning of D (1) [t]. Then, to get some insight on D (1) [t], let us write down its explicit form given in (B.3), which implies an IR renormalon pole only at t = 1. On the other hand, we can expand n(p) · n(−p) (1) in the limit p k ∼ m before the arc integral in (6.4) is carried out and it gives at the zeroth order 1 n(p) · n(−p) (1) in which the middle expression is exactly of the same form as σ 2 up to a prefactor in the large N theory. Indeed this is just the contribution from the σ 2 loop diagram in Figure  (4). Considering the form of (6.6), we see that this is identical to D (1) [t] at the lowest order in m 2 /p 2 . This observation signals that D (1) [t] can be associated to the condensate σ 2 . Indeed, taking the prefactor into account, the ambiguity given by choosing a contour above or below the pole at t = 1 is identical to that found in (5.5) using the more indirect method of considering the leading UV divergence. A similar consistency check may be done for the other operator ambiguities calculated in Section 5. Note that by expanding the side propagators (p 2 + m 2 ) −2 in powers of m 2 , we see that D (1) [t] must also contribute to the higher order operators σ 2n in the OPE and this is closely related to the factorization of the VEVs of the operators in the large N limit. This is also consistent with the presentation of the background field method diagramatically in Figure  3. The expansion of the side propagators is equivalent to arbitrary insertions of m 2 stubs in the ϕ propagator, while D (1) itself corresponds to a closed σ 2 0 loop, which is consistent with its origin in terms of the σ 2 loop diagram in Figure 4.
Of course, the reason we have been considering the asymptotic expansion (6.6) in the first place is to find the ambiguities of the coefficient functions in the OPE {C (1) } which are needed to cancel with the ambiguities of the operators given in (5.10). These coefficient function ambiguities come instead from the quantities B (n) in (6.6), where only relevant terms related to the coefficient functions are listed and the overall factor of (p 2 + m 2 ) was also expanded in terms of m 2 . The quantities proportional to the zeroth order VEVs 1, m 2 just represent the coefficient functions C +,0 , C +,2 . Note that although we are considering the cancellation of ambiguities up to dimension four in powers of m we do not need to include the coefficient function C +,4 since the renormalon ambiguity always introduces at least one extra power of m 2 . Using the explicit expression for B (n) in (B.3) the ambiguities are found to be,

±iπ) and {C
(1) In the case of C +,0 (i.e. the coefficient function of the identity operator), the ambiguity comes from the pole located at t = 2 in B (0) while in the case of C +,2 (i.e. the coefficient function of σ 2 0 ) the ambiguity includes the pole at t = 1 from both B (0) and B (1) . It can be easily seen from (6.10) and (5.10) that the two contributions do indeed cancel with each other. The meaning of this cancellation is slightly different from the cancellation we show in Section 6.4. In the present case, the vanishing of ambiguities is involves the usual notion of the OPE and can be shown order by order in the series in m 2n /p 2n+2 while the all-order proof below mixes contributions from different operator dimensions in the OPE due to the overall prefactor (p 2 + m 2 ) −1 in the definition of the quantities A, B, D.

Cancellation of IR renormalons
Now let us demonstrate how the manifest cancellation of the IR renormalons shows up at all orders. To start with, consider the functions listed in (B.2) and (B.3) and notice some of these functions have singularities at t = ±n 0 for n 0 ∈ N. We know that the ambiguities on the positive real axis correspond to the IR renormalons while the UV renormalons are related to the poles on the negative real axis. 2 The IR renormalons are what lead to ambiguities so we will focus on these. First note that for t > 0, no singularity appears in A (n) [t] for any n, so it will not lead to any ambiguities at all.
As for the poles in B (n) [t] and D (n) [t], to start with, notice that B (n) [t] has poles for all positive integers whereas D (n) [t] only has poles for positive integers less than or equal to n. However these poles can be paired in a one-to-one way where the pole of B (n) at t = n 0 is associated with the pole of D (n+n 0 ) at t = n 0 . It is easy to see that the residues of these poles 3 are of the same order in m and so have the potential to cancel. To be clear, the residue of B (n) [t → n 0 ] will lead to something proportional to which is clearly the same order as D (n+n 0 ) [t → n 0 ] which lacks the exponential factor.
2 These functions also have poles at t = 0, but it corresponds to no renormalon effect but the logarithmic UV divergence and can be canceled by renormalization [7]. Indeed, as transformed back to the series of g 2 (p), 1/t contributes as a constant term because t = 0 in e −4πt/N g 2 brings us no g 2 factor. 3 The residue and the ambiguity ±iπ are one and the same thing since they all emerge from the notion of the contour integration around a (simple) pole.
What remains here is to see that B (n) [t] and D (n+n 0 ) [t] do in fact lead to opposite contributions. To this end, we find the residue of B (n) [t] at t = n 0 is where we used the fact that 1/Γ(−n) vanishes for n a non-negative integer, so in the last line the sum of k is only from 0 to n. For simplicity let us denote what we get in (6.11) as R n 0 and the total contribution is We see that the IR renormalon poles indeed cancel order by order in the exact two-point correlator of n a . As an aside, the poles on the negative real axis do not cancel between these terms. For example, we can look at t = −n − , but we have to consider two different cases: n > n − and n ≤ n − . Then, for the former case, the residue in which again have the same cancellation as the case of poles on the positive real axis. 4 On the other hand, B [n] [t] again has a non-zero singularity at n ≤ n − but there is no term to cancel this pole. To see this is the case, we go back to the dimensional argument that B (n) [t] at n ≤ n − is accompanied by a factor in which the power is negative and D (n) [t] does not pair with any negative power of −m 2 /p 2 in (6.6).
ψψ ψψ ψψ Figure 5: The leading order correction to the Green function of ψ field in large N limit. All lines represent the same fields respectively as those in Fig. 2.

Fermionic propagator
Now, let us turn our attention to the fermionic correlation function. To proceed, the first order correction to the fermion two-point function comes from Fig. 5 and takes the form where I 1 is the exact same integral that was defined for the bosonic correlation function in (6.5). So there should be an identical expansion in terms of quantities A (n) , B (n) , D (n) , and since the tadpole term proportional to σ 0 has no ambiguity as indicated in Section 5, the renormalon ambiguities cancel to all orders in this OPE as well. However, although they are unchanged, the quantities A (n) , B (n) , D (n) will have a different interpretation in terms of operators in the fermionic case. The terms leading to ambiguities in the fermionic OPE up to O (m 3 ) are, Comparing this with (5.4) we see that the expression D (1) is associated to the operator condensate for both for σ 2 0 (as in the bosonic case) but also to the particular combination of dimension-3 operators mσ 2 0 and u 0ū0 at next order. This leads to no inconsistency since the ambiguity in (5.11) was indeed found to be proportional to ¡ p + m. Similarly the coefficient functions of both the identity operator and the operator σ 0 are both expressed in terms of B (0) [t], and the ambiguities to lowest order in m are given by the pole at t = 1. So using the same results for the bosonic case at dimension-2 we can find the renormalon ambiguities of the fermionic coefficient functions and indeed these cancel with the operator ambiguities found in (5.11).
7 Further discussions and conclusions

Factorial divergence of perturbation theory in various models
There is a profound difference between the perturbative expansion, say, for the energy eigenvalues of an anharmonic oscillator and in problems arising in asymptotically free field theories at strong coupling. In the former case a coupling constant is well-defined. If anharmonicity is small, g 1, perturbative series in g which are usually plagued by factorial divergences in high orders can be made well-defined based on the quasiclassical data obtained in the complex plane. In this way one arrives at trans-series, including both regularized perturbation theory and exponential terms e −c/g in a systematic manner (for an introductory review see [29]).
In Yang-Mills theory there is no dimensionless coupling constant to form a perturbative expansion in the strict sense of this word. If we ignore quarks for the time being then the only parameter of the theory is the dynamical scale Λ, due to dimensional transmutation. The only expansion parameter appearing in the 't Hooft limit [30] is 1/N where N is the number of colors. Quantitative methods for construction of perturbative series in 1/N have not yet been developed although a number of qualitative observations exist. If such a series could be obtained, say, for the mass of the lightest glueball, we would arrive at where c j are purely numerical dimensionless coefficients depending only on the quantum numbers of the glueball under consideration. Needless to say, the question of convergence in (7.1) at large j will arise. Exponential terms of the type ∼ exp(−cN ) will appear but we will not discuss the 1/N series in full in this work limiting ourselves to the leading and the first subleading terms. Passing from Yang-Mills to QCD with massless quarks aggravates the situation. As is well-known, spontaneous breaking of the continuous chiral symmetry (χSB) is not seen in perturbation theory in the gauge coupling, whatever this coupling might mean. There is a range of questions in which QCD perturbation theory is widely used, however.
In QCD and similar theories it is quite common to add external sources to use them as tools and consider various correlation functions at large Euclidean momenta p. Thus we acquire a large external parameter p/Λ and can develop a perturbation theory in α(p) ∼ (log p/Λ) −1 . Here α(p) ≡ g 2 4π . Even in the problems with a large parameter p/Λ the α(p) expansion cannot be made closed -i.e. cannot be continued to any desirable accuracy, as is the case in the quantummechanical trans-series. The problem is as follows: the α(p) expansion is intrinsically illdefined [31] because even if α(p) 1 any Feynman diagram saturated at k ∼ p still contains contributions from virtual momenta k ∼ Λ (see Fig. 6). In this domain the coupling α is not defined, simply because the Lagrangian formulated in terms of quarks and gluons ceases to exist.
The resurgence and trans-series program (such as in quantum mechanics) can not be fully successful in QCD-like theories at strong coupling. Underlying dynamics in confining theories at large distances in no way reduces to expansion in α even being supplemented by additional quasiclassic analyses. Confinement of the Nambu-Mandelstam-'t Hooft type k p Figure 6: The bubble-chain diagrams for the Adler function D representing renormalons. Wavy lines stand for an external "photon". Solid lines denote quark propagators, while dashed lines are for gluons. The quark bubbles are also to be added to the gluon bubbles.
was demonstrated [32,33] to emerge from the dual Meißner effect -a very special nonperturbative feature of the Yang-Mills vacuum -and so is χSB. Both are crucial at distances Λ −1 and leave no trace in perturbation theory.
In the absence of the solution of strong coupling Yang-Mills theories in 4D, the best one can do is the OPE in the form described in detail in the reviews [9] (see also references therein) which requires a new delimiting parameter µ, Above it was referred to as the normalization point. It is expected that the coefficient functions saturated by virtual momenta larger than µ can be calculated analytically while the operators saturated in the infrared (IR) domain below µ are either parametrized or obtained numerically. This procedure is in one-to-one correspondence with the Wilsonian renormalization group approach [34,27]. The parameter µ drops out from measurable correlation functions, for instance, the Adler function D(p 2 ),

4)
J µ is a conserved vector current of a massless quark and p is the external momentum in Euclidean space.

Generalities of OPE and renormalon conspiracy
Wilson's OPE for D(p 2 ) has the form × (possible logs) (7.5) where d j is the dimension of the operator O j normalized at µ 2 . The left-hand side of (7.5) is µ independent, and so should be the right-hand side. This is quaranteed by µ conspiracy between the coefficient functions and the operators in OPE. In the O(N ) model this fact was explicitly checked in [6]. Note that if µ is kept finite in the interval (7.2) there are no renormalons in the coefficient functions and the OPE operators are unambiguous.

Renormalons
The issue of renormalons as a source of factorially divergent perturbative series in Yang-Mills theories was raised by 't Hooft [31] (see also [35], detailed reviews can be found in [36,10] and the second paper in [9]). There is a formal assumption in 't Hooft's original consideration which cannot be justified, however in confining theories of the type of QCD, see below. 5 It was noted that a special class of the so-called bubble diagrams (Fig. 6) formally lead to a factorial divergence of the coefficients of (g 2 ) j at large j. It is worth explaining this in more detail.
After one integrates over the loop momentum of the "large" fermion loop and the angles of the gluon momentum k in Fig. 6 one arrives at where the function F (k 2 ) calculated in [37] has the IR limit Next, assume one combines the above expression with the standard formula for the running coupling constant 4π log(p 2 /k 2 ) (7.8) (β 1 is the first coefficient of the β function, β 1 = 11 3 N − 2 3 N f , α(p 2 ) is considered to be fixed) and perform integration over k 2 (7.9) expanding the denominator in (7.8) in powers of α(p 2 ). Note that formally integration runs from k 2 = 0 implying that the denominator in (7.8) hits zero at k 2 = k 2 * , and at k 2 < k 2 * the denominator is negative. Equation (7.9) can be identically rewritten as The series (7.10) is asymptotic and diverges at n → ∞ but this is obviously an artifact of using Eq. (7.8) in the domain of k 2 where it miserably fails. The confining regime at large distances in Yang-Mills theory at strong coupling implies that Eq. (7.8) is totally inappropriate at k 2 < ∼ Λ 2 since α(k 2 ) is not defined in this domain. Equation (7.8) is valid provided 4π log p 2 /k 2 < 1 (7.11) which in turn requires k 2 > µ 2 , µ 2 = cΛ 2 , c 1 .
The IR cut off of this type must be imposed in bona fide QCD OPE, see e.g. the review papers [9]. Then the factorial growth of the coefficients is cut off at a critical value 13) and the series in α(p 2 ) in the OPE coefficients is well-defined, i.e. convergent (see Fig. 9 in the second paper in Ref. [9].) The contribution of the domain k 2 < µ 2 goes into the matrix element of the gluon operator which is not calculable. The situation drastically changes in the exactly solvable models in which the two-point function under consideration is explicitly known. It satisfies all physical requirements. In principle, there is no need at all to expand it. If, however, we do so, this might prompt us to discover what is happening in 4D QCD.

Exactly solvable models
In this paper we have focused on 2D supersymmetric O(N ) models in the limit of large N . They can be solved in the leading and, to an extent, next-to leading order in 1/N . The exact formula for the two-point function of the n fields is unambiguous. One can readily establish that to the leading order in 1/N where m is the mass of the elementary excitation and M 0 is the UV regulator mass. 6 The exact formula for n a (−p) n a (p) in (6.4) depends only on the ratio p/m. Since this correlation function is known exactly and is µ independent, we can can consider OPE in the limit µ → 0. In this limit the separation of soft and hard contributions in Wilson's OPE becomes the separation between perturbation theory and non-perturbative effects. While the exact expression is unambiguous the above separation introduces ambiguities both in the coefficient functions and operators in OPE (e.g. at µ = 0 the 't Hooft IR renormalons indeed appear). The ambiguities must cancel each other. Dependence on p/m in the exact solution appears in a two-fold way. The exact formula contains logarithms of the type log p 2 /m 2 and powers of m 2 /p 2 . Of course, in mathematical sense m 2 /p 2 is just the exponent of log p 2 /m 2 . However, it is instructive to keep a double expansion n a (−p) n a (p) = j, where the first factor represents coefficients while the second matrix element in the limit µ → 0. In other words, from the exact answer we can define what can be called the "coupling constant" (to the leading 1/N order) At p 2 Λ 2 the definition in (7.16) coincides with the standard perturbative one in the leading 1/N order. We have shown that the sum (7.15) is well defined, which is guaranteed through the conspiracy in the renormalon cancellation between the coefficient functions and matrix elements.

A brief summary of our results
Now after this discussion on the relevance of our results, let us return to what we have concretely shown in this paper, which is summarized in more detail in Section 1.1. The problem of explicitly finding the operator product expansion in the large-N limit of the supersymmetric O(N ) sigma model in two-dimensions was considered. The related problem for the non-supersymmetric (bosonic) O(N ) sigma model has been considered long ago [5,6,7], but as far as we are aware this is the first time the operator product expansion of twopoint functions in the supersymmetric O(N ) model has been explicitly found, and we have demonstrated the cancellation of infrared renormalon poles between the coefficient functions and operators to all orders. What is more, the background field method developed here in Sections 3 to 5 gives arguably a cleaner and more direct expression for the operators in the OPE even for the bosonic case which has been thoroughly explored in the past.
The background field method herein is in principle entirely perturbative in g, and could be applied to find the OPE and coefficient functions even if we were ignorant of infrared theory and the values of the VEVs. An example of a calculation of a coefficient function in this manner is shown below in Appendix A. Of course this calculation does indeed take advantage of all the simplifications vector-like large N theories in two dimensions provide.
However we have also considered this problem from the large N perturbation theory side in Section 6 and Appendix B, where we first wrote down the well-known exact two-point functions, and then found the asymptotic expansion (6.6). This in some sense gives us the entire OPE data to all orders in m 2 at the 1/N level. This approach is similar in principle to what was done for the bosonic O(N ) model [7]. As was done there, we verified the explicit cancellation of the IR renormalons between the condensates and the coefficient functions. Thus, in this particular model -supersymmetric O(N ) in two dimensions -the challenge of supersymmetry is successfully solved.
We emphasize that we have a more complete picture of the OPE than that approach alone since we can associate parts of the expansion in terms of m 2 in terms of VEVs of specific operators in the OPE. This was done for the first few operator dimensions in Section 5. But as was made clear there, it can quite easily be extended to any order since it amounts to an expansion of the propagators of the n and ψ fields in the diagrams of Figs. 4 and 5 in powers of m and loop momentum, while leaving the propagators of the Lagrange multiplier fields unexpanded.
The most straightforward conclusion we draw from all this is that at least in 2D SUSY O(N ) there is really no tension between supersymmetry and renormalon poles in the OPE. When we extend the bosonic model to the supersymmetric case there are indeed some VEVs which must vanish, but there are new fields introduced so we have a wider set of operators available. And as shown here, these operators do indeed conspire to cancel with lower dimension renormalon poles. However, the challenge remains in four-dimensional super-Yang-Mills.
A Appendix: Coefficient functions from ordinary perturbation theory The asymptotic expansion of the exact correlation function to subleading order in the large N expansion given above in Section 6 is powerful in the sense that encodes information about the entire OPE. However in generic theories we would not have the exact correlation function at our disposal, and if we did, there would be little need for OPE methods. In the present section we will calculate a coefficient function directly from perturbation theory in g, in a manner that does not explicitly require knowledge of the IR behavior of the theory. Given the full background field Lagrangian in Section 4.2, the coefficient function of any operator could be calculated by including higher order corrections in g in addition to the background field insertions considered in Section 5. But here for simplicity we will consider only the identity coefficient C and we may also simplify the correlation function (3.5) The second term in the correlation function involving σ ϕ is necessary to cancel the IR divergences in the first term involving the N − 1 components ϕ a [24,25]. Strictly speaking ϕ may still be thought of as a UV field defined up to some arbitrary IR cutoff µ > m, beyond which it is more appropriate to consider the IR field n 0 via some non-perturbative method. The cancellation of IR divergences simply means µ can be taken arbitrarily small in a way which can be compared to asymptotic methods such as those in Section 6. Since it will be taken arbitrarily small anyway, for simplicity we will modify the IR cutoff to be a soft cutoff given by an ad-hoc mass term µ 2 ϕ 2 as is usual in perturbative treatments of the O(N ) model. Once the factors of σ ϕ in the interaction term of (A.1) are expanded as a power series in ϕ, there may be arbitrary powers of g 2 ϕ 2 on either side of the Laplacian ∂ 2 . It is convenient to represent this Laplacian in Feynman diagram notation as a dotted line that has "propagator" −p 2 /(2g 2 ), where p is the net momentum of the σ ϕ factors. A similar notation is used in [24,25], where the pairing of each factor of ϕ a ϕ a is also indicated explicitly in diagrams. Here this will not be necessary since large N considerations drastically restrict the relevant diagrams.
As in the discussion of the large N limit of the OPE in Section 3.2, since each power of g 2 ϕ 2 in the expansion of σ ϕ comes with a factor of g 2 ∼ N −1 the ϕ 2 factor must be contracted with itself as a tadpole to provide a compensating factor of N . If instead two separate factors of g 2 ϕ 2 are contracted with each other as a connected "bubble" then there is only one compensating factor of N so this is unfavorable. However, the g −2 factor in the dotted line propagator may compensate a single unfavorable bubble contraction, so we may form "bubble chains" in the large N limit as in Fig 7. Purely from considering factors of N in this manner the tadpole-like bubble chain on the left of Fig 7 would be expected to contribute at leading order in the large N limit. However there is no net momentum flowing through the dotted lines of the bubble chain, so in fact these diagrams vanish. These bubble chains may be routed in a loop in order to introduce a net momentum through the chain, but this gives up the factor of N associated to the contracted ϕ 2 at the end of the chain, so these diagrams first appear at subleading order in large N . The two distinct ways to form a bubble chain loop at this order are also shown in  Reserving one factor of :ϕ 2 : in each σ ϕ in the interaction term of (A.1) to either form a bubble or connect to the external legs and contracting the rest into tadpoles, we have the effective interaction term where the tadpole contributions naturally lead to the renormalized coupling constant The first diagram to consider is the arc diagram as in Note that this diagram originally had a power law divergence that was exactly canceled by the O (N −1 ) contribution of the Jacobian factor, which must be considered for regularizations with a hard cutoff M (see e.g. [24] for discussion of this Jacobian factor in the lattice regularization case). Now that the Jacobian factor is fully accounted for, let us proceed to consider the arc diagrams that involve at least one bubble in the chain, which lead to the correction  Even after fully accounting for the Jacobian, this still has a power law divergence coming from the region of integration where k is large, but this divergence is exactly canceled by the correction from the bubble chain diagrams on the right of . (A.5) Now finally we must also consider the corrections arising from the term g −2 σ ϕ (x)σ ϕ (0) in (A.2). There are corrections where each ϕ 2 term in the expansion of the σ ϕ is contracted with itself with or without possible vertex insertions, leading to g −2 1 − g 2 ϕ 2 .
When combined with the ϕ a (x)ϕ a (0) term in (A.2), this serves to ensure the constraint n 2 = 1 is satisfied.
Then there are corrections coming from singling out two factors of :ϕ 2 : in the expansion of g −2 σ ϕ (x)σ ϕ (0) and forming a bubble chain much as in Figure 7. The bubble chain can either connect a factor from σ ϕ (x) to one from σ ϕ (0) or it may connect two factors from the same σ ϕ in which case the diagrams associated to the factors σ ϕ (x) and σ ϕ (0) are disconnected as in the previous paragraph. Once again the disconnected contribution cancels with the connected contribution as x → 0, and serves to enforce the constraint. The connected contribution from the bubble chains connecting σ ϕ (x) to σ ϕ (0) in momentum space is .
(A. 6) A.2 Equivalence to the large-N theory So far we have found the leading and subleading contribution in large N to the identity coefficient function of two-point function in the bosonic O(N ) model from the ordinary g perturbation theory. In fact, this result coincides with the one found by expanding about the large N saddle point as in [21,7]. To see this is the case, let us first look at (A.4) and (A.5) and the lowest order expansion of their sum I BC in the limit p 2 m 2 . Note that the regularization I BC is done by applying a UV cutoff M (and we set the IR cutoff to be µ) , where I is the tadpole integral defined as This is in fact exactly what we get with the approach of Section 6 and Appendix B. The identity coefficient function that we are considering here can be found from the p m limit of (6.6) where only the lowest order term survives and takes a quite simple form n(p) · n(−p) (1) 0 ≈ − (A.11) Notice that the terms marked by underbrace in n(p) · n(−p) (1) 0 are identical to (A.10) and indeed match with our previous interpretation that the terms with the prefactor e −4πt/N g 2 originate from the UV regime as the coefficient functions should be. The additional terms which do not match should not cause us too much concern since we have only been focusing on terms in the identity coefficient which may lead to renormalons in this subsection. For instance in the bosonic case we began with above I BC does not include the contributions from (A.3) and (A.6) which do appear in the coefficient function but do not lead to renormalon poles.
To wrap up this section, we elaborate the relation between the location of the ambiguities and the g perturbation series with different normalization points. First, in the original derivation of the asymptotic expansion in Section B, we do not rely on the g perturbation theory or any sliding scale to find the Borel representation. If we adopt the procedure presented in the early part of this subsection, the different normalization point means the change of the integration cut, say, from 0 to µ and from µ to ∞. This is not as transparent as taking p as the sliding scale such that the asymptotic series can be written in a simple g series with some well-defined special functions. However, we can still perform the Borel transform first before carrying out the integration of the coefficients of g j (µ ) and study the poles along t-integral by analytic continuation. This eventually indicates that the location of the ambiguities does not change.

B Appendix: Details of asymptotic expansion
In this section, we present the derivation of the expansion of I 1 (p) given in (6.6). First in (6.5), we already resolved the complicated logarithm in the original integral and let us proceed to calculate the k-integral. Namely, where dimensional regularization is adopted with dimension 2 − 2 . Note that 2 F 1 (a, b, c, z) is the hypergeometric function and the linear transform formula is applied to analytically and the n-th order coefficient of the series can be found by taking the n-th power of m 2 /p 2 after the sum of residues for each j. We collect the full result in (B.3). Then, for the contribution in the other sector, we have to consider the residues at s = −1 − j, j ∈ N ∪ {0} with − Γ(2j + 3)Γ(−t)Γ(t + 1) Γ(2 + j + t)Γ(−t + j + 2) where we applied the identity Γ(t − j − 1) = (−) j+2 Γ(−t)Γ(t + 1) Γ(−t + j + 2) .
Again the n-th order coefficient can be extracted by considering the appropriate power of m 2 /p 2 after the residues are all summed. Lastly, we list the three functions in (6.6) found from the analytic method given above. That is, for n = 0, and for n ≥ 1, Γ(2n − 2k + 1)Γ(−t)Γ(t + 1) Γ(t + n − k + 1)Γ(−t + n − k + 1) Note that the regular part of D (0) [t] can take a different form if we adopt different subtraction schemes.