Renormalon Subtraction in OPE by Dual Space Approach: Nonlinear Sigma Model and QCD

It is becoming more important to subtract renormalons efficiently from perturbative calculations, in order to achieve high precision QCD calculations. We propose a new framework ``Dual Space Approach"for renormalon separation, which enables subtraction of multiple renormalons simultaneously. Using a dual transform which suppresses infrared renormalons, we derive a one-parameter integral representation of a general observable. We investigate systematically how renormalons emerge and get canceled in the entire operator product expansion (OPE) of an observable, by applying the expansion-by-regions (EBR) method to this one-parameter integral expression. In particular we investigate in detail OPEs in a solvable model, the 2-dimensional $O(N)$ nonlinear $\sigma$ model, by the dual space approach. A nontrivial mechanism of renormalon cancellation in this model can be understood from an integration identity on which the EBR method is founded. We demonstrate that the dual space approach can be useful by a simulation study imitating the QCD case. Application of this method to QCD calculations is also discussed.


Introduction
After the discovery of the Higgs boson, particle physics entered the era of the high-precision LHC experiments and the Belle II experiment, in which it is indispensable to improve accuracy of the theoretical prediction of QCD, to meet demands for the precision frontier physics.Due to the asymptotic freedom of QCD, theoretical predictions for high-scale observables can be made more precise by higher-order perturbative calculations.On the other hand, for systems with a scale of about O(1 − 10) GeV, theoretical uncertainties caused by "renormalons" can limit the accuracy of perturbative calculations seriously.
Renormalon [1][2][3][4] is a notion originating from the infrared (IR) region of certain loop diagrams, which is known to cause perturbative coefficients to diverge factorially.It makes perturbative expansion an asymptotic series, indicating that it is impossible to calculate the true value of an observable by perturbative expansion alone.Renormalons limit achievable accuracy of a perturbative calculation to be order (Λ QCD /Q) n with an integer n, for an observable characterized by a typical energy scale Q( Λ QCD ∼ 300 MeV).For a system of the electroweak scale, Λ QCD /Q ∼ 10 −3 would be negligible at present, while for a system of the bottom or charm quark, Λ QCD /Q ∼ 0.1 jeopardizes the precision of the prediction considerably.In high-precision flavor physics, it is desirable to remove the uncertainty due to renormalons from perturbative calculations.
Today, it has become a standard prescription to use a short-distance quark mass (e.g.MS mass) in various perturbative QCD calculations, and this is a way to remove renormalons.For instance, it was known that the perturbative expansion of the static QCD potential shows a drastically divergent behavior due to O(Λ QCD ) renormalon [5].When the QCD potential is combined with twice of the quark pole mass, the sum constitutes the main part of the total energy of a heavy quarkonium.The perturbative series of the pole mass, expressed in terms of a short-distance mass, also shows a divergent behavior due to O(Λ QCD ) renormalon [6,7].These divergent behaviors cancel out in the combination as the total energy [8][9][10].Besides heavy quarkonium observables, the O(Λ QCD ) renormalon in the B meson (semileptonic) inclusive decay width is also known to cancel out by rewriting the pole mass by a short-distance mass [7,11,12].The cancellation of renormalons has improved convergence of the perturbative series significantly.These aspects have been applied successfully in accurate determinations of fundamental physical constants, such as the masses of the heavy quarks [13][14][15][16][17], some of the Cabibbo-Kobayashi-Maskawa matrix elements [18][19][20], and the strong coupling constant α s [21].
In order to improve accuracy of perturbative calculations further, it is necessary to eliminate renormalons beyond O(Λ QCD ), which requires a theoretical framework that systematically incorporates nonperturbative QCD effects.The operator product expansion [22] (OPE) is suited for this purpose.The OPE of an observable (with scale Q) is given by an expansion in 1/Q, and the expansion coefficients are factorized into Wilson coefficients and nonperturbative matrix elements.Each Wilson coefficient can be calculated as perturbative series, while each matrix element needs to be evaluated nonperturbatively and has the size of an integer power of Λ QCD .Thus, the OPE provides a systematic double expansion in Λ QCD /Q and α s (Q 2 ).Since, however, renormalons in the perturbative series of the Wilson coefficients induce order (Λ QCD /Q) n uncertainties, they obscure systematic improvement of theoretical accuracy by the double expansion.It is believed that the uncertainties caused by renormalons cancel out within the framework of OPE, where each matrix element is believed to cancel (absorb) the corresponding renormalon in the Wilson coefficients [23].Hence, by subtracting renormalons from the Wilson coefficients and absorbing them into the matrix elements, each term of the OPE can be made well defined.
There have appeared analyses of the OPEs including cancellation of renormalons beyond the O(Λ QCD ) renormalon of the pole mass.In the OPE of the QCD potential, the O(Λ QCD ) and O(Λ 3 QCD ) renormalons in the leading-order (LO) Wilson coefficient were sub-tracted and this OPE was used to determine α s (M 2 Z ) precisely [24][25][26].Also the O(Λ 4 QCD ) renormalon contained in the lattice plaquette action was subtracted and absorbed into the local gluon condensate [27].Many different methods to achieve subtraction of renormalons from the LO Wilson coefficient have been developed recently [26][27][28][29][30][31][32][33][34][35][36].Among them, the FTRS method [32][33][34] utilizes the general property that IR renormalons of the LO Wilson coefficient can be suppressed by taking the Fourier transform of the Wilson coefficient.
In this paper we propose another renormalon subtraction method by replacing the Fourier transform in the FTRS method by the Laplace transform.In this way we can improve convergence property compared to the FTRS method without losing its merits.We call this method "dual space renormalon subtraction" (DSRS) method, by regarding the Laplace transform and its inverse transform as transformation between dual spaces (in analogy to the Fourier transform as transformation between the position and momentum spaces).
Subtraction of renormalons is generally scheme dependent.The DSRS and FTRS methods give the results in the same scheme as in most other methods, the principal value (PV) scheme [26][27][28][29][30][31].These methods provide one-parameter integral representations of the Wilson coefficient, which are equivalent to the conventional formula for the PV scheme expressed by the Borel integral.A major characteristics, as compared to other methods, is the ability to simultaneously separate the contributions of multiple renormalons, i.e., effects of different powers in Λ QCD /Q, by a simple procedure.In principle, the other methods can also separate multiple renormalon contributions, but the separation of renormalons corresponding to higher powers of Λ QCD /Q requires multiple steps of calculation to estimate the magnitudes of the renormalons, which complicates the calculation procedure.Another unique property of the DSRS and FTRS methods is that the Landau singularity of the running coupling constant, which is an unphysical singularity, is regularized together with subtraction of renormalons.It is a consequence of contour deformation in the dual space to avoid the Landau singularity.This property stabilizes the calculation of the OPE.
In this paper we also elucidate how the renormalon cancellation takes place in the entire OPE systematically.To our knowledge this is a first attempt in this direction for a general observable.(Most previous studies of renormalon subtraction focused on the analysis of the LO Wilson coefficient.)For this purpose the dual-space integral representation of an observable provides a useful framework.We can apply a generalization of the expansion-byregions (EBR) method [37][38][39] to this integral and investigate the mechanism of renormalon cancellation.
The EBR method is a technique for evaluating an integral by expanding the integrand in small parameters.The method is widely used in evaluation of Feynman loop integrals or in evaluation of Wilson coefficients (matching coefficients) in various effective field theories (EFTs).In this method, first an integral is regularized by analytic continuation (dimensional regularization), then the integral is separated into contributions from different regions of the integration variables.The correct expansion of the original integral is obtained in a sophisticated way, order by order in powers of a small parameter (ratio of external mass scales).
To investigate the properties of the DSRS method in detail, we apply the method to a solvable model, the O(N ) nonlinear σ model in two dimensions in the large N limit.It can be regarded as a toy model which mimics QCD.In this model, various observables can be computed exactly, and their OPEs can also be computed.We show in some examples that the dual transform removes all the IR renormalons as well as IR divergences of the Wilson coefficients.This leads to a systematic understanding of renormalon cancellation.
We demonstrate by a simulation study imitating the QCD case that the DSRS method can be used to determine nonperturbative parameters accurately.Finally we discuss how to apply the DSRS method to the QCD case.We discuss possible complications as well as solutions to them.Nevertheless, by construction of the DSRS method, we expect the application to be mostly straightforward.
The paper is organized as follows.In Sec. 2 we explain our basic idea how we understand the OPE and cancellation of renormalons.In Sec. 3 we explain the framework of the dual space approach.In Sec. 4 we apply the DSRS method to the O(N ) nonlinear σ model and study the OPE of some observables.In Sec. 5 we discuss application of the DSRS method to QCD and a remaining problem.Sec.6 is devoted to summary and conclusions.Details are explained in Appendices.In App.A we show the relation between the DSRS and FTRS methods.In App.B we explain how to resum UV renormalons.The proof of the EBR method in the dual space approach is given in App. C. Details of the calculations for the O(N ) nonlinear σ model are given in App.D.

Basic idea
To begin with, we explain our basic idea on how renormalons appear and cancel in the OPE, taking the static QCD potential as an example.
Consider the static QCD potential defined nonperturbatively, e.g., as computed by lattice QCD.It is well defined nonperturbatively both in the position space and momentum space.V (r) is related to V (q) by Fourier transform as We focus on the short-distance region 1/r Λ QCD , while we neglect effects from the ultrasoft scale for simplicity.We divide Eq. (2.1) in an EBR manner by employing dimensional regularization with d = 3 − 2 .
On the right-hand side (RHS), we divide V (r) into the contribution from the hard region, q = | q| > ∼ 1/r, and that from the soft region, q < ∼ Λ QCD 1/r.These contributions are denoted as V hard (r) and V soft (r), respectively.In the hard region, since q Λ QCD , we replace V (q) by its high-energy expansion, namely by the OPE, [ V (q)] OPE .In the soft region, instead we expand the factor e i q• r in q ( 1/r), while V (q) should be evaluated nonperturbatively.We obtain a series expansion in r, V soft (r) = n d n r n , with the expansion coefficients d n defined nonperturbatively.As is the case of the usual EBR of a Feynman integral, we expect that V hard (r) and V soft (r) are given by expansions in r (with log r corrections in the former).We also expect that the expansion coefficients in each of V hard (r) and V soft (r) contain IR and UV divergences as 1/ poles, while in the sum V hard (r) + V soft (r) the divergences cancel since V (r) is finite as → 0. In this way, the EBR gives the OPE of V (r) in r.
Similarly to the appearance of 1/ poles, there also appear imaginary part in the form ±iN n (Λ MS r) n /r in each of V hard (r) and V soft (r) for n = 1, 3, 5, . . . .The imaginary part in the Wilson coefficients of the OPE originate from renormalons.(This feature is based on the large-β 0 approximation, and for simplicity we neglect corrections to this approximation in this section. 1) The imaginary part should cancel in the sum V hard (r)+V soft (r) since V (r) is real. 2n short, the imaginary part of V hard (r) appears as follows.IR renormalons in the leading order (LO) Wilson coefficient of [ V (q)] OPE are suppressed, as shown in ref. [40].Therefore, the Wilson coefficient in q space can be calculated accurately in series expansion in α s (q 2 ).Then the IR renormalons in the regularized LO Wilson coefficient of V hard (r) can be separated as imaginary part, using the contour integral surrounding the Landau singularity of α s (q 2 ) in the complex q-plane (FTRS method) [32][33][34].The renormalon subtraction in this way agrees with the conventional PV prescription in the Borel resummation method in the all-order limit of the perturbative expansion.
The higher-order operators and associated Wilson coefficients in the OPE of V (r) or V (q) can be obtained using the potential nonrelativistic QCD EFT [48,49] by matching to full QCD.Hence, we expect that we can apply the same technique also to the subleading Wilson coefficients in [ V (q)] OPE (although we need to solve problems caused by the ultrasoft scale).
Cancellation of the imaginary part is a general feature which is guaranteed by validity of the generalized EBR method.However, explicit confirmation of the cancellation is difficult in the case of QCD, since V soft (r) is difficult to calculate in dimensional regularization.The cancellation can be checked explicitly in a solvable model such as the 2D O(N ) nonlinear σ model.
The EBR method used here and throughout this paper is different from the ordinary EBR method used for expanding loop integrals.We note that we use the EBR to the most external momentum integral only, and we do not apply EBR to the construction of V (q).We use the terminology "hard" and "soft" with respect to this most external momentum in Eq. (2.2).We will give a detailed explanation in the following sections.
In the following sections we will apply the above idea to general observables using the DSRS method.We will replace the Fourier transform by the dual transform, by which we separate renormalons more efficiently.We will explain how the imaginary part arise in the conventional approach and also in the dual space approach.

Framework of dual space approach
In this section we explain the framework of the dual space approach based on the Laplace transform.In Sec.3.1, we review the conventional treatment of IR renormalons and the OPE.In Sec.3.2 we introduce the dual space in which IR renormalons are suppressed.We use it to derive an expression of the renormalon-subtracted LO Wilson coefficient.In Sec.3.3 we construct the OPE of a general observable in the dual space approach and explain how renormalons emerge and are canceled.

OPE, IR renormalons and regularized Wilson coefficient
We consider a dimensionless observable S(Q 2 ) which is renormalization group (RG) invariant and characterized by a single hard scale Q ( Λ QCD ).Its OPE is expressed as where C 0 , C 1 are the Wilson coefficients and O 1 is the matrix element of an operator O 1 of mass dimension d 1 .The Wilson coefficients can be computed by perturbative calculation, and in particular, we express the LO Wilson coefficient as where α s (µ 2 ) is the coupling constant renormalized at the renormalization scale µ; the coefficients In the second equality we set µ = Q using the RG invariance.
In asymptotic-free theories such as QCD, naive perturbative series has IR renormalons, namely, Eq. (3.2) is a divergent series where the coefficient grows factorially, c n (0) ∼ n!.The conventional method to regularize such a divergent series is to use the Borel resummation method.The Borel transform of Eq. (3.2) is defined as  by 1/n!, this series in u has a finite radius of convergence, and the IR renormalons are seen as singularities on the positive real axis in the complex u-plane. 4In the leadinglogarithmic (LL) or large-β 0 approximation, each renormalon singularity is a pole, while including subleading logarithmic corrections, it becomes a branch point.We denote the positions of the IR renormalon singularities as u = u * (> 0), which make up an infinite sequence of integers or half-integers that depends on the observable.In particular, the most dominant renormalon, which produces the fastest divergent growth, corresponds to the singularity closest to u = 0.The regularized LO Wilson coefficient is constructed as where 0 ± i0 du. 5 The IR renormalon singularities are avoided by deforming the integration contour infinitesimally above/below the positive real axis.
Since C 0 (Q 2 ) Borel ± are complex conjugate to each other, their difference 2iδC 0 ≡ Borel − is pure imaginary.The contribution of each renormalon to δC can be calculated in terms of the residue of the renormalon or the integral around the branch cut: C * is shown in Fig. 1.In the second line, the form of δC u * 0 is determined assuming that the imaginary part is canceled by the corresponding term of OPE, Eq. (3.1) (in the case that there is a unique operator with dimension 2u * ) [23].Apart from the overall normalization N u * , the parameters b 0 , ν u * , γ 0 and the expansion coefficients s n can be calculated in perturbation theory; Λ MS denotes the integration constant of the running coupling constant in the MS scheme.
We neglect effects of singularities other than IR renormalons.For example, instanton effects are also known to generate singularities on the positive real axis (located far from the origin compared to the LO IR renormalon), but we neglect their effects.
For a convergent series without IR renormalons, Eq. (3.4) is simply the inverse Borel transform.
Thus, renormalons generate imaginary part to the regularized Wilson coefficient, and its sign depends on how the renormalon singularities are avoided along the integral contour in the u plane: where PV stands for the principal value integration or the average of the C ± integrals.
There is scheme dependence in how to subtract renormalons from the Wilson coefficient.A conventional prescription, the "principal value (PV) prescription," is to define a renormalon-subtracted Wilson coefficient by [C 0 ] PV .By definition, the renormalon contributions are minimally subtracted from C 0 (Q 2 ).

Dual transform of Wilson coefficient and renormalon suppression
Next we present an alternative expression of the regularized Wilson coefficient using a dual transform.Our argument is based on the formula of the inverse Laplace transform, We introduce two parameters (a, u ) and define the LO Wilson coefficient in the dual space as C0 (p 2 ) = 1 2πi It transforms the α s expansion in Q 2 space Eq.(3.2) to the α s expansion in a dual space (p 2 space).This is regarded as a dual transformation, similarly to the Fourier transform as a transformation from the coordinate space to the momentum space.With a proper choice of (a, u ), the factorial growth of the series due to IR renormalons can be suppressed in the dual space.This can be seen as follows.We can apply the (regularized) Borel resummation to both sides of Eq. (3.8) with respect to the series expansions in α s (µ 2 ).C 0 on the RHS is replaced by C 0 Borel ± , and the left-hand side (LHS) is replaced by C0 Borel ± .Then we take the difference of the ± contributions and express the RHS as the sum of the contributions from all the singularities u * .Thus, by replacing are of the form (Λ 2 MS /Q 2 ) n/a−u .Thus, IR renormalons of C0 (p 2 ) can be suppressed.It is possible to suppress more than one IR renormalons simultaneously with an appropriate choice of (a, u ) [32][33][34].In principle, it is also possible to take into account the higher-order α s (Q 2 ) corrections to the simple power of Λ 2 MS /Q 2 and anomalous dimension effects order by order.(See the end of this subsection.)Hence, the perturbative expansion of C0 can be regarded as a convergent series.
While C0 PT is a convergent series, C 0 PT is not, due to the IR renormalons.The difference emerges as follows.
where, using the RG invariance, we have rewritten Eq. (3.2) as with Integrating over t we have where cn (0)'s are computed by comparing α s (µ 2 ) expansion (using the expansion in Ĥ) of the first and second lines.We can see that the coefficients of the α s (µ 2 ) expansion of Eq. (3.14) have log n (µ 2 /p 2a ) dependence.The Laplace integral of such logarithms in the region p 2a < ∼ Λ 2 MS gives a factorial behavior ∼ n!, which is the origin of the IR renormalons in C 0 .Noting that C0 is RG invariant, we can resum the logarithms by setting µ 2 = p 2a in Eq. (3.14): As a result, the higher power logarithms are turned into the behavior of the running coupling constant α s (p 2a ) around the Landau singularity, which becomes an alternative source of renormalons upon considering the Laplace integral.The expansion of C0 (p 2 ) in α s (p 2a ) is also a convergent series for a sufficiently small α s (p 2a ).In practical calculations, the original series Eq. (3.12) is known only up to a certain order n ≤ k.Then we compute cn (0) also up to n = k.Thus, using the running coupling constant of the N k LL approximation and truncating C0 (p 2 ) at order α s (p 2a ) k+1 , one can compute a good approximation (at the level of N k LL) of the true C0 (p 2 ), which is defined by summing up the series to all orders.
The original Wilson coefficient C 0 (Q 2 ) can be reconstructed by the Laplace transform of C0 (p 2 ).As a result of the resummation described above, the singularity of the running coupling constant arises in the integral path along the real axis of the Laplace integral.Choosing the same integration contour C ± as Eq.(3.4), the regularized Wilson coefficient C 0 (Q 2 ) ± is given by where PV stands for the real part of the integration.It can be proven that C 0 (Q 2 ) ± given here agrees with that defined via the Borel resummation in the previous subsection.(See below.)Keeping the integral value, the integral path C ∓ can be deformed away from the Landau singularity in the right-half plane, such that α s (p 2a ) is sufficiently small along the entire integral path. 6Hence, we obtain a good estimate of C 0 (Q 2 ) PV and δC 0 (Q 2 ) using the series expansion of C0 (p 2 ).In particular, we call this method of computing C 0 (Q 2 ) PV dual space renormalon subtraction (DSRS) method.
As compared to the FTRS method [32][33][34], the DSRS method is more efficient.The Fourier transform in the FTRS method generates artificial UV renormalons, which worsen the convergence of the series in the Fourier space.In contrast the dual transform suppresses IR renormalons without generating artificial UV renormalons.It turns out that the DSRS method is equivalent to the FTRS method after resummation of the artificial UV renormalons in the latter.We give the proof in App. A. Once this is understood, we can prove the equivalence of C 0 (Q 2 ) ± calculated by the Borel resummation and by the DSRS method through the equivalence of the former and that by the FTRS method which is shown in App.C of ref. [33].Similarly, we can incorporate the corrections to the integer power forms proportional to (Λ QCD /Q) n of the imaginary part, given by the anomalous dimension and expansion in α s (Q 2 ), into the DSRS method.This can be achieved by using the corresponding formula for the FTRS method, given in App.B of ref. [33].(See, however, the discussion at the end of Sec.5.)

OPE in dual space approach
We explain how the OPE of an observable is derived in the dual space approach.Consider an observable which is calculated nonperturbatively (using lattice calculation, for example).Since it does not include perturbative series, it includes no renormalon ambiguities.Hence, its dual transform and inverse transform can be calculated without regularization of the integrals.
Below we assume that the above integrals are well defined.
In the case Q Λ QCD , we obtain the OPE of S(Q 2 ) by evaluating the p 2 integral in Eq. (3.19) by expansion in 1/Q using a generalized EBR method.We divide the integral into the contribution from the hard region, p 2 > ∼ t −1 = Q 2/a , and that from the soft region, p 2 t −1 , after regularizing the integral.
In the hard region, S(p 2 ) is replaced by its OPE, namely expansion in 1/p 2 and α s (p 2a ).In the soft region, e −tp 2 is expanded in p 2 .The IR and UV divergences which originate from the expansions are regularized by introducing δ analogously to the dimensional regularization, while the Landau singularity included in S(p 2 ) OPE is regularized by deforming the integral contour to C ∓ .(Or, we can deform the contour further away from the Landau singularity without changing the integral values.)The same contour C ∓ is chosen in both hard and soft contributions.Each term of S(p 2 ) OPE = i Ci (p 2 ) O i /p n i is uniquely determined by the dual transform of each term of S(Q 2 ) OPE given by Eq. (3.1). 7Hence, the Wilson coefficients Ci (p 2 ) can be calculated in perturbative expansion following the same calculational procedure as in the previous subsection, and there are no IR renormalons in the perturbative series.In this derivation of S(Q 2 ) OPE , we evaluate Ci (p 2 ) in expansion in α s (p 2a ), which generates the Landau singularity in the integrand.Convergence of the first term of Eq. (3.21) follows from the fact that the integral contour can be deformed, such that it is sufficiently away from the Landau singularity and the integrand converges as we include higher order terms of perturbative expansion in α s (p 2a ).
Since the original integral in Eq. (3.19) has no singularities and is well defined, the IR and UV divergences (regularized by δ) as well as imaginary part must cancel in the sum of the hard and soft contributions.This property is guaranteed by the validity of the generalized EBR method.We will prove the validity in the cases of the O(N ) nonlinear σ model and QCD (in N k LL approximation) in the following sections.We will see that the sign of the imaginary part of δ must be correlated with the choice of the regularized contour C ∓ in order to guarantee the validity of the EBR method.
The evaluation of the soft contribution, the second term of Eq. (3.21), involves nonperturbative evaluation of S(p 2 ), which is possible only in limited cases.We will confirm desired properties for some examples in the O(N ) nonlinear σ model.On the other hand, we can always reduce the nonperturbative information to a few nonperturbative parameters (matrix elements of local operators) if we use only the first few terms of the OPE.This aspect is useful in practical applications of the OPE for precise theoretical predictions of observables.

OPEs in nonlinear σ model
In this section we demonstrste usefulness of the dual space approach through a study of the OPEs in the 2D O(N ) nonlinear σ model.In Sec.4.1 we briefly explain the model and fix notation.In Secs.4.2 and 4.3 we analyze the OPE of a simple observable in detail.We confirm emergence and cancellation of renormalons in the dual space approach (Sec.4.2).It is followed by a simulation study, where we see that accuracy of theoretical prediction can be improved by the DSRS method (Sec.4.3).In Sec.4.4 we investigate the OPE of the self-energy of σ as another example.

The model
We consider the O(N ) nonlinear σ model in two dimensions in the large-N limit [41].The partition function of the model is given by where g 0 is the bare coupling constant.The Lagrange multiplier α(x) is introduced to take into account the constraint N a=1 σ a σ a = N .This model is solvable: it is an asymptotically free theory and a mass gap is generated nonperturbatively.(The O(N ) symmetry is unbroken.)The mass gap is given by the vacuum expectation value of α as where g(µ 2 ) denotes the coupling constant in the MS scheme renormalized at the renormalization scale µ.For convenience, we use the rescaled coupling constant below: We investigate the OPEs of observables in this model using the dual space approach.

Local condensate of δα regularized by gradient flow
The first example is a vacuum expectation value of a composite operator of δα(x) regularized by gradient flow, where the flow time t has mass dimension −2.I(t) is dimensionless, UV and IR finite (p 2 integral is well defined), and RG invariant. 8In this example, we give its dual transform by with the choice (a, u ) = (1, −2).Recall Eq. (3.20).
The OPE of Ĩ(p 2 ) (the expansion in m 2 /p 2 ) is given by IR renormalons are absent in all the Wilson coefficients in the dual space because the perturbative series for the Wilson coefficients terminate at finite orders when expressed with ḡ(p 2 ) [44].
In comparison, the OPE of I(t) has renormalons.The perturbative expansion in ḡ(µ 2 ) of the LO Wilson coefficient of I(t) is given by the Laplace transform of that of Ĩ(p 2 ), and its Borel transform is expressed by 9 This shows that the IR renomalons are located at u = 2, 3, 4, • • • and the leading renormalon uncertainty is O(t 2 m 4 ).There are no UV renormalons due to the gradient-flow regularization.These IR renormalons stem from the p 2 -integral of log(µ 2 /p 2 ) n encoded in ḡ(p 2 ).They should eventually cancel since I(t) is defined unambiguously.
In the dual space approach, we apply a dual transform to an observable under consideration whose perturbative expansion has IR renormalons.We therefore suppose that I(t) is an observable of our interest and we apply the dual transform as Eq.(4.6) (although 8 The OPE of a conceptually similar observable, δα(x) 2 with a UV cutoff regularization, has been studied in refs.[41][42][43].The gradient flow allows us to construct, in a simple and more sophisticated way, a UV and IR finite quantity which has IR renormalons in its perturbative evaluation.A similar observable in a nonlinear σ model has been considered in ref. [47] in the gradient-flow regularization, which discusses the renormalon structure of the CP N theory on R × S 1 .
the above procedure may look reversed).The parameters (a, u ) = (1, −2) in Eq. (4.6) are chosen to suppress the renormalons at u = 2, 3, 4, • • • .It is notable that this is an example where all the Wilson coefficients in the dual-space OPE are free of renormalons, which is the situation assumed in Sec.3.3.(This feature may seem obvious in this example, and we consider another example in Sec.4.4 to compensate for this drawback.)Since D α (p2 ) is known exactly even in the IR region, we can explicitly examine how renormalon cancellation takes place in the dual space approach and can check the validity of this framework in detail, as discussed in this and the following subsections.

Expansion-by-regions and renormalon cancellation
The OPE is obtained by EBR of the p 2 -integral in Eq. (4.7): where The EBR realizes the cancellation of IR renormalons.IR and UV divergences also appear, which are regularized by δ and canceled in I(t) OPE .
For a while we focus on the first few terms of the OPE of I(t): Explicitly they are given by and (s) 2 = 1.The Wilson coefficients are denoted as C (h,s) i according to their origins.We note that we applied EBR to the p 2 -integral of the Laplace transform and not to loop integrals associated with Feynman diagrams.That is why we have a soft-origin Wilson coefficient Each term of the OPE (4.14) has ambiguities, which stem from the regions where p 2 -integrals have divergences.The regions of p 2 which cause ambiguities are as follows: The hard quantities have ambiguities associated with the Landau pole singularity (p 2 ∼ m 2 ) and IR divergences (p 2 ∼ 0), whereas the soft quantities (or condensate) have UV divergences (p 2 ∼ ∞).In other words, the ambiguities originate from the region where their original expansions are not valid.By focusing on the contributions up to O(t 2 m 4 ) in the OPE, these imaginary part and divergences are given by where 0 < p 2 1 < m 2 < p 2 2 < 1/t.In the above equations, the dots "• • • " represent real finite contributions.In evaluating the integrals for C (h) 0,2 around the Landau pole we expand the exponential factor e −tp 2  1.The ambiguities found in C (h) 0 | Landau are the same as the renormalon uncertainties derived from Eq. (4.9).It turns out that C (h) 1 | Landau has higher order imaginary part which we neglect here.Also the higher order terms in ḡ(p 2 ) in the integrand for C (h) 2 | Landau give higher order imaginary part.We note that, in evaluating the integral for δα 2 in the region p 2 ∼ ∞, it is justified to expand Ĩ(p 2 ) in m 2 /p 2 as in Eq. (4.8).We therefore have the same integrands in the condensate (4.23) as those in the Wilson coefficients (4.20) and (4.22).In the integrand of Eq. (4.23), we only show the terms which give imaginary part or UV divergences.The imaginary part in Eq. (4.23) arises as [4] ∞ and the appearance of the log δ ± can be understood by a similar calculation to Eq. ( 4 The above cancellation indicates the validity of the EBR method.In fact the validity of the EBR relation Eqs.(4.10)-(4.12)can be proved following the strategy of ref. [39].(See App.C.) The crucial identity which guarantees the validity is given by10 where The cancellation explicitly confirmed above is a consequence of the particular identities which lead to ±iπm 4 ∓ iπm 4 = 0 and ∓2iπ − 2 log (−δ ± ) + 2 log δ ± = 0, respectively.Eq. (4.25) also explains the simultaneous absence of the imaginary part associated with the Landau pole 2 and the UV divergence in the condensate, which could be caused by either of them is zero, the other should also be zero.
Using the relation (4.25), we can discuss in a general manner the cancellation of the imaginary part in I ± h (t) and I ± s (t) which are related to renormalons.As seen explicitly in the first few terms of the OPE, the imaginary part arise from three regions: (i) p 2 → 0, (ii) p 2 ∼ m 2 , (iii) p 2 → ∞.The imaginary part from (i) and (ii) are included in the Wilson coefficients of I ± h , while that from (iii) is included in the matrix elements of I ± s . 11In fact, the above relation ensures the following identity: where n, k, ∈ Z ≥0 .n and k correspond to the expansion orders of Ĩ(p 2 ) in m 2 /p 2 and ḡ(p 2 ), respectively, i.e., the O(ḡ(p 2 ) k+1 /p 2(n−1) ) term of Ĩ(p 2 ), whereas is an expansion order of e −tp 2 in tp 2 .Note that in the regions where imaginary part or IR/UV divergences arise it is justified to expand Ĩ(p 2 ) in m 2 /p 2 and ḡ(p 2 ), and e −tp 2 in tp 2 ; for this reason it is sufficient to assume the form of the integrand of Eq. (4.29).The integrals on both sides are the contributions to O(t +2 ) term of I(t).For notational convenience, we have taken the < l a t e x i t s h a 1 _ b a s e 6 4 = " K m d a P n l j 0 k O 5 m C I A m o n P y c j e V a c = " Contour C L surrounding the Landau pole in the complex p 2 plane.
integral path on the positive real axis and instead shifted the Landau pole infinitesimally above or below the axis, m 2 ± ≡ m 2 ± i0.The imaginary part of the LHS of Eq. (4.29) is calculated as The first term is from the hard contribution (in I ± h ) and the second is from the soft contribution (in I ± s ). 12 Note that each contribution is independent of p 2 1 , p 2 2 , as it should be.The RHS of Eq. (4.29) contains the imaginary part due to the Landau pole of the running coupling constant ḡ(p 2 ), which is calculated as where C L is the closed path surrounding the Landau pole counterclockwise in the complex p 2 plane (Fig. 2).These imaginary part from different regions appear to be unrelated at a first glance, but the EBR relation shows that there is a one-to-one correspondence among them at each order term of the OPE and coupling expansion, namely before summing over n, k, .Using this relation, we can show that the imaginary part of I ± h (t) + I ± s (t) is identically zero, consistently with the fact that I(t) is real: where d n,k (0) represent perturbative coefficients appearing in Ĩ(p 2 ).See Eq. (4.39).The upper bound of the summation in k is equal to n because d n,k (0) = 0 for k > n.The above conclusion does not depend on these coefficients.This shows how the imaginary part cancel in the entire OPE of I(t).Although different reasonings justify the expansion in n, k, in the individual regions in the p 2 space (e.g., to extract the UV behavior of the soft expansion), essentially the factorization scale independence between the neighboring regions underlies the cancellation mechanism.
There is another nontrivial relation regarding the imaginary part.When summed over n and k, the imaginary part due to the Landau pole vanishes at each order of t.Namely, the coefficient of (m 2 t) vanishes as for arbitrary ∈ Z ≥0 .We used the explicit result for d n,k given below [Eq.(4.38)].The last equality is shown using an identity N n=0 N C n (−1) n n j = 0 for j < N .This feature can be understood as follows.If we take the contour C L sufficiently away from the Landau pole in the right half plane, the integrand becomes a convergent series everywhere along the contour C L in Eq. (4.31).Since before the OPE Ĩ(p 2 ) is regular and has no Landau pole inside the contour, the integral (residue) should converge to zero by summation over n and k.It is interesting to note that due to the identity (4.31) this argument shows that also the sum of the imaginary part from the other two regions vanishes.Therefore the absence of the Landau pole in Ĩ(p 2 ) guarantees the fact that I(t) is real.
We make a supplementary comment on the cancellation mechanism of the Landau pole singularity.The cancellation of the Landau pole, which is a physical requirement, is realized not within perturbation theory (the expansion in ḡ(p 2 )) but within the double expansion in ḡ(p 2 ) and m 2 /p 2 in the dual-space.The clear observation of the cancellation is allowed by the feature that the double expansion, in particular the perturbative expansion, is well defined in the dual space due to the absence of renormalons.If the perturbative expansion were asymptotic, as in the original t space, it would be difficult to examine the cancellation of the Landau pole in the double expansion.
We showed that the cancellation of log δ terms and of imaginary part take place in this example at any order of the OPE.This understanding leads us to a scheme appropriate for DSRS of the Wilson coefficients: (1) Obtain a converging series in the dual space, corresponding to the dual-space OPE.Perturbative coefficients to all orders in the OPE Below we present explicit formulas of the perturbative coefficients to all orders in the OPE.(See App.D for details of the calculations.)So far, we have treated the p 2 -integral with the integrands expressed with the running coupling constant ḡ(p 2 ).Here, we mainly study expressions with the fixed-scale coupling constant ḡ(µ 2 ).The OPE of Ĩ(p 2 ) is expressed in the form with the Wilson coefficients in the dual space given by (4.36)The Borel transform of [ Ĩ(p 2 )] OPE is given by The Borel transform of each Wilson coefficient is given by where we define an n-th degree polynomial of u as (a) j = Γ(a + j)/Γ(a) denotes the Pochhammer symbol.We confirm that B Ĩ (u) has no singularities on the positive real u axis.(In fact, it has no singularities in the entire u plane.)Hence, [ Ĩ(p 2 )] OPE does not have IR renormalons, even though it is rather obvious from the definition of I(t).We can also calculate d n,k (L p ) from Eq. (4.38).
The Borel transform of I ± h (t) is given by where the Borel transform of each Wilson coefficient is expressed in terms of X n (u) as We find that B C (h) As already stated, this is equivalent to Eq. (4.11).For example, we obtain the LO Wilson coefficient as from Eq. (4.41) for n = 0, as given by eq.(4.9).We can insert 1/ḡ(µ 2 ) = log(µ 2 /m 2 ), express Γ(2 − u) as C ∓ dx x 1−u e −x , transform the integral variable from x to p 2 = x/t, and integrate over u.In this way we derive the dual space expression

.44)
Although IR divergence appears in C (h) n ± for n ≥ 2 as multiple poles in δ at fixed order as mentioned above, by Borel resummation they turn to ∼ log δ.In the Borel integral, this originates from the u ∼ 0 region, where the integrand behaves as ∼ 1/(u + δ).In the dual space integral, the same log δ divergence stems from the p 2 ∼ 0 region.
The soft contribution can be expressed as by Eq. (4.12) and the definition of D α (p 2 ).All the matrix elements have UV divergences, δα −δ ± δα ± ∼ log δ, which stem from the region p 2 → ∞ and cancel the IR divergences in I ± h .Thus, The Wilson coefficients of the soft contribution are t independent and finite as δ → 0.

Thus, the Wilson coefficients of the hard contribution, [C (h)
n ] ± , have imaginary part whose sign depends on the regularization.In addition [C (h) n ] ± for n ≥ 2 contains an IR divergence which behaves as log δ.This is canceled by the UV divergence contained in the matrix element δα n−2−δ ± δα ± .The cancellation of log δ takes place in the same order of t n .On the other hand, the cancellation of the imaginary part is ensured at the level of the entire OPE, since renormalons shift the order counting of t.

DSRS method and fit of nonperturbative parameters (Simulation)
Using the OPE of I(t) obtained above, let us simulate practical application of the DSRS method to QCD.In that case, we can only calculate the Wilson coefficients of [I(t)] OPE using a finite number of perturbative coefficients.The nonperturbative matrix elements need to be determined by comparison of [I(t)] OPE with experimental data or with nonperturbative calculations by lattice QCD.As we showed in the previous section, the EBR method guarantees that the IR divergences and imaginary part in [C (h) n ] ± are canceled by the matrix elements.Therefore, it is possible to determine the well-defined nonperturbative matrix elements in a systematic approximation by subtracting the imaginary part and the IR divergences from the computation of the Wilson coefficients.
Specifically, we attempt to determine the nonperturbative parameters of the OPE by comparing the (simulated) experimental values with theoretical calculation of I h (t) which is defined by subtracting the IR divergence and imaginary part from I ± h (t) (Eq.(4.11)).The experimental value here is the exact I(t) defined by Eqs.(4.4) and (4.5), and there are no statistical uncertainties in this analysis.Only systematic uncertainties due to truncation of t expansion and perturbative evaluation of the Wilson coefficients are involved and studied.We perform an analysis of the OPE up to O(t 2 ), and I h (t) is given by

.51)
In the above formulas, we have subtracted the imaginary part and IR divergence from each Wilson coefficient by taking the principal value of the integral and calculating the explicit form of the IR divergence.The difference from I ± h (t) is given analytically by < l a t e x i t s h a 1 _ b a s e 6 4 = " X 8 A h Y 2 F b I l B l p z H 9 K X 7 X 7 e c K 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 g U f 5 o N e u q N z e q J b u q A H e v 2 2 V h j X i L z s 8 a x X t c L b b D v q W n j 5 V W X z L L H z r v r R s 0 Q R I 7 F X k 7 1 7 M R P d w q j q y / v H T w t j 8 7 1 h H / 2 j R / Z / S v d 0 w z d w y s / G / z k x f 4 I U f 4 D 6 + b m / g q W B v D q U H 5 w b z I 5 P J F / R i G 7 0 I M f v P Y x x T G E W B T 7 3 E J e 4 w r W S U k g Z V k a r q U p N o u n E h 1 A m 3 g C F 2 5 c r < / l a t e x i t > I(t) th up to O(t 1 ) < l a t e x i t s h a 1 _ b a s e 6 4 = " l z / r k + W F v o m w 9 h g P 4 x z p 1 X d g I which is RG invariant.The Wilson coefficients in the above formulas are expressed in expansion in ḡ(µ 2 ).In the analysis below, we will set µ 2 ∝ t −1 or µ 2 ∝ p 2 ; the coupling expansion of the integrand will be truncated at a finite order and the truncation error will be estimated.In order to calculate the IR divergence explicitly in terms of δ, C (h) 2 (t) ren PV is separated into two contributions by introducing a cutoff Λ IR , although the cutoff dependence cancels out between them.The IR divergence behaves as log δ, which is quite different from the typical form of IR divergences.The finite part that remains after subtracting the IR divergence is RG invariant.
The soft contribution I s (t) can also be computed in this model.Its expression up to O(t 2 ) is obtained as where Ei(z) = − ∞ −z dt e −t /t, and

.56)
To explicitly compute the UV divergence, we introduce a cutoff Λ UV , although δα 2 ren PV also has no overall cutoff dependence.Comparing Eqs.(4.52) and (4.53), we can see that the IR and UV divergences and the imaginary part cancel explicitly.δα 2 ren  PV can be evaluated by numerical integration to be As a cross check, we compare the exact (experimental) value of I(t) and I(t) th = I h (t) + δα 2 ren PV t 2 using α = m 2 and Eq.(4.57).The difference I(t) th − I(t) is shown in Fig. 3.We set µ 2 = p 2 in the integrand (in the dual space) and calculate I(t) th up to O(t 1 ) and O(t 2 ), which correspond to the black and red lines, respectively.The typical size of the coupling constant in this figure is We see that the exact value of I(t) can be predicted systematically, by increasing the expansion order of t for I ± h (t) and I ± s (t) and by the DSRS method.
We examine whether we can determine the nonperturbative parameters of the OPE correctly by a fit to an experimental data set.Here, we choose the position of the Landau pole, α = m 2 , and δα 2 ren PV given by Eq. ( 4.54), as the fit parameters.We consider the case where the coupling constant (or equivalently, the position of the Landau pole) needs to be determined by the fit as well.We compare a finite number of data points from theoretical computation of I h (t) and from exact (experimental) values of I(t).The experimental data is given by r defines a scale unit, which is equal to the ratio of m and the position of the Landau pole.(In this case r = 1.)On the other hand, I h (t) is obtained using the Wilson coefficients computed from a finite number of perturbative coefficients.We calculate the Wilson coefficients in two different ways, in order to check whether the dual space approach is more useful than the ordinary perturbative approach.The first fit function is calculated by the RG-improvement in the original t-space, which is given by 13 where n in ḡ(µ 2 ). 13In order to subtract the IR divergence in the same prescription as in the DSRS method, we proceed as follows.The IR divergence takes the form ∼ 1/(u + δ) in B C (h) 2 of Eq. (4.41), and the Borel integral of 1/(u + δ) gives the log δ singularity.Hence, we subtract the contribution from B C (h) 2 corresponding to log δ, expand in u and take the Borel integral.We adjust the finite (δ independent) term such that the same divergence is subtracted as in the DSRS method.
The second one is calculated using the DSRS method and given by  at k 1 = 3 and k 2 = 3, respectively, and fix the scale at s = 1 in both calculations, for simplicity. 14In the analysis with Eq. (4.61), we assign another systematic uncertainty from (supposedly) unknown O(t 3 ) nonperturbative corrections in the OPE.The uncertainty is estimated by the difference of the central values of the parameters r, A and B determined by the fits, using I PV fit and I PV fit + Cm 6 t 3 , in which C is an additional fit parameter.This effect is negligible for the analysis with Eq. (4.59) because the renormalon effect of order m 4 t 2 is large.
The results of the fits are shown in Fig. 4. The blue dots with error bars represent the parameters determined with different truncation orders k 0 .The red dashed lines show the exact values of the parameters.In the RG-improvement method, the accuracy seemingly improves up to k 0 ∼ 13 for all the parameters, but the uncertainty increases beyond that    < l a t e x i t s h a 1 _ b a s e 6 4 = " I s L H r T r C / v v W i 7 u r t z S m V P + / k < l a t e x i t s h a 1 _ b a s e 6 4 = " j X u 7 h q 6 e R a r 0 Q y a r Y q M a u q m n V d k R + i a I T K u 5 u o i b 9 l C r i q 6 y C m V z c 5 + r i 5 s R z O N P b d h i W J V L h v a o a b K L l P p S o l K k R g l y I t o L 0 j 6 I A Y / t s 3 I L f Z x A B M q a q h C w I D L W I c M h 1 s B S R A s 5 o p o M m c z 0 r x 9 g R M E W V v j L M E Z M r M V H s u 8 K v i s w e t O T c d T q 3 y K z t 1 m Z R R x e q I 7 a t M j 3 d M L f f x a q + n V 6 H h p 8 K x 0 t c I q h U 9 n 0 + / / q q o 8 u z j 6 U v 3 p 2 c U h V j 2 v G n u 3 P K Z z C 7 W r r x 9 f t N N r u / H m A l 3 T K / u / o h Y 9 8 A 2 M + p t 6 s y N 2 L x H k D 0 j + f O 5 e k F 1 K J J c T q Z 1 U b H 3 D / 4 p h z G E e i / z e K 1 j H F r a R 4 X P L O M M 5 L g L P U k i a l m a 6 q V L A 1 0 z h W 0 j R T 0 l P i p c = < / l a t e x i t > k < l a t e x i t s h a 1 _ b a s e 6 4 = " j X u 7 h q 6 e R a r 0 Q y a r Y q M a u q m n V d k R + i a I T K u 5 u o i b 9 l C r i q 6 y C m V z c 5 + r i 5 s R z O N P b d h i W J V L h v a o a b K L l P p S o l K k R g l y I t o L 0 j 6 I A Y / t s 3 I L f Z x A B M q a q h C w I D L W I c M h 1 s B S R A s 5 o p o M m c z 0 r x 9 g R M E W V v j L M E Z M r M V H s u 8 K v i s w e t O T c d T q 3 y K z t 1 m Z R R x e q I 7 a t M j 3 d M L f f x a q + n V 6 H h p 8 K x 0 t c I q h U 9 n 0 + / / q q o 8 u z j 6 U v 3 p 2 c U h V j 2 v G n u 3 P K Z z C 7 W r r x 9 f t N N r u / H m A l 3 T K / u / o h Y 9 8 A 2 M + p t 6 s y N 2 L x H k D 0 j + f O 5 e k F 1 K J J c T q Z 1 U b H 3 D / 4 p h z G E e i / z e K 1 j H F r a R 4 X P L O M M 5 L g L P U k i a l m a 6 q V L A 1 0 z h W 0 j R T 0 l P i p c = < / l a t e x i t > k 0 order.Since the u = 2 renormalon is included in C (h) 0 , we expect that the perturbative series converges up to k 0 ≈ k * 0 ≡ 2/ḡ ≈ 13.8, which is consistent with the result.While r and A are obtained with reasonable accuracy around k 0 = k * 0 even under the influence of the u = 2 renormalon, B is affected by an order 100 per cent uncertainty in the determination.This feature is also consistent with the effects of the u = 2 renormalon.The determined values of r, A and B at k 0 = 13 ≈ k * 0 are given by r − 1 = 4(6) × 10 −6 , A − 1 = 5(6) × 10 −3 , B = −3.1(2.9),(4.64) which shows the limit of accuracy of the parameter determination in the presence of the renormalon.In general, if we use data points from a smaller (larger) t region, we can determine r and A with better (worse) accuracy, while B is always subject to an order 100 per cent uncertainty.In the analysis using the DSRS method, the perturbative uncertainties decrease even beyond k 0 = k * 0 .This is because all the renormalons are removed from C (h) 0 .Moreover, the convergence is even accelerated before reaching the order k * 0 , where the original series is already dominated by the LO renormalon contribution.On the other hand, the total uncertainty never becomes smaller than a certain bound due to the systematic uncertainties by the higher-order of the OPE, namely the t expansion.The determined values of r, A and B at k 0 = 13 ≈ k * 0 are given by r − 1 = 4(4) × 10 −8 , A − 1 = 4(4) × 10 ).The results show that the DSRS method is useful to determine the nonperturbative parameters accurately by eliminating renormalons in the entire OPE.The theory prediction by the OPE can be made accurate using the determined parameters.

Self-energy of σ
The inverse propagator of σ with momentum q µ is given by (q 2 + m 2 ) at LO of the 1/N expansion.The second example we consider is the self-energy of σ at O(1/N ), denoted as Π σ (q 2 ).The self-energy is included universally in calculations of correlation functions of σ and is given by the diagram in Fig. 5. (We omit the tadpole diagram which is absorbed by renormalization of α .)Cancellation of the renormalons in the entire OPE of the selfenergy has been shown in ref. [45].Here, we show that the cancellation can be understood systematically in the dual space approach.The dual transform indeed eliminates all the IR renormalons.We also see that the DSRS method is useful to improve theoretical accuracy in the OPE.
In order to remove the UV divergences of Π σ (q 2 ), we consider the second derivative of the self-energy, It is defined to be dimensionless and RG invariant.The above integral is well defined nonperturbatively, without UV or IR divergence.We can construct the OPE of ρ(q 2 ) in the dual space approach similarly to the previous example.We define the dual space observable ρ(p 2 ) by < l a t e x i t s h a 1 _ b a s e 6 4 = " S z m 9 S A b U 9 W m Q a u s j N j J S s B a G R m s = " > A A A C b 3 i c h V H L S g M x F D 0 d 3 / X R q g s F Q Y p F 0 U 1 J p a i 4 E t 2 4 b N W q Y E U y 0 1 S H p j P D T F q o x R / w A 3 T h w g e I i J / h x h 9 w 0 Using the generalized EBR, we can express the OPE of ρ(q 2 ) as the sum of the hard and soft contributions: where the integrals are regularized in the same way as in the previous example.We can check that the Wilson coefficients of [ρ(p 2 )] OPE do not have IR renormalons or IR divergences.In fact, the Borel transform of [ρ(p 2 )] OPE is given by It can destabilize convergence of the series expansion in the dual space.Since the UV renormalons are Borel summable, we stabilize the series by Borel resummation.This is possible in the dual space approach, see App.B for details.B ρ has no pole at u = 0 and is independent of any regularization parameter.Hence, [ρ(p 2 )] OPE is also free from IR divergences.This aspect is the same as that of the previous example Ĩ(p 2 ).
which provides a framework to analyze properties of the OPEs including renormalons.Secondly, using the framework, we show a method to analyze cancellation of renormalons systematically at arbitrary order of the expansion in the OPE [using Eqs.(4.25) and (4.29)], for a general one-scale observable.Thirdly, our method goes beyond diagrammatic analyses of renormalons, since our analysis uses only the dependence on the most external single scale p 2 , without resolving detailed kinematical information of individual diagrams.

Application to QCD and discussion
In this section we discuss how to subtract multiple renormalons simultaneously from the OPE of an observable in QCD.We consider a procedure which can be used in practical applications.We also explain an unsolved question.
We first discuss the following points in addition to the procedure used for the O(N ) nonlinear σ model in Sec.4: (A) How to incorporate the effects of the fact that the form of each renormalon can deviate from a simple integer power of Λ QCD .(B) How to incorporate the effects of the beta function beyond one loop order (b 1 , b 2 , • • • = 0).
(A) Renormalons are generally not in the simple power form of Λ QCD /Q but include corrections in the form of the anomalous dimension and expansion in α s (Q 2 ), as given by Eq. (3.5).A method for incorporating these corrections is given in App.B of ref. [33] in the context of the FTRS method.Due to equivalence of the DSRS method to the FTRS method combined with the resummation of the artificial UV renormalons, we can also use that method in the dual space approach.
(B) Including b 1 , b 2 , • • • = 0 in the dual space approach is straightforward, and how this is done is already explained in Sec. 3 in principle.Note that Λ MS in e.g.Eq. (3.5) is defined by the k-loop beta function with k not restricted to be one.The points which may be worth notifying are as follows.(i) In practical calculations, it is convenient to numerically obtain the running coupling constant in the complex p 2 plane by solving the RG equation along the contour C ± .(ii) IR divergences occur from the integral over p 2 from the IR region p 2 ∼ 0, beyond a certain order in expansion in 1/p.By regularizing the p 2 integral similarly to the nonlinear σ model, the IR divergences appear as ∼ log δ .This is the case, even including NLL and beyond, since the coupling constant behaves as as p 2 → 0, and only the leading term causes the divergence.For instance we can subtract the log δ term as in Sec. 4.
Let us add some more comments on the properties of the dual space approach in QCD.
(1) The conventional picture of renormalon cancellation in the OPE is that the imaginary part of the (regularized) Wilson coefficients originating from the IR renormalons are canceled by the imaginary part contained in nonperturbative matrix elements of higher dimensional operators.In the dual space approach, the cancellation can be understood as follows.The cancellation takes place at each order of the (simultaneous) soft-and-hard expansion in the p 2 integral.In the integrand, the soft expansion generates the expansion first through the anomalous dimension and then order by order in α s (Q 2 ) expansion, is only an asymptotic approximation in this case.

Summary and conclusions
We proposed the DSRS method as a new method to separate and subtract IR renormalons from Wilson coefficients in the OPE of a general one-scale observable S(Q 2 ).We constructed a dual transform from S(Q 2 ) to S(p 2 ) such that the OPE of S(p 2 ) is free from renormalons.The DSRS method has the following properties.
(a) The DSRS method is equivalent to the conventional PV prescription for subtracting IR renormalons, which is based on the (regularized) Borel resummation technique.
(b) The DSRS method is equivalent to the FTRS method after resumming the artificially generated UV renormalons in the latter method.In the DSRS method (unlike the FTRS method) no UV renormalons are additionally generated by the dual transform.
(c) The DSRS method can remove multiple IR renormalons simultaneously.
(d) In the DSRS method it is not required to evaluate the normalization constants of renormalons.The real and imaginary part of the OPE are separately calculable by the given procedure.
(e) Wilson coefficients are calculated by integrating series expansions in α s (p 2a ) along the deformed contour in the p 2 plane.As a result, this method is insensitive to instabilities due to the Landau singularity of the running coupling constant.(Note that, with the current construction of the dual transform, this property is limited to the case where the renormalons are given by integer powers of Λ QCD /Q.) (f) The DSRS method can give the PV result approximately with a given finite number of perturbative coefficients.The accuracy improves as more perturbative coefficients are included.This feature is of use in practical applications to QCD.
(g) Using the DSRS method, we can determine nonperturbative matrix elements accurately by a fit to experimental data or lattice results.
(2) In the entire OPE, IR and UV divergences cancel between the contributions from the hard and soft regions.The imaginary part cancel among the contributions from the hard, soft and Landau singularity regions.The cancellation can be understood at each order of the (simultaneous) soft-and-hard expansion, irrespective of the details of the observable.
We checked the above features by applying the DSRS method to some observables in the 2D O(N ) nonlinear σ model in the large N limit.In addition to their detailed examinations, we find that IR divergences in the Wilson coefficients are removed in the dual space, together with IR renormalons.As a result, we can subtract the IR divergences in the original space in an RG invariant manner using the p 2 integral.
We discussed application of the DSRS method to QCD.In principle, the method can be applied straightforwardly.The equality that guarantees cancellation of renormalons is given by Eq. (5.2), which includes subleading logarithmic corrections of the running coupling constant.We note that the FTRS method has already been applied to several QCD observables and good agreement with theoretical expectation has been confirmed [33].We anticipate similar features for the DSRS method.
In conclusion, the DSRS method gives approximation of the OPE from a finite number of perturbative expansions of Wilson coefficients, with subtraction of multiple renormalons.Hence, the method is expected to be simple and efficient in phenomenological applications, for improving accuracy of theoretical calculations.
There remains an unsolved problem, which is not specific to the DSRS method, when the correction to the power behavior (Λ QCD /Q) n of each IR renormalon is incorporated by anomalous dimension and expansion in α s (Q 2 ).Such a correction can be an asymptotic series and we do not know how to resum the series systematically.This is a subdominant effect and at present it would be a less serious question phenomenologically.Nevertheless, it would be a subject that needs to be solved toward realization of high precision theoretical calculations in the future.
The regularized S(Q 2 ) is obtained by the inverse Fourier transform of S (τ ) with contour deformation.Hence, the principal value integral that eliminates renormalons are given by In the third equality the variable of integration is changed, and in the fourth equality we perform the y integration.Finally, we define τ = p, r 2 = t for clarity.This is just in the form of the Laplace transform from t-space to p 2 -space and is an expression equivalent to the DSRS method.See Eqs.(3.15) and (3.18).

B Resummation of UV renormalons in DSRS method
In this appendix, we give a resummation formula for UV renormalons in the case that they are given by simple poles in the Borel plane and in the LL approximation.The formula is used in the analysis of Sec.4.4.In the formula, we introduce a new parameter ū to control the resummation.Let us take the LO Wilson coefficient in the dual space C0 (p 2 ) and explain how to resum renormalons.C0 (p 2 ) can be expressed by the following one-parameter integral form, Then C 0 ± with resummation of UV renormalons is calculated by Recall that t = Q −2/a .In the third line we have changed the variable of integration from p 2 to p2 = p 2 ζ −1 .

C Proof of EBR method in Jantzen's manner
In this appendix we provide a proof of the EBR method which we use in this paper.We follow the strategy which was used by Jantzen in ref. [39].
Let us take the observable I(t) studied in Sec.4.2 as an example: where T x [f (x)] represents the Taylor expansion of f (x) in x.In the last line, we expand the integrand using hierarchy of the variables satisfied in each region.([ Ĩ(p 2 ) ± ] OPE represents double expansion in 1/p 2 and 1/ log(p 2 /m 2 ± ).)We further rewrite the integral regions:  (C.3) In the last expression, the first two terms are equal to I ± s (t) and I ± h (t), respectively.The last term is the sum of terms in the form of Eq. (4.25) or (5.2), which are zero by analytic continuation.This completes the proof.

D Calculations for OPEs in nonlinear σ model
In this appendix we present details of the calculation for the OPEs of the observables in the nonlinear σ model.

D.1 Factorization formula for Borel transform
First we present a formula which is useful in calculating Borel transform of a perturbative series in the LL approximation, namely in the case that the one-loop running is exact, such as in the nonlinear σ model.
Consider an observable given as an integral form S(q 2 ) = We assume both S(q 2 ) and G(p 2 ) are RG invariant.
We define an integral where G(p 2 ) is replaced by a simple power of x n ∞ k=0 d n,k (0) ḡ(p 2 ) k+1 , (D.9) To compute the dual space observable [ρ(p 2 )] OPE , the region (iii) does not contribute.This is because, it gives a series expansion with only integer powers of q 2 (without log q 2 terms), whose dual transformation vanishes.The sum of the contributions from the regions (i) and (ii) can be calculated together as follows.
ρ(q 2 ) (i,ii) reg = q 2 ∂ ∂q 2 where We set D = 2 − 2δ.Hence, using the factorization formula, we obtain the Borel transform of ρ ± h (q 2 ) as We can rewrite 2 F 1 (a, b; c; − q 2 m 2 ) in terms of 2 F 1 (a , b ; c ; − m 2 q 2 ) using a functional identity, by which one can readily find the all order formula of F (n+u−1) in expansion in m 2 /q 2 .All the IR renormalons and IR divergences of ρ ± h (q 2 ) originate from the factor Γ(2 − n − u − δ) in F (n + u − 1).We can take the limit δ → 0 for the other part, which gives the formula Eq. (4.74).
We can take the dual transform of Eq. (4.74) using Eq.(4.69), order by order in expansion in 1/q 2 .The dual transform essentially generates 1/Γ(2 − M − u − δ).Since Γ(2−n−u−δ)/Γ(2−M−u−δ) for n ≤ M is a polynomial of u, all the singularities corresponding to the IR renormalons and IR divergences are canceled by the dual transform.After the cancellation, we can safely take the limit δ → 0. Thus, we obtain B ρ given by Eq. (4.71).It is straightforward to obtain [ρ(p 2 )] OPE from B ρ at each order of the expansion in 1/p 2 and ḡ(p 2 ).We expect that the obtained result for B ρ(u) or [ρ(p 2 )] OPE is independent of the regularization scheme we employ at intermediate stage of the calculation, due to the absence of IR renormalons and IR divergences. 16  16 The algorithm for deriving [ρ(p 2 )]OPE is partly cut short, for simplicity.The true algorithm is as follows.
We first determine the regularized dual transform which cancels all the IR renormalons and IR divergences of ρ ± h (q 2 ), and afterwards we remove the regularization.Once [ρ(p 2 )]OPE is determined, ρ ± h (q 2 ) can be redetermined regularization dependently.The regularization dependence cancels against that of ρ ± s (q 2 ).

. 3 )
Here, b 0 denotes the one-loop coefficient of the beta function.(We adopt the convention where b 0 = (11 − 2n f /3)/(4π) in the case of QCD.) Due to the acceleration of convergence < l a t e x i t s h a 1 _ b a s e 6 4 = " Y u o 2 1 1 4 0

Figure 1 .
Figure 1.Contour C * for the imaginary part of the regularized Wilson coefficient in Eq. (3.5).

( 2 )
Apply the Laplace transform to the dual-space OPE, which provides the Wilson coefficients of the hard contribution.(3) In the Laplace transform, IR divergences are regularized by p −2δ ± and the log δ terms are subtracted minimally by counter terms.(4) Adapt the PV prescription, i.e., remove the imaginary part.(5) The nonperturbative matrix elements are treated as parameters.(If necessary, one can resum UV renormalons in addition, see Sec. 4.4.)This DSRS scheme naturally encodes the cancellation between IR and UV divergences and the renormalon cancellation.(It adds renormalization of IR divergences to the conventional PV prescription.) and so on.(For n < 2 we can safely set δ to zero.)Thus, I ± h (t) has IR renormalons.In addition, expanding B C (h) n in u, we see that the perturbative expansion of C (h) n in ḡ(µ 2 ) for n ≥ 2 includes IR divergences as multiple poles in δ.This originates from the u ∼ 0 region, by expansion of B C (h) n ∼ 1/(u + δ) in u.I ± h (t) is constructed by Borel resummation as 7 e b C A 3 e y c w 8 8 8 z 7 v P P M j G L p m i M Y a w a k j s 6 u 7 p 5 g b 6 r t J t u d z e 7 0 y b V N O 7 + g I M T S S M i 4 U e 4 + A M O f o I 4 k r g 4 e H e 7 i S B 4 J z P z z D P v 8 8 4 z M 4 q l a 4 5 g 7 D E g t b V 3 d H Y F u 0 M 9 v X 3 9 A + H B o S 3 H r N g q z 6 q m b t o 7 i u x w X T N 4 V m h C 5 z u W z e W y o v N t p b T s 7 m 9 X u e 1 o p r E p a h b f n M 8 a t H Z r O p 5 a p V N 0 6 j Y p I 5 h g D + y K v b B 7 d s 2 e 2 P u v t e p e D d d L j W a l p e V W f u B k Z O P t X 1 W Z Z o H i p + p P z w K H S H l e N f J u e Y x 7 C 7 W l r x 6 d v m w s r k / U Y + y C P Z P / c / b I 7 u g G R v V V b a 7 x 9 T O E 6 A M S 3 5 / 7 J 9 i a j S f m 4 8 m 1 Z D S 9 5 H 9 F E G M Y x y S 9 9 w L 0 5 V 6 Z 5 l B p L o v k a p r a v / 0 9 D Y l G p u a f 3 b l m 7 v W A r c k m + I g u F a r r + i a 4 G w T E c U p C k t s e L 5 Q r N 1 S y z r u 5 P R / n J Z + I H p O o t y z x P r t r b t m E X T 0 C R T m + n s d E 7 2 b 4 R y p 3 I Q r v l 2 p u Q d S D f C h m Z l Z i o 5 u a H 2 c x b l K Y 7 M V 6 A m I I s k Z t 3 0 G d a w B R c G S r A h 4 E A y t q A h 4 L Y K F Q S P u X W E z P m M z H h f o I I U a 0 u c J T h D Y 3 a X x 2 1 e r S a s w + u o Z h C r D T 7 F 4 u 6 z M

Figure 3 .
Figure3.Log-log plot of the difference between I(t) and I(t) th = I h (t) + δα 2 ren PV t 2 .The black (red) line corresponds to I(t) th computed up to O(t 1 ) (O(t 2 )).The slope of each line represents the power behavior of t.
.60) and A = α /m 2 and B = δα 2 ren PV /m 4 stand for the fit parameters.[C (h) n ] PT denotes the (naive) perturbative expansion of C (h) ) RG , which would give an order 100 per cent uncertainty in the determination accuracy of B. On the other hand, in Eq. (4.61), the Wilson coefficients are free of renormalons.Hence, it is expected that we can extract the exact values of the parameters,r = 1, A = 1, B = −1.386• • • ,(4.63)accurately by the fit.The details of the fits are as follows.The data points used for the fits correspond to t ≡ m 2 t = n/10 4 for n = 1, 2, • • • , 10.The maximal size of the coupling constant in this region is ḡ = 1/log(10 3 ) ≈ 0.14.The parameters r, A and B are determined simultaneously by the method of least squares.For each Wilson coefficient, only k 0 , the truncation order of C (h) 0 , is varied to examine the accuracy of the parameter determination.The systematic uncertainty is evaluated from the scale variation of C (h) 0 .The central value of the parameters is taken from the calculation with s = 1 in Eqs.(4.59) and (4.61), and the systematic uncertainty is assigned from the difference as we vary s ∈ {1/2, 1, 2}.In this analysis, the truncation errors of perturbative expansions of the Wilson coefficients C (h) nfor n ≥ 1 are assumed to be negligible.Hence, we fix the truncation orders of C

2 > 2 >
t ∈ [10 -4 ,10 -3 ], Δt=10 -4 t ∈ [10 -4 ,10 -3 ], Δt=10 -4 10 -4 ,10 -3 ], Δt=10 -4 < l a t e x i t s h a 1 _ b a s e 6 4 = " t H C g w f b 5 5 5 n 3 e e m Z E M V b F s x r o e Y W R 0 b H z C O + m b 8 k / P B I K z c 3 l L b 5 o y z 8 m 6 q p t F S b S 4 q m g 8 Z y u 2 y o u G y c W G 5 5 5 n 3 e e m Z E M V b F s x r o e Y W R 0 b H z C O + m b 8 k / P B I K z c 3 l L b 5 o y z 8 m 6 q p t F S b S 4 q m g 8 Z y u 2 y o u G y c W G r D A n t g d 6 7 B H d s + e 2 f u v t e p e j a 6 X M 5 qV n p Z b x d D F T P r t X 5 V O s 4 v T T 9 W f n l 2 U s O F 5 F e T d 8 p j u L d S e v n Z + 2 U l v p h b q i 6 z F X s h / k 7 X Z A 9 3 A q L 2 q N / s 8 d Y 0 g f U D 8 + 3 P / B I e r s f h a L L G f i C a 3 / K 8 Y x h z m s U T v v Y 4 k d r G H D J 1 r 4 w p N t A I d a V a K S P O 9 V C n g a 6 bx J a T l D + C w j v M = < / l a t e x i t > RG imp.
r D A n t g d 6 7 B H d s + e 2 f u v t e p e j a 6 X M 5 qV n p Z b x d D F T P r t X 5 V O s 4 v T T 9 W f n l 2 U s O F 5 F e T d 8 p j u L d S e v n Z + 2 U l v p h b q i 6 z F X s h / k 7 X Z A 9 3 A q L 2 q N / s 8 d Y 0 g f U D 8 + 3 P / B I e r s f h a L L G f i C a 3 / K 8 Y x h z m s U T v v Y 4 k d r G H D J 1 r 4 w p N t A I d a V a K S P O 9 V C n g a 6 bx J a T l D + C w j v M = < / l a t e x i t > RG imp.

Figure 4 .
Figure 4. Results of the determination of the nonperturbative parameters r, A, B by the fit.The left (right) panels correspond to the Wilson coefficients C (h) n calculated by the RG-improvement (DSRS method).k 0 (horizontal axis) represents the truncation order of the perturbative expansion of C (h) 0 .The blue dots with error bars represent the results of the fit by varying the scale from the central value (s = 1) by factors 1/2 and 2. The red dashed lines represent the exact values of the parameters.
.67) < l a t e x i t s h a 1 _ b a s e6 4 = " l Q R z 9 p o k F C V Y + 6 M K 3 m 1 h E q 1 d 7 R s = " > A A A C a X i c h V F N L w N B G H 6 6 v q o + W l y E S 2 N T c W q m I o i T c H H U V k u C y O 4 a N b p f 2 Z 0 2 o f E H n N w E J x I R 8 T N c / A G H / g R x J H F x 8 O 5 2 E6 H B O 5 m Z Z 5 5 5 n 3 e e m d F d U / i S s W Z M 6 e j s 6 u 6 J 9 y b 6 + g c G k 6 m h 4 b L v 1 D y D l w z H d L w N X f O 5 K W x e k k K a f M P 1 u G b p J l / X q 8 v B / n q d e 7 5 w 7 D V 5 6 P J t S 6 v Y Y k 8 Y 7 Y 6 / s k d 2 z Z / b x a 6 1 G W C P w c k i z 3 t J y d y d 5 M l p 8 / 1 d l 0 S y x / 6 X 6 0 7 P E H u Z D r 4 K 8 u y E T 3 M J o 6 e t H Z 6 / F h U K m M c m u 2 Q v 5 v 2 J N 9 k A 3 s O t v x k 2 e F y 6 R o A / I / X z u d l C e z u Z m s

/ l a t e x i t > ↵Figure 5 .
Figure 5. Self-energy diagram of σ at O(1/N ).The tadpole diagram is omitted.
Bm 4 t 2 , IR exactly.)Inboth formulas, we subtract the IR divergence of the O(t 2 ) term in the OPE, and each Wilson coefficient is calculated by truncation of the perturbative expansion.In Eq. (4.59), the Wilson coefficients contain renormalons which cause uncertainties in the fit.The major uncertainty comes from the u = 2 renormalon in C (19)B = −1.406(19),(4.65)whichare consistent with the exact values of the parameters, Eq. (4.63).We can see that the error sizes of r and A in Eq. (4.65) are orders of magnitude smaller than those in Eq. (4.64).This is because the major source of systematic uncertainties in Eq. (4.65) is the O(t 3 ) OPE corrections, while that in Eq. (4.64) originates from the u = 2 renormalon in C