aainstitutetext: Department of Physics, Tohoku University,
Sendai, 980-8578 Japanbbinstitutetext: Institute of Particle and Nuclear
studies, KEK,
Tsukuba, 305-0801, Japan

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

Yuuki Hayashi a Go Mishima a Yukinari Sumino b HiromasaTakaura
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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.

TU–1182

KEK–TH–2501

##  1 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 ${\cal O}(1-10)~{}{\rm GeV}$, theoretical
uncertainties caused by “renormalons” can limit the accuracy of perturbative
calculations seriously.

Renormalon Gross:1974jv ; Lautrup:1977hs ; `t Hooft:1977jv ; Beneke:1998ui  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 $(\Lambda_{\rm QCD}/Q)^{n}$ with an integer $n$, for
an observable characterized by a typical energy scale $Q(\gg\Lambda_{\rm
QCD}\sim 300~{}{\rm MeV})$. For a system of the electroweak scale,
$\Lambda_{\rm QCD}/Q\sim 10^{-3}$ would be negligible at present, while for a
system of the bottom or charm quark, $\Lambda_{\rm QCD}/Q\sim 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. $\overline{\rm 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 ${\cal O}(\Lambda_{\rm QCD})$ renormalon
Aglietti:1995tg . 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 ${\cal
O}(\Lambda_{\rm QCD})$ renormalon Beneke:1994sw ; Bigi:1994em . These
divergent behaviors cancel out in the combination as the total energy
Pineda:id ; Hoang:1998nz ; Beneke:1998rk . Besides heavy quarkonium
observables, the ${\cal O}(\Lambda_{\rm 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 Bigi:1994em ; Neubert:1994wq ;
Ball:1995wa . 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 Penin:2014zaa ; Kiyo:2015ufa ;
Beneke:2016oox ; Bazavov:2018omf ; Peset:2018ria , some of the Cabibbo-
Kobayashi-Maskawa matrix elements Hoang:1998hm ; Hoang:1998ng ;
Alberti:2014yda , and the strong coupling constant $\alpha_{s}$ Bazavov:2012ka
.

In order to improve accuracy of perturbative calculations further, it is
necessary to eliminate renormalons beyond ${\cal O}(\Lambda_{\rm QCD})$, which
requires a theoretical framework that systematically incorporates
nonperturbative QCD effects. The operator product expansion Wilson:1969zs
(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 $\Lambda_{\rm QCD}$. Thus, the OPE provides a systematic double
expansion in $\Lambda_{\rm QCD}/Q$ and $\alpha_{s}(Q^{2})$. Since, however,
renormalons in the perturbative series of the Wilson coefficients induce order
$(\Lambda_{\rm 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 Mueller:1984vh . 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 ${\cal O}(\Lambda_{\rm QCD})$ renormalon of the pole mass. In the
OPE of the QCD potential, the ${\cal O}(\Lambda_{\rm QCD})$ and ${\cal
O}(\Lambda_{\rm QCD}^{3})$ renormalons in the leading-order (LO) Wilson
coefficient were subtracted and this OPE was used to determine
$\alpha_{s}(M_{Z}^{2})$ precisely Takaura:2018lpw ; Takaura:2018vcy ;
Ayala:2020odx . Also the ${\cal O}(\Lambda_{\rm QCD}^{4})$ renormalon
contained in the lattice plaquette action was subtracted and absorbed into the
local gluon condensate Ayala:2020pxq . Many different methods to achieve
subtraction of renormalons from the LO Wilson coefficient have been developed
recently Lee:2002sn ; Ayala:2019uaw ; Ayala:2019hkn ; Ayala:2020odx ;
Ayala:2020pxq ; Takaura:2020byt ; Hayashi:2020ylq ; Hayashi:2021vdq ;
Hayashi:2021ahf ; Benitez-Rathgeb:2022yqb ; Benitez-Rathgeb:2022hfj . Among
them, the FTRS method Hayashi:2020ylq ; Hayashi:2021vdq ; Hayashi:2021ahf
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 Lee:2002sn ; Ayala:2019uaw ; Ayala:2019hkn ;
Ayala:2020odx ; Ayala:2020pxq ; Takaura:2020byt . 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 $\Lambda_{\rm 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 $\Lambda_{\rm 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-by-regions
(EBR) method Beneke:1997zp ; Smirnov:2002pj ; Jantzen:2011nz  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 $\sigma$ 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 $\sigma$ 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 $\sigma$ model are given in App. D.

##  2 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 $\widetilde{V}(q)$ by Fourier
transform as

$\displaystyle V(r)$
$\displaystyle=\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\,e^{i\vec{q}\cdot\vec{r}}\,\,\widetilde{V}(q)\,.$
(1)

We focus on the short-distance region $1/r\gg\Lambda_{\rm QCD}$, while we
neglect effects from the ultra-soft scale for simplicity. We divide Eq. (1) in
an EBR manner by employing dimensional regularization with $d=3-2\epsilon$.

$\displaystyle V(r)$ $\displaystyle=\lim_{d\to
3}\Biggl{\\{}\underbrace{\int\\!\\!\\!\frac{d^{d}\vec{q}}{(2\pi)^{d}}\,e^{i\vec{q}\cdot\vec{r}}\,\Bigl{[}\widetilde{V}(q)\Bigr{]}_{\rm
OPE}}_{\begin{array}[]{c}{\displaystyle V_{\rm hard}(r)}\\\
\text{\rotatebox{90.0}{$\in$}}\end{array}}+\underbrace{\int\\!\\!\\!\frac{d^{d}\vec{q}}{(2\pi)^{d}}\,\Bigl{[}1+i\vec{q}\\!\cdot\\!\vec{r}+\frac{1}{2}(i\vec{q}\\!\cdot\\!\vec{r})^{2}+\dots\Bigr{]}\widetilde{V}(q)}_{\begin{array}[]{c}{\displaystyle
V_{\rm soft}(r)}\\\ \text{\rotatebox{90.0}{$\in$}}\end{array}}\Biggr{\\}}\,.$
(6) $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\pm
iN_{n}(\Lambda_{\rm\overline{MS}}r)^{n}/r~{}~{}~{}~{}~{}~{}~{}~{}\stackrel{{\scriptstyle\rm
cancel}}{{\Longleftrightarrow}}~{}~{}~{}~{}~{}~{}~{}~{}\mp
iN_{n}(\Lambda_{\rm\overline{MS}}r)^{n}/r$

On the right-hand side (RHS), we divide $V(r)$ into the contribution from the
hard region, $q=|\vec{q}|\hbox{ \raise 3.0pt\hbox
to0.0pt{$>$}\raise-3.0pt\hbox{$\sim$} }1/r$, and that from the soft region,
$q\hbox{ \raise 3.0pt\hbox to0.0pt{$<$}\raise-3.0pt\hbox{$\sim$} }\Lambda_{\rm
QCD}\ll 1/r$. These contributions are denoted as $V_{\rm hard}(r)$ and $V_{\rm
soft}(r)$, respectively. In the hard region, since $q\gg\Lambda_{\rm QCD}$, we
replace $\widetilde{V}(q)$ by its high-energy expansion, namely by the OPE,
$[\widetilde{V}(q)]_{\rm OPE}$. In the soft region, instead we expand the
factor $e^{i\vec{q}\cdot\vec{r}}$ in $q$ $(\ll 1/r)$, while $\widetilde{V}(q)$
should be evaluated nonperturbatively. We obtain a series expansion in $r$,
$V_{\rm soft}(r)=\sum_{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_{\rm hard}(r)$ and $V_{\rm 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_{\rm hard}(r)$ and $V_{\rm
soft}(r)$ contain IR and UV divergences as $1/\epsilon$ poles, while in the
sum $V_{\rm hard}(r)+V_{\rm soft}(r)$ the divergences cancel since $V(r)$ is
finite as $\epsilon\to 0$. In this way, the EBR gives the OPE of $V(r)$ in
$r$.

Similarly to the appearance of $1/\epsilon$ poles, there also appear imaginary
part in the form $\pm iN_{n}(\Lambda_{\rm\overline{MS}}r)^{n}/r$ in each of
$V_{\rm hard}(r)$ and $V_{\rm soft}(r)$ for $n=1,3,5,\dots$. The imaginary
part in the Wilson coefficients of the OPE originate from renormalons. (This
feature is based on the large-$\beta_{0}$ approximation, and for simplicity we
neglect corrections to this approximation in this section.111 Namely, we
neglect corrections to the integer power behavior in $\Lambda_{\rm QCD}$ of
the renormalons. ) The imaginary part should cancel in the sum $V_{\rm
hard}(r)+V_{\rm soft}(r)$ since $V(r)$ is real.222 The $r$-independent
imaginary part (by the ${\cal O}(\Lambda_{\rm QCD})$ renormalon) in the
leading Wilson coefficient of $V_{\rm hard}(r)$ also cancels against the
imaginary part in $V_{\rm soft}(r)$ in the nonperturbative definition of
$V(r)$. This may look different from the cancellation of the ${\cal
O}(\Lambda_{\rm QCD})$ renormalon in the perturbative expansion of $V(r)$
against the corresponding renormalon in twice of the pole mass (see Sec. 1).
Note, however, that the nonperturbative definition of $V(r)$ via Wilson loop
includes the self-energies of the static charges, which are evaluated to be
zero to all orders of perturbative expansion but nonzero nonperturbatively.

In short, the imaginary part of $V_{\rm hard}(r)$ appears as follows. IR
renormalons in the leading order (LO) Wilson coefficient of
$[\widetilde{V}(q)]_{\rm OPE}$ are suppressed, as shown in ref. Sumino:2020mxk
. Therefore, the Wilson coefficient in $q$ space can be calculated accurately
in series expansion in $\alpha_{s}(q^{2})$. Then the IR renormalons in the
regularized LO Wilson coefficient of $V_{\rm hard}(r)$ can be separated as
imaginary part, using the contour integral surrounding the Landau singularity
of $\alpha_{s}(q^{2})$ in the complex $q$-plane (FTRS method) Hayashi:2020ylq
; Hayashi:2021vdq ; Hayashi:2021ahf . 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 $\widetilde{V}(q)$ can be obtained using the potential
nonrelativistic QCD EFT Brambilla:1999qa ; Brambilla:2004jw  by matching to
full QCD. Hence, we expect that we can apply the same technique also to the
subleading Wilson coefficients in $[\widetilde{V}(q)]_{\rm 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_{\rm 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
$\sigma$ 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 $\widetilde{V}(q)$. We use the terminology “hard” and
“soft” with respect to this most external momentum in Eq. (6). 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.

##  3 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.

###  3.1 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$
($\gg\Lambda_{\rm QCD}$). Its OPE is expressed as333 For many observables of
our interest, the LO nonperturbative matrix element is one, $\langle{{\cal
O}_{0}}\rangle=1$, or equivalently, the LO Wilson coefficient is identical to
perturbative evaluation of the observable. We consider those cases in this
paper.

$\Bigl{[}S(Q^{2})\Bigr{]}_{\rm
OPE}=C_{0}(Q^{2})+C_{1}(Q^{2})\frac{\langle{{\cal
O}_{1}}\rangle}{Q^{d_{1}}}+\cdots,$ (7)

where $C_{0},C_{1}$ are the Wilson coefficients and $\langle{{\cal
O}_{1}}\rangle$ is the matrix element of an operator $\mathcal{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

$\displaystyle\big{[}C_{0}(Q^{2})\big{]}_{\rm
PT}=\sum_{n=0}^{\infty}\alpha_{s}(\mu^{2})^{n+1}c_{n}(L_{Q})=\sum_{n=0}^{\infty}\alpha_{s}(Q^{2})^{n+1}c_{n}(0)\,,$
(8)

where $\alpha_{s}(\mu^{2})$ is the coupling constant renormalized at the
renormalization scale $\mu$; the coefficients $c_{n}(L_{Q})$ are polynomials
in $L_{Q}=\log\big{(}\mu^{2}/Q^{2}\big{)}$. In the second equality we set
$\mu=Q$ using the RG invariance.

In asymptotic-free theories such as QCD, naive perturbative series has IR
renormalons, namely, Eq. (8) is a divergent series where the coefficient grows
factorially, $c_{n}(0)\sim n!$. The conventional method to regularize such a
divergent series is to use the Borel resummation method. The Borel transform
of Eq. (8) is defined as

$B_{C_{0}}(u)=\sum_{n=0}^{\infty}\frac{c_{n}(L_{Q})}{n!}(u/b_{0})^{n}\,.$ (9)

Here, $b_{0}$ denotes the one-loop coefficient of the beta function. (We adopt
the convention where $b_{0}=(11-2n_{f}/3)/(4\pi)$ in the case of QCD.) Due to
the acceleration of convergence 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.444 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.  In the leading-logarithmic (LL) or large-$\beta_{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

$\big{[}C_{0}(Q^{2})\big{]}_{\pm}^{\rm
Borel}=\frac{1}{b_{0}}\int_{C_{\pm}}du\,B_{C_{0}}(u)\,e^{-u/(b_{0}\alpha_{s}(\mu^{2}))},$
(10)

where $\int_{C_{\pm}}du=\int_{0_{\pm}i0}^{\infty_{\pm}i0}du$.555 For a
convergent series without IR renormalons, Eq. (10) is simply the inverse Borel
transform.  The IR renormalon singularities are avoided by deforming the
integration contour infinitesimally above/below the positive real axis.

Since $\big{[}C_{0}(Q^{2})\big{]}_{\pm}^{\rm Borel}$ are complex conjugate to
each other, their difference $2i\delta
C_{0}\equiv\big{[}C_{0}(Q^{2})\big{]}_{+}^{\rm
Borel}-\big{[}C_{0}(Q^{2})\big{]}_{-}^{\rm Borel}$ is pure imaginary. The
contribution of each renormalon to $\delta C_{0}$ can be calculated in terms
of the residue of the renormalon or the integral around the branch cut:

$\displaystyle\delta C_{0}^{u_{*}}$
$\displaystyle=\frac{1}{2ib_{0}}\int_{C_{*}}du\,B_{C_{0}}(u)\,e^{-u/(b_{0}\alpha_{s}(\mu^{2}))}$
$\displaystyle=\frac{\pi\,N_{u_{*}}\,{u_{*}^{1+\nu_{u_{*}}}}}{b_{0}\,\Gamma(1+\nu_{u_{*}})}\left(\frac{\Lambda_{\rm\overline{MS}}^{2}}{Q^{2}}\right)^{\\!u_{*}}\\!\\!\\!(b_{0}\alpha_{s}(Q^{2}))^{\gamma_{0}/b_{0}}\bigg{(}1+\sum_{n=0}^{\infty}{s}_{n}\alpha_{s}(Q^{2})^{n+1}\bigg{)}\,.$
(11)

$C_{*}$ is shown in Fig. 1. In the second line, the form of $\delta
C_{0}^{u_{*}}$ is determined assuming that the imaginary part is canceled by
the corresponding term of OPE, Eq. (7) (in the case that there is a unique
operator with dimension $2u_{*}$) Mueller:1984vh . Apart from the overall
normalization $N_{u_{*}}$, the parameters $b_{0}$, $\nu_{u_{*}}$, $\gamma_{0}$
and the expansion coefficients $s_{n}$ can be calculated in perturbation
theory; $\Lambda_{\rm\overline{MS}}$ denotes the integration constant of the
running coupling constant in the $\overline{\rm MS}$ scheme.

Figure 1: Contour $C_{*}$ for the imaginary part of the regularized Wilson
coefficient in Eq. (11).

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:

$\big{[}C_{0}\big{]}_{\pm}^{\rm Borel}=\big{[}C_{0}\big{]}_{\rm PV}\pm i\delta
C_{0}\,,~{}~{}~{}~{}~{}[C_{0}]_{\rm PV}=\frac{1}{b_{0}}\int_{0,\rm
PV}^{\infty}du\,e^{-u/(b_{0}\alpha_{s})}B_{C_{0}}(u)\,,$ (12)

where PV stands for the principal value integration or the average of the
$C_{\pm}$ 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}]_{\rm PV}$. By definition, the renormalon contributions are minimally
subtracted from $C_{0}(Q^{2})$.

###  3.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,

$\frac{1}{2\pi
i}\int_{t_{0}-i\infty}^{t_{0}+i\infty}\\!\\!\\!\\!dt\,e^{tp^{2}}t^{u}=\frac{1}{(p^{2})^{1+u}}\frac{1}{\Gamma(-u)}~{}~{}~{};~{}~{}~{}t_{0}>0\,.$
(13)

We introduce two parameters $(a,u^{\prime})$ and define the LO Wilson
coefficient in the dual space as

$\tilde{C}_{0}(p^{2})=\frac{1}{2\pi
i}\int_{t_{0}-i\infty}^{t_{0}+i\infty}\\!\\!\\!\\!dt\,e^{tp^{2}}t^{au^{\prime}}C_{0}(t^{-a})\,,~{}~{}~{}~{}~{}~{}~{}t=Q^{-2/a}\,.$
(14)

It transforms the $\alpha_{s}$ expansion in $Q^{2}$ space Eq. (8) to the
$\alpha_{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^{\prime})$, 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.
(14) with respect to the series expansions in $\alpha_{s}(\mu^{2})$. $C_{0}$
on the RHS is replaced by $\big{[}C_{0}\big{]}_{\pm}^{\rm Borel}$, and the
left-hand side (LHS) is replaced by $\big{[}\tilde{C}_{0}\big{]}_{\pm}^{\rm
Borel}$. Then we take the difference of the $\pm$ contributions and express
the RHS as the sum of the contributions from all the singularities $u_{*}$.
Thus, by replacing $C_{0}(Q^{2})$ by $\delta
C_{0}^{u_{*}}(Q^{2})\approx\big{(}\Lambda_{\rm\overline{MS}}^{2}/Q^{2}\big{)}^{u_{*}}$
in Eq. (14), we obtain

$\displaystyle\delta\tilde{C}_{0}^{u_{*}}(p^{2})$ $\displaystyle\approx$
$\displaystyle\frac{1}{2\pi
i}\int_{t_{0}-i\infty}^{t_{0}+i\infty}\\!\\!\\!\\!\\!\\!dt\,e^{tp^{2}}t^{au^{\prime}}\big{(}\Lambda_{\rm\overline{MS}}^{2}\,t^{a}\big{)}^{u_{*}}$
(15) $\displaystyle=$
$\displaystyle\frac{1}{(p^{2})^{1+au^{\prime}}}\frac{1}{\Gamma(-a(u_{*}+u^{\prime}))}\left(\frac{\Lambda_{\rm\overline{MS}}^{2}}{p^{2a}}\right)^{u_{*}}\,.$
(16)

Eq. (16) is zero if $u_{*}+u^{\prime}=n/a$ for $n\in\mathbb{Z}_{\geq 0}$.
Namely, we can choose $(a,u^{\prime})$ such that there no longer exist
singularities in the Borel plane which induce factorial growth of
$\tilde{C}_{0}(p^{2})$ if the imaginary part of
$\big{[}C_{0}(Q^{2})\big{]}_{\pm}^{\rm Borel}$ are of the form
$(\Lambda_{\rm\overline{MS}}^{2}/Q^{2})^{n/a-u^{\prime}}$. Thus, IR
renormalons of $\tilde{C}_{0}(p^{2})$ can be suppressed. It is possible to
suppress more than one IR renormalons simultaneously with an appropriate
choice of $(a,u^{\prime})$ Hayashi:2020ylq ; Hayashi:2021vdq ; Hayashi:2021ahf
. In principle, it is also possible to take into account the higher-order
$\alpha_{s}(Q^{2})$ corrections to the simple power of
$\Lambda_{\rm\overline{MS}}^{2}/Q^{2}$ and anomalous dimension effects order
by order. (See the end of this subsection.) Hence, the perturbative expansion
of $\tilde{C}_{0}$ can be regarded as a convergent series.

While $\big{[}\tilde{C}_{0}\big{]}_{\mathrm{PT}}$ is a convergent series,
$\big{[}C_{0}\big{]}_{\mathrm{PT}}$ is not, due to the IR renormalons. The
difference emerges as follows.

$\big{[}\tilde{C}_{0}(p^{2})\big{]}_{\rm PT}=\frac{1}{2\pi
i}\int_{t_{0}-i\infty}^{t_{0}+i\infty}\\!\\!\\!\\!dt\,e^{tp^{2}}t^{au^{\prime}}\big{(}{\mu^{2}}t^{a}\big{)}^{\hat{H}}\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\mu^{2})^{n+1}\,,$
(17)

where, using the RG invariance, we have rewritten Eq. (8) as

$\displaystyle\big{[}C_{0}(Q^{2})\big{]}_{\rm
PT}=\bigg{(}\frac{\mu^{2}}{Q^{2}}\bigg{)}^{\hat{H}}\sum_{n=0}^{\infty}\alpha_{s}(\mu^{2})^{n+1}c_{n}(0)\,,$
(18)

with

$\hat{H}=-\beta(\alpha_{s})\,\frac{\partial}{\partial\alpha_{s}(\mu^{2})}\,,~{}~{}~{}~{}~{}~{}~{}~{}\beta(\alpha_{s})=\mu^{2}\frac{d}{d\mu^{2}}\alpha_{s}(\mu^{2})=-\sum_{i=0}^{\infty}b_{i}\alpha_{s}^{i+2}(\mu^{2})\,.$
(19)

Integrating over $t$ we have

$\displaystyle\big{[}\tilde{C}_{0}(p^{2})\big{]}_{\rm PT}$
$\displaystyle=\frac{1}{(p^{2})^{1+au^{\prime}}}\,e^{\hat{H}\,\log(\mu^{2}/p^{2a})}\,\frac{1}{\Gamma(-a(\hat{H}+u^{\prime}))}\,\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\mu^{2})^{n+1}$
$\displaystyle=\frac{1}{(p^{2})^{1+au^{\prime}}}\,e^{\hat{H}\,\log(\mu^{2}/p^{2a})}\,\sum_{n=0}^{\infty}\tilde{c}_{n}(0)\alpha_{s}(\mu^{2})^{n+1}\,,$
(20)

where $\tilde{c}_{n}(0)$’s are computed by comparing $\alpha_{s}(\mu^{2})$
expansion (using the expansion in $\hat{H}$) of the first and second lines. We
can see that the coefficients of the $\alpha_{s}(\mu^{2})$ expansion of Eq.
(20) have $\log^{n}(\mu^{2}/p^{2a})$ dependence. The Laplace integral of such
logarithms in the region $p^{2a}\hbox{ \raise 3.0pt\hbox
to0.0pt{$<$}\raise-3.0pt\hbox{$\sim$} }\Lambda_{\rm\overline{MS}}^{2}$ gives a
factorial behavior $\sim n!$, which is the origin of the IR renormalons in
$C_{0}$. Noting that $\tilde{C}_{0}$ is RG invariant, we can resum the
logarithms by setting $\mu^{2}=p^{2a}$ in Eq. (20):

$\big{[}\tilde{C}_{0}(p^{2})\big{]}_{\rm
PT}=\frac{1}{(p^{2})^{1+au^{\prime}}}\,\sum_{n=0}^{\infty}\tilde{c}_{n}(0)\alpha_{s}(p^{2a})^{n+1}\,.$
(21)

As a result, the higher power logarithms are turned into the behavior of the
running coupling constant $\alpha_{s}(p^{2a})$ around the Landau singularity,
which becomes an alternative source of renormalons upon considering the
Laplace integral. The expansion of $\tilde{C}_{0}(p^{2})$ in
$\alpha_{s}(p^{2a})$ is also a convergent series for a sufficiently small
$\alpha_{s}(p^{2a})$. In practical calculations, the original series Eq. (18)
is known only up to a certain order $n\leq k$. Then we compute
$\tilde{c}_{n}(0)$ also up to $n=k$. Thus, using the running coupling constant
of the NkLL approximation and truncating $\tilde{C}_{0}(p^{2})$ at order
$\alpha_{s}(p^{2a})^{k+1}$, one can compute a good approximation (at the level
of NkLL) of the true $\tilde{C}_{0}(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 $\tilde{C}_{0}(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_{\pm}$ as Eq. (10), the regularized Wilson
coefficient $\big{[}C_{0}(Q^{2})\big{]}_{\pm}$ is given by

$\displaystyle\big{[}C_{0}(Q^{2})\big{]}_{\pm}$
$\displaystyle=t^{-au^{\prime}}\int_{C_{\mp}}dp^{2}\,e^{-tp^{2}}\tilde{C}_{0}(p^{2})$
(22) $\displaystyle=\big{[}C_{0}(Q^{2})\big{]}_{\rm PV}\pm\,i\,\delta
C_{0}(Q^{2}),$ (23) $\big{[}C_{0}(Q^{2})\big{]}_{\rm
PV}=t^{-au^{\prime}}\int_{0,\rm
PV}^{\infty}\\!\\!\\!\\!dp^{2}\,e^{-tp^{2}}\tilde{C}_{0}(p^{2}),\quad\delta
C_{0}(Q^{2})=\frac{t^{-au^{\prime}}}{2i}\int_{C_{-}-C_{+}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!dp^{2}\,e^{-tp^{2}}\tilde{C}_{0}(p^{2})\,,$
(24)

where PV stands for the real part of the integration. It can be proven that
$\big{[}C_{0}(Q^{2})\big{]}_{\pm}$ given here agrees with that defined via the
Borel resummation in the previous subsection. (See below.) Keeping the
integral value, the integral path $C_{\mp}$ can be deformed away from the
Landau singularity in the right-half plane, such that
$\bigl{|}\alpha_{s}(p^{2a})\bigr{|}$ is sufficiently small along the entire
integral path.666 We neglect effects of the singularities which possibly exist
in the right-half plane (${\rm Re}\,p^{2}>0$) other than the Landau
singularity. This would be reasonable for observables defined in Euclidean
momentum region (e.g., Adler function, QCD potential). At the level of NkLL
approximation in QCD ($k\leq 4$), singularities which affect our present
argument do not exist; see App. C of ref. Hayashi:2021vdq  for details.
Hence, we obtain a good estimate of $\big{[}C_{0}(Q^{2})\big{]}_{\rm PV}$ and
$\delta C_{0}(Q^{2})$ using the series expansion of $\tilde{C}_{0}(p^{2})$. In
particular, we call this method of computing $\big{[}C_{0}(Q^{2})\big{]}_{\rm
PV}$ dual space renormalon subtraction (DSRS) method.

As compared to the FTRS method Hayashi:2020ylq ; Hayashi:2021vdq ;
Hayashi:2021ahf , 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
$\big{[}C_{0}(Q^{2})\big{]}_{\pm}$ 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. Hayashi:2021vdq . Similarly, we can
incorporate the corrections to the integer power forms proportional to
$(\Lambda_{\rm QCD}/Q)^{n}$ of the imaginary part, given by the anomalous
dimension and expansion in $\alpha_{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. Hayashi:2021vdq . (See, however, the discussion at the end of
Sec. 5.)

###  3.3 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.

$\displaystyle
S(Q^{2})=t^{-au^{\prime}}\int_{0}^{\infty}dp^{2}\,e^{-tp^{2}}\tilde{S}(p^{2})\,,~{}~{}~{}~{}~{}~{}~{}t=Q^{-2/a}\,,$
(25) $\displaystyle\Longleftrightarrow~{}~{}~{}$
$\displaystyle\tilde{S}(p^{2})=\frac{1}{2\pi
i}\int_{t_{0}-i\infty}^{t_{0}+i\infty}\\!\\!\\!\\!dt\,e^{tp^{2}}t^{au^{\prime}}S(t^{-a})\,,~{}~{}~{}~{}~{}~{}~{}t_{0}>0\,.$
(26)

Below we assume that the above integrals are well defined.

In the case $Q\gg\Lambda_{\rm QCD}$, we obtain the OPE of $S(Q^{2})$ by
evaluating the $p^{2}$ integral in Eq. (25) by expansion in $1/Q$ using a
generalized EBR method. We divide the integral into the contribution from the
hard region, $p^{2}\hbox{ \raise 3.0pt\hbox
to0.0pt{$>$}\raise-3.0pt\hbox{$\sim$} }t^{-1}=Q^{2/a}$, and that from the soft
region, $p^{2}\ll t^{-1}$, after regularizing the integral.

$\displaystyle\Bigl{[}S(Q^{2})\Bigr{]}_{\rm OPE}=\lim_{\delta\to
0}\biggl{\\{}$ $\displaystyle
t^{-au^{\prime}}\int_{C_{\mp}}dp^{2}\,p^{-2\delta}\,e^{-tp^{2}}\left[\tilde{S}(p^{2})\right]_{\rm
OPE}$
$\displaystyle+t^{-au^{\prime}}\int_{C_{\mp}}dp^{2}\,p^{-2\delta}\,\left[1-tp^{2}+\frac{1}{2}(tp^{2})^{2}-\dots\right]\tilde{S}(p^{2})\biggr{\\}}\,.$
(27)

In the hard region, $\tilde{S}(p^{2})$ is replaced by its OPE, namely
expansion in $1/p^{2}$ and $\alpha_{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 $\delta$
analogously to the dimensional regularization, while the Landau singularity
included in $\bigl{[}\tilde{S}(p^{2})\bigr{]}_{\rm OPE}$ is regularized by
deforming the integral contour to $C_{\mp}$. (Or, we can deform the contour
further away from the Landau singularity without changing the integral
values.) The same contour $C_{\mp}$ is chosen in both hard and soft
contributions.

Each term of $\bigl{[}\tilde{S}(p^{2})\bigr{]}_{\rm
OPE}=\sum_{i}\tilde{C}_{i}(p^{2})\langle O_{i}\rangle/p^{n_{i}}$ is uniquely
determined by the dual transform of each term of
$\bigl{[}S(Q^{2})\bigr{]}_{\rm OPE}$ given by Eq. (7).777 We assume that the
OPE is equivalent to the high-energy expansion of $S(Q^{2})$ or
$\tilde{S}(p^{2})$.  Hence, the Wilson coefficients $\tilde{C}_{i}(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 $\bigl{[}S(Q^{2})\bigr{]}_{\rm
OPE}$, we evaluate $\tilde{C}_{i}(p^{2})$ in expansion in
$\alpha_{s}(p^{2a})$, which generates the Landau singularity in the integrand.
Convergence of the first term of Eq. (27) 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 $\alpha_{s}(p^{2a})$.

Since the original integral in Eq. (25) has no singularities and is well
defined, the IR and UV divergences (regularized by $\delta$) 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 $\sigma$ model and QCD
(in NkLL approximation) in the following sections. We will see that the sign
of the imaginary part of $\delta$ must be correlated with the choice of the
regularized contour $C_{\mp}$ in order to guarantee the validity of the EBR
method.

The evaluation of the soft contribution, the second term of Eq. (27), involves
nonperturbative evaluation of $\tilde{S}(p^{2})$, which is possible only in
limited cases. We will confirm desired properties for some examples in the
$O(N)$ nonlinear $\sigma$ 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.

##  4 OPEs in nonlinear $\sigma$ model

In this section we demonstrste usefulness of the dual space approach through a
study of the OPEs in the 2D $O(N)$ nonlinear $\sigma$ 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 $\sigma$ as another example.

###  4.1 The model

We consider the $O(N)$ nonlinear $\sigma$ model in two dimensions in the
large-$N$ limit Bardeen:1976zh . The partition function of the model is given
by

$Z=\int{\cal D}\sigma{\cal D}\alpha\exp\bigg{[}-\frac{1}{2g_{0}}\int
d^{2}x\bigg{\\{}\partial_{\mu}\sigma^{a}\partial_{\mu}\sigma^{a}+{\alpha}\bigg{(}\sigma^{a}\sigma^{a}-N\bigg{)}\bigg{\\}}\bigg{]},$
(28)

where $g_{0}$ is the bare coupling constant. The Lagrange multiplier
$\alpha(x)$ is introduced to take into account the constraint
$\sum_{a=1}^{N}\sigma^{a}\sigma^{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 $\alpha$ as

$\langle\alpha\rangle=m^{2}\equiv\mu^{2}\,\exp\big{[}-4\pi/g(\mu^{2})\big{]}\,,$
(29)

where $g(\mu^{2})$ denotes the coupling constant in the $\overline{\rm MS}$
scheme renormalized at the renormalization scale $\mu$. For convenience, we
use the rescaled coupling constant below:

$\bar{g}(\mu^{2})=\frac{g(\mu^{2})}{4\pi}=\frac{1}{\log(\mu^{2}/m^{2})}\,.$
(30)

It diverges at the Landau pole $\mu^{2}=m^{2}$. The propagater of
$\delta\alpha\equiv\alpha-\langle\alpha\rangle$ is given by

$D_{\alpha}(p^{2})=\int
d^{2}x\,e^{-ip\cdot(x-y)}\,\langle\delta\alpha(x)\delta\alpha(y)\rangle=\frac{1}{N}\,\frac{4\pi\sqrt{p^{2}(p^{2}+4m^{2})}}{\log\bigg{[}\frac{\sqrt{p^{2}+4m^{2}}+\sqrt{p^{2}}}{\sqrt{p^{2}+4m^{2}}-\sqrt{p^{2}}}\bigg{]}}\,.$
(31)

We investigate the OPEs of observables in this model using the dual space
approach.

###  4.2 Local condensate of $\delta\alpha$ regularized by gradient flow

The first example is a vacuum expectation value of a composite operator of
$\delta\alpha(x)$ regularized by gradient flow,

$\displaystyle I(t)$ $\displaystyle\equiv
Nt^{2}\cdot\langle\delta\alpha(x)\,e^{t\,\Box_{x}}\delta\alpha(x)\rangle=Nt^{2}\int_{0}^{\infty}{dp^{2}}\,\frac{D_{\alpha}(p^{2})}{4\pi}\,e^{-tp^{2}}\,,$
(32)

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.888 The OPE
of a conceptually similar observable, $\langle\delta\alpha(x)^{2}\rangle$ with
a UV cutoff regularization, has been studied in refs. Bardeen:1976zh ;
David:1982qv ; David:1985xj . 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 $\sigma$ model has been considered in ref. Ishikawa:2019tnw  in the
gradient-flow regularization, which discusses the renormalon structure of the
$CP^{N}$ theory on $\mathbb{R}\times S^{1}$.  In this example, we give its
dual transform by

$\displaystyle\tilde{I}(p^{2})=\frac{1}{2\pi
i}\int_{t_{0}-i\infty}^{t_{0}+i\infty}\\!\\!\\!\\!dt\,e^{tp^{2}}\frac{I(t)}{t^{2}}=N\frac{D_{\alpha}(p^{2})}{4\pi}\,,$
(33) $\displaystyle
I(t)=t^{2}\int_{0}^{\infty}{dp^{2}}\,\tilde{I}(p^{2})\,e^{-tp^{2}}\,,$ (34)

with the choice $(a,u^{\prime})=(1,-2)$. Recall Eq. (26).

The OPE of $\tilde{I}(p^{2})$ (the expansion in $m^{2}/p^{2}$) is given by

$\displaystyle\tilde{I}(p^{2})=p^{2}\Biggl{[}\bar{g}(p^{2})+\frac{m^{2}}{p^{2}}\left(2\bar{g}(p^{2})-2\bar{g}(p^{2})^{2}\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(\frac{m^{2}}{p^{2}}\right)^{2}\left(-2\bar{g}(p^{2})-\bar{g}(p^{2})^{2}+4\bar{g}(p^{2})^{3}\right)+\cdots\Biggr{]}\,.$
(35)

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 $\bar{g}(p^{2})$ Novikov:1984ac .

In comparison, the OPE of $I(t)$ has renormalons. The perturbative expansion
in $\bar{g}(\mu^{2})$ of the LO Wilson coefficient of $I(t)$ is given by the
Laplace transform of that of $\tilde{I}(p^{2})$, and its Borel transform is
expressed by999 Here, we define the Borel transform in the nonlinear $\sigma$
model by replacing $\alpha_{s}(\mu^{2})$ with $\bar{g}(\mu^{2})$ and $b_{0}$
with $1$ in Eqs. (8) and (9).

$(\mu^{2}t)^{u}\,\Gamma(2-u).$ (36)

This shows that the IR renomalons are located at $u=2,3,4,\cdots$ and the
leading renormalon uncertainty is $\mathcal{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(\mu^{2}/p^{2})^{n}$ encoded in
$\bar{g}(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. (33) (although the above procedure may look reversed). The
parameters $(a,u^{\prime})=(1,-2)$ in Eq. (33) are chosen to suppress the
renormalons at $u=2,3,4,\cdots$. 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_{\alpha}(p^{2})$ 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. (34):

$\displaystyle\big{[}I(t)\big{]}_{\rm OPE}=\lim_{\delta\to
0}\left\\{I^{\pm}_{h}(t)+I^{\pm}_{s}(t)\right\\}\,,$ (37) $\displaystyle
I^{\pm}_{h}(t)=t^{2-\delta}\int_{C_{\mp}}dp^{2}\,(p^{2})^{-\delta_{\pm}}\left[\tilde{I}(p^{2})\right]_{\rm
OPE}\,e^{-tp^{2}}\,,$ (38) $\displaystyle
I^{\pm}_{s}(t)=t^{2-\delta}\int_{C_{\mp}}dp^{2}\,(p^{2})^{-\delta_{\pm}}\tilde{I}(p^{2})\,\left\\{1-tp^{2}+\frac{1}{2}(tp^{2})^{2}-\dots\right\\}\,,$
(39)

where

$\delta_{\pm}=\delta\pm
i0\,,~{}~{}~{}~{}~{}\delta\in\mathbb{R}\backslash\\{0\\}\,.$ (40)

The EBR realizes the cancellation of IR renormalons. IR and UV divergences
also appear, which are regularized by $\delta$ and canceled in
$\big{[}I(t)\big{]}_{\rm OPE}$.

For a while we focus on the first few terms of the OPE of $I(t)$:

$\big{[}I(t)\big{]}_{\rm
OPE}=C_{0}^{(h)}(t)+C_{1}^{(h)}(t)\langle\alpha\rangle
t+C_{2}^{(h)}(t)\langle\alpha\rangle^{2}t^{2}+C_{2}^{(s)}\langle\delta\alpha^{2}\rangle
t^{2}+{\cal O}(t^{3}).$ (41)

Explicitly they are given by

$\displaystyle C_{0}^{(h)}(t)$
$\displaystyle=t^{2}\int_{C_{\mp}}dp^{2}p^{2}\bar{g}(p^{2})e^{-tp^{2}},$ (42)
$\displaystyle C_{1}^{(h)}(t)$
$\displaystyle=t\int_{C_{\mp}}dp^{2}\,\left(2\bar{g}(p^{2})-2\bar{g}(p^{2})^{2}\right)e^{-tp^{2}},$
(43) $\displaystyle C_{2}^{(h)}(t)$
$\displaystyle=t^{-\delta}\int_{C_{\mp}}\frac{dp^{2}}{(p^{2})^{1+\delta_{\pm}}}\,\left(-2\bar{g}(p^{2})-\bar{g}(p^{2})^{2}+4\bar{g}(p^{2})^{3}\right)e^{-tp^{2}},$
(44)

and

$\langle\delta\alpha^{2}\rangle=\int_{C_{\mp}}dp^{2}(p^{2})^{-\delta_{\pm}}\tilde{I}(p^{2})\,,$
(45)

with $C_{2}^{(s)}=1$. The Wilson coefficients are denoted as $C_{i}^{(h,s)}$
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
$C_{2}^{(s)}$. There appear two nonperturbative matrix elements,
$\langle\alpha(x)\rangle=m^{2}$ and $\langle\delta\alpha(x)^{2}\rangle={\cal
O}(m^{4})$, up to ${\cal O}(t^{2})$. We omit $\delta_{\pm}$ when not necessary
(in the absence of IR/UV divergences).

Each term of the OPE (41) has ambiguities, which stem from the regions where
$p^{2}$-integrals have divergences. The regions of $p^{2}$ which cause
ambiguities are as follows:

$\displaystyle C_{0}^{(h)}:$ $\displaystyle~{}p^{2}\sim m^{2},$ $\displaystyle
C_{1}^{(h)}:$ $\displaystyle~{}p^{2}\sim m^{2},$ $\displaystyle C_{2}^{(h)}:$
$\displaystyle~{}p^{2}\sim m^{2},\,p^{2}\sim 0,$
$\displaystyle\langle\delta\alpha^{2}\rangle:$
$\displaystyle~{}p^{2}\sim\infty.$ (46)

The hard quantities have ambiguities associated with the Landau pole
singularity ($p^{2}\sim m^{2}$) and IR divergences ($p^{2}\sim 0$), whereas
the soft quantities (or condensate) have UV divergences ($p^{2}\sim\infty$).
In other words, the ambiguities originate from the region where their original
expansions are not valid. By focusing on the contributions up to
$\mathcal{O}(t^{2}m^{4})$ in the OPE, these imaginary part and divergences are
given by

$\displaystyle C_{0}^{(h)}|_{\rm Landau}=t^{2}\int_{p_{1}^{2}\mp
i0}^{p_{2}^{2}\mp i0}dp^{2}p^{2}\bar{g}(p^{2})=\pm i\pi t^{2}m^{4}+\cdots\,,$
(47) $\displaystyle C_{2}^{(h)}|_{\rm Landau}=\int_{p_{1}^{2}\mp
i0}^{p_{2}^{2}\mp i0}\frac{dp^{2}}{p^{2}}\,\left(-2\bar{g}(p^{2})\right)=\mp
2i\pi+\cdots\,,$ (48) $\displaystyle C_{2}^{(h)}|_{\rm
soft}=t^{-\delta}\int_{0}^{p_{1}^{2}}\frac{dp^{2}}{(p^{2})^{1+\delta_{\pm}}}\,\left(-2\bar{g}(p^{2})\right)=2\int_{0}^{\infty}d\alpha\left(\frac{p_{1}^{2}}{m^{2}}\right)^{\alpha}\frac{1}{\alpha-\delta_{\pm}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}=-2\log{(-\delta_{\pm}})+\cdots\,,$
(49) $\displaystyle\langle\delta\alpha^{2}\rangle|_{\rm hard}$
$\displaystyle=\int_{p_{2}^{2}}^{\infty}dp^{2}(p^{2})^{-\delta_{\pm}}\left[p^{2}\bar{g}(p^{2})+\frac{m^{4}}{p^{2}}(-2\bar{g}(p^{2}))\right]$
$\displaystyle=\mp i\pi m^{4}+2m^{4}\log{\delta_{\pm}}+\cdots\,,$ (50)

where $0<p_{1}^{2}<m^{2}<p_{2}^{2}<1/t$. In the above equations, the dots
“$\cdots$” represent real finite contributions. In evaluating the integrals
for $C_{0,2}^{(h)}$ around the Landau pole we expand the exponential factor
$e^{-tp^{2}}\simeq 1$. The ambiguities found in $C_{0}^{(h)}|_{\rm Landau}$
are the same as the renormalon uncertainties derived from Eq. (36). It turns
out that $C_{1}^{(h)}|_{\rm Landau}$ has higher order imaginary part which we
neglect here. Also the higher order terms in $\bar{g}(p^{2})$ in the integrand
for $C_{2}^{(h)}|_{\rm Landau}$ give higher order imaginary part. We note
that, in evaluating the integral for $\langle\delta\alpha^{2}\rangle$ in the
region $p^{2}\sim\infty$, it is justified to expand $\tilde{I}(p^{2})$ in
$m^{2}/p^{2}$ as in Eq. (35). We therefore have the same integrands in the
condensate (50) as those in the Wilson coefficients (47) and (49). In the
integrand of Eq. (50), we only show the terms which give imaginary part or UV
divergences. The imaginary part in Eq. (50) arises as Beneke:1998ui

$\displaystyle\int_{p_{2}^{2}}^{\infty}dp^{2}(p^{2})^{1-\delta_{\pm}}\bar{g}(p^{2})$
$\displaystyle=\int_{0}^{\infty}d\alpha\,m^{2\alpha}\int_{p_{2}^{2}}^{\infty}dp^{2}(p^{2})^{1-\delta_{\pm}-\alpha}$
$\displaystyle=\int_{0}^{\infty}d\alpha\,(p_{2}^{2})^{2-\alpha-\delta_{\pm}}m^{2\alpha}\frac{1}{\alpha+\delta_{\pm}-2}=\mp
i\pi m^{4}+\cdots\,,$ (51)

and the appearance of the $\log{\delta_{\pm}}$ can be understood by a similar
calculation to Eq. (49). One can see that the above imaginary part and
$\log{\delta}$ terms found in Eqs. (47)–(50) cancel in the OPE (41), as
expected.

The above cancellation indicates the validity of the EBR method. In fact the
validity of the EBR relation Eqs. (37)–(39) can be proved following the
strategy of ref. Jantzen:2011nz . (See App. C.) The crucial identity which
guarantees the validity is given by101010 This identity can be shown by
analytic continuation of $a_{\pm}$: Change the integral variable to
$z=\log(p^{2}/m^{2})$, and close the contour in the upper/lower half $z$-plane
for $a_{R}=-1$, $a_{I}>0$. Note that it is crucial to take the sign of
$C_{\pm}$ and that of $a_{\pm}$ to be the same.

$\int_{C_{\pm}}dp^{2}\,\frac{(p^{2})^{a_{\pm}}}{\\{\log(p^{2}/m^{2})\\}^{n}}=0\,,$
(52)

where

$n\in\mathbb{Z}_{>0}\,,~{}~{}~{}~{}a_{\pm}\equiv a_{R}\pm
ia_{I}~{}~{}\text{with}~{}~{}a_{R},a_{I}\in\mathbb{R},~{}~{}a_{I}>0\,.$ (53)

The cancellation explicitly confirmed above is a consequence of the particular
identities

$\int_{C_{\mp}}dp^{2}(p^{2})^{-\delta_{\pm}}p^{2}\bar{g}(p^{2})=0,$ (54)
$\int_{C_{\mp}}dp^{2}(p^{2})^{-1-\delta_{\pm}}(-2\bar{g}(p^{2}))=0,$ (55)

which lead to $\pm i\pi m^{4}\mp i\pi m^{4}=0$ and $\mp
2i\pi-2\log{(-\delta_{\pm})}+2\log{\delta_{\pm}}=0$, respectively. Eq. (52)
also explains the simultaneous absence of the imaginary part associated with
the Landau pole $C_{1}^{(h)}|_{\rm Landau}=\int_{p_{1}^{2}\mp
i0}^{p_{2}^{2}\mp
i0}dp^{2}(p^{2})^{-\delta_{\pm}}\,\left(2\bar{g}(p^{2})-2\bar{g}(p^{2})^{2}\right)$
and the UV divergence in the condensate, which could be caused by
$\int_{p_{2}^{2}}^{\infty}dp^{2}(p^{2})^{-\delta_{\pm}}\,\left(2\bar{g}(p^{2})-2\bar{g}(p^{2})^{2}\right)$;
if either of them is zero, the other should also be zero.

Using the relation (52), we can discuss in a general manner the cancellation
of the imaginary part in $I_{h}^{\pm}(t)$ and $I_{s}^{\pm}(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}\to 0$, (ii) $p^{2}\sim
m^{2}$, (iii) $p^{2}\to\infty$. The imaginary part from (i) and (ii) are
included in the Wilson coefficients of $I_{h}^{\pm}$, while that from (iii) is
included in the matrix elements of $I_{s}^{\pm}$.111111 These features are
somewhat confusing, since they stem from “IR region of the hard contribution”
and “UV region of the soft contribution,” respectively.  In fact, the above
relation ensures the following identity:

$\bigg{(}\int_{0}^{p_{1}^{2}}+\int_{p_{2}^{2}}^{\infty}\bigg{)}dp^{2}\,\bigg{(}\frac{m^{2}}{p^{2}}\bigg{)}^{n}\frac{(p^{2})^{1+\ell-\delta_{\pm}}}{\log(p^{2}/m^{2}_{\pm})^{k+1}}=-\int_{p_{1}^{2}}^{p_{2}^{2}}dp^{2}\,\bigg{(}\frac{m^{2}}{p^{2}}\bigg{)}^{n}\frac{(p^{2})^{1+\ell-\delta_{\pm}}}{\log(p^{2}/m^{2}_{\pm})^{k+1}},$
(56)

where $n,k,\ell\in\mathbb{Z}_{\geq 0}$. $n$ and $k$ correspond to the
expansion orders of $\tilde{I}(p^{2})$ in $m^{2}/p^{2}$ and $\bar{g}(p^{2})$,
respectively, i.e., the ${\cal O}(\bar{g}(p^{2})^{k+1}/p^{2(n-1)})$ term of
$\tilde{I}(p^{2})$, whereas $\ell$ 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 $\tilde{I}(p^{2})$ in $m^{2}/p^{2}$ and
$\bar{g}(p^{2})$, and $e^{-tp^{2}}$ in $tp^{2}$; for this reason it is
sufficient to assume the form of the integrand of Eq. (56). The integrals on
both sides are the contributions to ${\cal O}(t^{\ell+2})$ term of $I(t)$. For
notational convenience, we have taken the integral path on the positive real
axis and instead shifted the Landau pole infinitesimally above or below the
axis, $m^{2}_{\pm}\equiv m^{2}\pm i0$. The imaginary part of the LHS of Eq.
(56) is calculated as

$\displaystyle{\rm
Im}\int_{0}^{\infty}\\!\\!\\!d\alpha\bigg{[}\int_{0}^{p_{1}^{2}}\\!\\!\\!dp^{2}\frac{(-\alpha)^{k}(p^{2})^{1+\ell-\delta_{\pm}}}{k!}\bigg{(}\frac{m^{2}_{\pm}}{p^{2}}\bigg{)}^{n-\alpha}\\!+\int_{p_{2}^{2}}^{\infty}\\!\\!\\!dp^{2}\frac{\alpha^{k}(p^{2})^{1+\ell-\delta_{\pm}}}{k!}\bigg{(}\frac{m^{2}_{\pm}}{p^{2}}\bigg{)}^{n+\alpha}\bigg{]}$
$\displaystyle=\pm\pi(m^{2})^{\ell+2}\frac{(\ell+2-n)^{k}}{k!}\big{(}\delta_{\ell+2-n\leq
0}+\delta_{\ell+2-n>0}\big{)}.$ (57)

The first term is from the hard contribution (in $I_{h}^{\pm}$) and the second
is from the soft contribution (in $I_{s}^{\pm}$).121212 In the last line of
Eq. (57), we chose $\delta>0$. The sum is unchanged if we choose $\delta<0$.
Note that each contribution is independent of $p_{1}^{2},p_{2}^{2}$, as it
should be. The RHS of Eq. (56) contains the imaginary part due to the Landau
pole of the running coupling constant $\bar{g}(p^{2})$, which is calculated as

$\displaystyle\mp\frac{1}{2i}\bigg{(}\int_{p_{1}^{2}+i0}^{p_{2}^{2}+i0}-\int_{p_{1}^{2}-i0}^{p_{2}^{2}-i0}\bigg{)}dp^{2}\,\bigg{(}\frac{m^{2}}{p^{2}}\bigg{)}^{n}\frac{(p^{2})^{1+\ell-\delta_{\pm}}}{\log(p^{2}/m^{2})^{k+1}}$
$\displaystyle=\pm\frac{1}{2i}\int_{C_{L}}dp^{2}\,\bigg{(}\frac{m^{2}}{p^{2}}\bigg{)}^{n}\frac{(p^{2})^{1+\ell-\delta_{\pm}}}{\log(p^{2}/m^{2})^{k+1}}=\pm\pi(m^{2})^{\ell+2}\frac{(\ell+2-n)^{k}}{k!}\,,$
(58)

where $C_{L}$ is the closed path surrounding the Landau pole counterclockwise
in the complex $p^{2}$ plane (Fig. 2).

Figure 2: Contour $C_{L}$ surrounding the Landau pole in the complex $p^{2}$
plane.

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,\ell$. Using this relation, we can show that the imaginary
part of $I_{h}^{\pm}(t)+I_{s}^{\pm}(t)$ is identically zero, consistently with
the fact that $I(t)$ is real:

$\displaystyle{\rm Im}\bigg{[}I_{h}^{\pm}(t)+I_{s}^{\pm}(t)\bigg{]}$
$\displaystyle=\pm\pi
m^{4}t^{2}\sum_{\ell=0}^{\infty}\frac{(-m^{2}t)^{\ell}}{\ell!}\sum_{n=0}^{\infty}\sum_{k=0}^{n}\frac{d_{n,k}(0)}{k!}(\ell+2-n)^{k}\big{(}\delta_{\ell+2-n\leq
0}+\delta_{\ell+2-n>0}\big{)}$ $\displaystyle~{}~{}~{}\mp\pi
m^{4}t^{2}\sum_{\ell=0}^{\infty}\frac{(-m^{2}t)^{\ell}}{\ell!}\sum_{n=0}^{\infty}\sum_{k=0}^{n}\frac{d_{n,k}(0)}{k!}(\ell+2-n)^{k}$
$\displaystyle=0\,,$ (59)

where $d_{n,k}(0)$ represent perturbative coefficients appearing in
$\tilde{I}(p^{2})$. See Eq. (66). 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,\ell$ 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)^{\ell}$ vanishes as

$\displaystyle\sum_{n=0}^{\infty}\sum_{k=0}^{n}\frac{d_{n,k}(0)}{k!}(\ell+2-n)^{k}=\sum_{n=0}^{\infty}X_{n}(\ell+2-n)$
$\displaystyle=\sum_{n=0}^{2\ell+4}(-1)^{n}\frac{2(6+2n^{2}+7\ell+2\ell^{2}-4n(2+\ell))}{(4+2\ell)(3+2\ell)}\,_{4+2\ell}C_{n}=0\,,$
(60)

for arbitrary $\ell\in\mathbb{Z}_{\geq 0}$. We used the explicit result for
$d_{n,k}$ given below [Eq. (65)]. The last equality is shown using an identity
$\sum_{n=0}^{N}{{}_{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. (58). Since before the OPE
$\tilde{I}(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 (58) 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 $\tilde{I}(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
$\bar{g}(p^{2})$) but within the double expansion in $\bar{g}(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{\delta}$ 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. (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\delta_{\pm}}$ and the $\log\delta$ 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.)

#### 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 $\bar{g}(p^{2})$. Here, we mainly study expressions
with the fixed-scale coupling constant $\bar{g}(\mu^{2})$. The OPE of
$\tilde{I}(p^{2})$ is expressed in the form

$\displaystyle\left[\tilde{I}(p^{2})\right]_{\rm
OPE}=p^{2}\sum_{n=0}^{\infty}\biggl{(}\frac{m^{2}}{p^{2}}\biggr{)}^{n}\,\tilde{C}_{n}(p^{2})\,,$
(61)

with the Wilson coefficients in the dual space given by

$\displaystyle\tilde{C}_{n}(p^{2})=\sum_{k=0}^{\infty}\bar{g}(\mu^{2})^{k+1}\,d_{n,k}(L_{p})\,,~{}~{}~{}~{}~{}L_{p}=\log\left(\frac{\mu^{2}}{p^{2}}\right)\,.$
(62)

They are obtained by expanding Eq. (31) in $m^{2}/p^{2}$ and
$1/\log(p^{2}/m^{2})=\bar{g}(p^{2})$, and then by re-expressing
$\bar{g}(p^{2})=\bar{g}(\mu^{2})/[1+\bar{g}(\mu^{2})\log(p^{2}/\mu^{2})]$ in a
geometric series of $\bar{g}(\mu^{2})$. Inserting Eq. (61) into Eq. (38) we
obtain the Wilson coefficients of the hard contribution,

$I^{\pm}_{h}=\sum_{n=0}^{\infty}(m^{2}t)^{n}\left[C_{n}^{(h)}(t)\right]_{\pm}\,,~{}~{}~{}~{}~{}\left[C_{n}^{(h)}(t)\right]_{\pm}\\!=t^{2-\delta-n}\\!\int_{C_{\mp}}\\!\\!\\!dp^{2}\,(p^{2})^{1-\delta_{\pm}-n}\tilde{C}_{n}(p^{2})\,e^{-tp^{2}}\,.$
(63)

The Borel transform of $[\tilde{I}(p^{2})]_{\rm OPE}$ is given by

$B_{\tilde{I}}(u)=\sqrt{p^{2}(p^{2}+4m^{2})}\left(\frac{\sqrt{p^{2}+4m^{2}}+\sqrt{p^{2}}}{2\mu}\right)^{-2u}=p^{2}\sum_{n=0}^{\infty}\biggl{(}\frac{m^{2}}{p^{2}}\biggr{)}^{n}\,B_{\tilde{C}_{n}}(u)\,.$
(64)

The Borel transform of each Wilson coefficient is given by

$\displaystyle B_{\tilde{C}_{n}}(u)$
$\displaystyle=\sum_{k=0}^{\infty}\frac{d_{n,k}(L_{p})}{k!}u^{k}=\biggl{(}\frac{\mu^{2}}{p^{2}}\biggr{)}^{u}\frac{(4u^{2}-2u-2n)\Gamma(-n-2u)}{\Gamma(n+1)\Gamma(2-2u-2n)}$
$\displaystyle=\biggl{(}\frac{\mu^{2}}{p^{2}}\biggr{)}^{u}X_{n}(u)\,,$ (65)

where we define an $n$-th degree polynomial of $u$ as

$\displaystyle X_{n}(u)$
$\displaystyle=\sum_{k=0}^{n}\frac{d_{n,k}(0)}{k!}u^{k}=\frac{2(-1)^{n}}{n!}(2u^{2}-u-n)(1+n+2u)_{n-2}\,.$
(66)

$(a)_{j}=\Gamma(a+j)/\Gamma(a)$ denotes the Pochhammer symbol. We confirm that
$B_{\tilde{I}}(u)$ has no singularities on the positive real $u$ axis. (In
fact, it has no singularities in the entire $u$ plane.) Hence,
$[\tilde{I}(p^{2})]_{\rm 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. (65).

The Borel transform of $I_{h}^{\pm}(t)$ is given by

$B_{I_{h}}(u)=\sum_{n=0}^{\infty}(m^{2}t)^{n}\,B_{C_{n}^{(h)}}(u)\,,$ (67)

where the Borel transform of each Wilson coefficient is expressed in terms of
$X_{n}(u)$ as

$B_{C_{n}^{(h)}}(u)=(\mu^{2}t)^{u}\,\Gamma(2-n-u-\delta)\,X_{n}(u)\,.$ (68)

We find that $B_{C_{0}^{(h)}}$ has poles at $u=2,3,4,\cdots$;
$B_{C_{1}^{(h)}}$ at $u=2,3,4,\cdots$; $B_{C_{2}^{(h)}}$ at
$u=1-\delta,2-\delta,3-\delta,\cdots$, and so on. (For $n<2$ we can safely set
$\delta$ to zero.) Thus, $I^{\pm}_{h}(t)$ has IR renormalons. In addition,
expanding $B_{C_{n}^{(h)}}$ in $u$, we see that the perturbative expansion of
$C_{n}^{(h)}$ in $\bar{g}(\mu^{2})$ for $n\geq 2$ includes IR divergences as
multiple poles in $\delta$. This originates from the $u\sim 0$ region, by
expansion of $B_{C_{n}^{(h)}}\sim 1/(u+\delta)$ in $u$.

$I^{\pm}_{h}(t)$ is constructed by Borel resummation as

$I^{\pm}_{h}(t)=\int_{C_{\pm}}du\,B_{I_{h}}(u)\,e^{-u/\bar{g}(\mu^{2})}\,.$
(69)

As already stated, this is equivalent to Eq. (38). For example, we obtain the
LO Wilson coefficient as

$\left[C_{0}^{(h)}(t)\right]_{\pm}=\int_{C_{\pm}}du\,(\mu^{2}t)^{u}\,\Gamma(2-u)\,e^{-u/\bar{g}(\mu^{2})}\,,$
(70)

from Eq. (68) for $n=0$, as given by eq. (36). We can insert
$1/\bar{g}(\mu^{2})=\log(\mu^{2}/m^{2})$, express $\Gamma(2-u)$ as
$\int_{C_{\mp}}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

$\left[C_{0}^{(h)}(t)\right]_{\pm}=t^{2}\int_{C_{\mp}}dp^{2}\,\frac{p^{2}\,e^{-tp^{2}}}{\log(p^{2}/m^{2})}\,.$
(71)

Although IR divergence appears in $\bigl{[}C_{n}^{(h)}\bigr{]}_{\pm}$ for
$n\geq 2$ as multiple poles in $\delta$ at fixed order as mentioned above, by
Borel resummation they turn to $\sim\log\delta$. In the Borel integral, this
originates from the $u\sim 0$ region, where the integrand behaves as $\sim
1/(u+\delta)$. In the dual space integral, the same $\log\delta$ divergence
stems from the $p^{2}\sim 0$ region.

The soft contribution can be expressed as

$I_{s}^{\pm}(t)=\sum_{\ell=0}^{\infty}\frac{t^{\ell-\delta+2}}{\ell!}\,\langle\delta\alpha\,\Box^{\ell-\delta_{\pm}}\delta\alpha\rangle_{\pm}\,,$
(72)

by Eq. (39) and the definition of $D_{\alpha}(p^{2})$. All the matrix elements
have UV divergences,
$\langle\delta\alpha\,\Box^{\ell-\delta_{\pm}}\delta\alpha\rangle_{\pm}\sim\log\delta$,
which stem from the region $p^{2}\to\infty$ and cancel the IR divergences in
$I_{h}^{\pm}$. Thus,

$I_{s}^{\pm}(t)=\sum_{\ell=0}^{\infty}C^{(s)}_{\ell+2}\,t^{\ell+2}\,\langle\delta\alpha\,\Box^{\ell-\delta_{\pm}}\delta\alpha\rangle_{\pm}\,,~{}~{}~{}~{}~{}C^{(s)}_{\ell+2}=\frac{1}{\ell!}\,.$
(73)

The Wilson coefficients of the soft contribution are $t$ independent and
finite as $\delta\to 0$.

Thus, the Wilson coefficients of the hard contribution, $[C_{n}^{(h)}]_{\pm}$,
have imaginary part whose sign depends on the regularization. In addition
$[C_{n}^{(h)}]_{\pm}$ for $n\geq 2$ contains an IR divergence which behaves as
$\log\delta$. This is canceled by the UV divergence contained in the matrix
element
$\langle\delta\alpha\,\Box^{n-2-\delta_{\pm}}\delta\alpha\rangle_{\pm}$. The
cancellation of $\log\delta$ 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$.

###  4.3 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)]_{\rm OPE}$ using a finite number of perturbative
coefficients. The nonperturbative matrix elements need to be determined by
comparison of $[I(t)]_{\rm 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_{n}^{(h)}]_{\pm}$ 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}^{\pm}(t)$ (Eq. (38)). The experimental value
here is the exact $I(t)$ defined by Eqs. (31) and (32), 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
${\cal O}(t^{2})$, and $I_{h}(t)$ is given by

$I_{h}(t)=\big{[}C^{(h)}_{0}(t)\big{]}_{\rm
PV}+\big{[}C^{(h)}_{1}(t)\big{]}_{\rm PV}\langle\alpha\rangle
t+\big{[}C^{(h)}_{2}(t)\big{]}_{\rm PV}^{\rm
ren}\,\langle\alpha\rangle^{2}t^{2}+{\cal O}(t^{3}),$ (74)
$\big{[}C^{(h)}_{0}(t)\big{]}_{\rm PV}=t^{2}\int_{0,\rm
PV}^{\infty}\\!\\!\\!\\!dp^{2}\,p^{2}\sum_{k=0}^{\infty}d_{0,k}(L_{p})\,\bar{g}(\mu^{2})^{k+1}e^{-tp^{2}},$
(75) $\big{[}C^{(h)}_{1}(t)\big{]}_{\rm PV}=t\int_{0,\rm
PV}^{\infty}\\!\\!\\!\\!dp^{2}\,\sum_{k=0}^{\infty}d_{1,k}(L_{p})\,\bar{g}(\mu^{2})^{k+1}e^{-tp^{2}},$
(76) $\big{[}C^{(h)}_{2}(t)\big{]}_{\rm PV}^{\rm ren}=G^{\rm
IR}(t,\Lambda_{\rm IR})+\int_{\Lambda_{\rm IR}^{2},\rm
PV}^{\infty}\\!\\!\frac{dp^{2}}{p^{2}}\,\sum_{k=0}^{\infty}d_{2,k}(L_{p})\,\bar{g}(\mu^{2})^{k+1}e^{-tp^{2}},$
(77) $\displaystyle G^{\rm IR}(t,\Lambda_{\rm IR})$ $\displaystyle=$
$\displaystyle\bigg{[}\int_{0}^{\Lambda_{\rm
IR}^{2}}\\!\\!\frac{dp^{2}}{p^{2+2\delta_{\pm}}}\,\big{[}-2{\bar{g}(p^{2})}-{\bar{g}(p^{2})}^{2}+4{\bar{g}(p^{2})}^{3}\big{]}+2\log(-\delta_{\pm})\bigg{]}_{\delta_{\pm}\to
0}$ (78) $\displaystyle=$ $\displaystyle-2\gamma_{E}-\frac{2}{L_{\rm
IR}^{2}}-\frac{1}{L_{\rm IR}}-2\log(L_{\rm IR}),\quad\quad L_{\rm
IR}=\log(m^{2}/\Lambda_{\rm IR}^{2}).$

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}^{\pm}(t)$ is given analytically by

$I_{h}^{\pm}(t)-I_{h}(t)=-(2\log\delta\,\pm i\pi)m^{4}t^{2}+{\cal
O}(t^{3})\,,$ (79)

which is RG invariant. The Wilson coefficients in the above formulas are
expressed in expansion in $\bar{g}(\mu^{2})$. In the analysis below, we will
set $\mu^{2}\propto t^{-1}$ or $\mu^{2}\propto 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 $\delta$, $\big{[}C^{(h)}_{2}(t)\big{]}_{\rm PV}^{\rm ren}$ is separated
into two contributions by introducing a cutoff $\Lambda_{\rm IR}$, although
the cutoff dependence cancels out between them. The IR divergence behaves as
$\log\delta$, which is quite different from the typical form of IR
divergences. The finite part that remains after subtracting the IR divergence
is RG invariant.

Figure 3: Log-log plot of the difference between $I(t)$ and
$I(t)^{th}=I_{h}(t)+\langle\delta\alpha^{2}\rangle_{\rm PV}^{\rm ren}t^{2}$.
The black (red) line corresponds to $I(t)^{th}$ computed up to ${\cal
O}(t^{1})$ (${\cal O}(t^{2})$). The slope of each line represents the power
behavior of $t$.

The soft contribution $I_{s}(t)$ can also be computed in this model. Its
expression up to ${\cal O}(t^{2})$ is obtained as

$I_{s}(t)^{\pm}=\langle\delta\alpha^{2}\rangle_{\rm PV}^{\rm
ren}\,t^{2}+(2\log\delta\pm i\pi)\,m^{4}t^{2}+{\cal O}(t^{3}),$ (80)
$\langle\delta\alpha^{2}\rangle_{\rm PV}^{\rm ren}=\int_{0}^{\Lambda_{\rm
UV}^{2}}\\!\\!\\!\\!dp^{2}\tilde{I}(p^{2})+G^{\rm UV}(\Lambda_{\rm UV}),$ (81)
$\displaystyle G^{\rm UV}(\Lambda_{\rm UV})$ $\displaystyle=$
$\displaystyle\bigg{[}\int_{\Lambda_{\rm
UV}^{2}}^{\infty}\\!\\!{dp^{2}}{p^{2-2\delta_{\pm}}}\,\sum_{n=0}^{2}(m^{2}/p^{2})^{n}\sum_{k=0}^{\infty}d_{n,k}(0)\bar{g}(p^{2})^{k+1}-2\log(\delta_{\pm})\bigg{]}_{\delta\to
0}$ (82) $\displaystyle=$ $\displaystyle 2\gamma_{E}-{\rm Ei}\big{(}2L_{\rm
UV}\big{)}+\frac{2-L_{\rm UV}\big{(}1+2(\Lambda_{\rm
UV}^{2}/m^{2})\big{)}}{L_{\rm UV}^{2}}+2\log\big{(}L_{\rm UV}\big{)},$

where ${\rm Ei}(z)=-\int_{-z}^{\infty}dt\,e^{-t}/t$, and

$L_{\rm UV}=\log(\Lambda_{\rm UV}^{2}/m^{2}).$ (83)

To explicitly compute the UV divergence, we introduce a cutoff $\Lambda_{\rm
UV}$, although $\langle\delta\alpha^{2}\rangle_{\rm PV}^{\rm ren}$ also has no
overall cutoff dependence. Comparing Eqs. (79) and (80), we can see that the
IR and UV divergences and the imaginary part cancel explicitly.
$\langle\delta\alpha^{2}\rangle_{\rm PV}^{\rm ren}$ can be evaluated by
numerical integration to be

$\frac{\langle\delta\alpha^{2}\rangle_{\rm PV}^{\rm
ren}}{m^{4}}=-1.386\cdots\,.$ (84)

As a cross check, we compare the exact (experimental) value of $I(t)$ and
$I(t)^{th}=I_{h}(t)+\langle\delta\alpha^{2}\rangle_{\rm PV}^{\rm ren}\,t^{2}$
using $\langle\alpha\rangle=m^{2}$ and Eq. (84). The difference
$I(t)^{th}-I(t)$ is shown in Fig. 3. We set $\mu^{2}=p^{2}$ in the integrand
(in the dual space) and calculate $I(t)^{th}$ up to ${\cal O}(t^{1})$ and
${\cal O}(t^{2})$, which correspond to the black and red lines, respectively.
The typical size of the coupling constant in this figure is
$\bar{g}(1/\sqrt{t})\in[0.1,\,0.2]$. We see that the exact value of $I(t)$ can
be predicted systematically, by increasing the expansion order of $t$ for
$I_{h}^{\pm}(t)$ and $I_{s}^{\pm}(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, $\langle\alpha\rangle=m^{2}$, and
$\langle\delta\alpha^{2}\rangle_{\rm PV}^{\rm ren}$ given by Eq. (81), 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

$I_{\rm exp}(t_{\rm exp};r)=t_{\rm
exp}^{2}\int_{0}^{\infty}dp^{2}\tilde{I}(p^{2})\,e^{-t_{\rm
exp}p^{2}},\quad\quad t_{\rm exp}=r^{2}\,t.$ (85)

$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
by131313 In 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 $\sim
1/(u+\delta)$ in $B_{C_{2}^{(h)}}$ of Eq. (68), and the Borel integral of
$1/(u+\delta)$ gives the $\log\delta$ singularity. Hence, we subtract the
contribution from $B_{C_{2}^{(h)}}$ corresponding to $\log\delta$, expand in
$u$ and take the Borel integral. We adjust the finite ($\delta$ independent)
term such that the same divergence is subtracted as in the DSRS method.

$\displaystyle I_{\rm fit}^{\rm RG}(t,s,\\{k_{i}\\};A,B)=$
$\displaystyle\big{[}C^{(h)}_{0}(t)\big{]}_{\rm
RG}^{(k_{0})}+\big{[}C^{(h)}_{1}(t)\big{]}_{\rm RG}^{(k_{1})}Am^{2}t$
$\displaystyle+\Big{[}\big{[}C^{(h)}_{2}(t)\big{]}_{\rm
RG}^{(k_{2})}+2\log\delta\Big{]}A^{2}m^{4}t^{2}+Bm^{4}t^{2},$ (86)

where

$\big{[}C^{(h)}_{n}(t)\big{]}_{\rm
RG}^{(k_{n})}=\big{[}C^{(h)}_{n}(t)\big{]}_{\rm PT}\bigg{|}_{{\rm
up\,to\,}{\cal O}(\bar{g}(\mu^{2})^{k_{n}+1})}\,{\rm
with}~{}~{}~{}\mu^{2}=s^{2}/t\,,$ (87)

and $A={\langle{\alpha}\rangle}/{m^{2}}$ and
$B={\langle\delta{\alpha}^{2}\rangle_{\rm PV}^{\rm ren}}/{m^{4}}$ stand for
the fit parameters. $[C_{n}^{(h)}]_{\rm PT}$ denotes the (naive) perturbative
expansion of $C_{n}^{(h)}$ in $\bar{g}(\mu^{2})$.

The second one is calculated using the DSRS method and given by

$\displaystyle I_{\rm fit}^{\rm PV}(t,s,\\{k_{n}\\};A,B)=$
$\displaystyle\big{[}C^{(h)}_{0}(t)\big{]}_{\rm
PV}^{(k_{0})}+\big{[}C^{(h)}_{1}(t)\big{]}_{\rm PV}^{(k_{1})}Am^{2}t$
$\displaystyle+\big{[}C^{(h)}_{2}(t)\big{]}_{\rm
PV}^{(k_{2})}A^{2}m^{4}t^{2}+Bm^{4}t^{2}\,,$ (88)

where we truncate the integrands of Eqs. (75)–(77) at ${\cal
O}(\bar{g}(\mu^{2})^{k_{n}+1})$ with

$\mu^{2}=s^{2}p^{2}.$ (89)

(We treat $G^{\rm IR}$ exactly.) In both formulas, we subtract the IR
divergence of the ${\cal O}(t^{2})$ term in the OPE, and each Wilson
coefficient is calculated by truncation of the perturbative expansion.

In Eq. (86), the Wilson coefficients contain renormalons which cause
uncertainties in the fit. The major uncertainty comes from the $u=2$
renormalon in $\big{[}C^{(h)}_{0}(t)\big{]}_{\rm RG}$, which would give an
order $100$ per cent uncertainty in the determination accuracy of $B$. On the
other hand, in Eq. (88), 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\cdots\,,$ (90)

accurately by the fit.

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.

The details of the fits are as follows. The data points used for the fits
correspond to $\hat{t}\equiv m^{2}t={n}/{10^{4}}$ for $n=1,\,2,\,\cdots,\,10$.
The maximal size of the coupling constant in this region is
$\bar{g}={1}/{\log(10^{3})}\approx 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. (86) and (88), and the systematic uncertainty is assigned from the
difference as we vary $s\in\\{1/2,1,2\\}$. In this analysis, the truncation
errors of perturbative expansions of the Wilson coefficients $C^{(h)}_{n}$ for
$n\geq 1$ are assumed to be negligible. Hence, we fix the truncation orders of
$C^{(h)}_{1}$ and $C^{(h)}_{2}$ at $k_{1}=3$ and $k_{2}=3$, respectively, and
fix the scale at $s=1$ in both calculations, for simplicity.141414 In the DSRS
method, it is equivalent to fixing the truncation orders at $k_{1}=1$ and
$k_{2}=2$, since $d_{1,k}(0)=0$ for $k\geq 2$ and $d_{2,k}(0)=0$ for $k\geq
3$.  In the analysis with Eq. (88), we assign another systematic uncertainty
from (supposedly) unknown ${\cal 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_{\rm fit}^{\rm
PV}$ and $I_{\rm fit}^{\rm PV}+Cm^{6}t^{3}$, in which $C$ is an additional fit
parameter. This effect is negligible for the analysis with Eq. (86) 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}\sim 13$ for
all the parameters, but the uncertainty increases beyond that order. Since the
$u=2$ renormalon is included in $C^{(h)}_{0}$, we expect that the perturbative
series converges up to $k_{0}\approx k_{0}^{*}\equiv 2/\bar{g}\approx 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\approx
k_{0}^{*}$ are given by

$r-1=4(6)\times 10^{-6},\quad A-1=5(6)\times 10^{-3},\quad B=-3.1(2.9),$ (91)

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\approx k_{0}^{*}$ are given by

$r-1=4(4)\times 10^{-8},\quad A-1=4(4)\times 10^{-5},\quad B=-1.406(19),$ (92)

which are consistent with the exact values of the parameters, Eq. (90). We can
see that the error sizes of $r$ and $A$ in Eq. (92) are orders of magnitude
smaller than those in Eq. (91). This is because the major source of systematic
uncertainties in Eq. (92) is the ${\cal O}(t^{3})$ OPE corrections, while that
in Eq. (91) originates from the $u=2$ renormalon in
$\big{[}C^{(h)}_{0}\big{]}_{\rm RG}$ which is order ${\cal O}(t^{2})$. 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.

###  4.4 Self-energy of $\sigma$

The inverse propagator of $\sigma$ with momentum $q^{\mu}$ is given by
$(q^{2}+m^{2})$ at LO of the $1/N$ expansion. The second example we consider
is the self-energy of $\sigma$ at ${\cal O}(1/N)$, denoted as
$\Pi_{\sigma}(q^{2})$. The self-energy is included universally in calculations
of correlation functions of $\sigma$ and is given by the diagram in Fig. 5.
(We omit the tadpole diagram which is absorbed by renormalization of
$\langle\alpha\rangle$.) Cancellation of the renormalons in the entire OPE of
the self-energy has been shown in ref. Beneke:1998eq . 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 $\Pi_{\sigma}(q^{2})$, we consider
the second derivative of the self-energy,

$\rho(q^{2})\equiv
Nq^{2}\,\Pi_{\sigma}^{\prime\prime}(q^{2})=q^{2}\int\frac{d^{2}k}{(2\pi)^{2}}\left(\frac{\partial}{\partial
q^{2}}\right)^{2}\frac{1}{(q+k)^{2}+m^{2}}\,N\,D_{\alpha}(k^{2})\,.$ (93)

It is defined to be dimensionless and RG invariant. The above integral is well
defined nonperturbatively, without UV or IR divergence.

Figure 5: Self-energy diagram of $\sigma$ at ${\cal O}(1/N)$. The tadpole
diagram is omitted.

We can construct the OPE of $\rho(q^{2})$ in the dual space approach similarly
to the previous example. We define the dual space observable
$\tilde{\rho}(p^{2})$ by

$\rho(q^{2})=\left(\frac{1}{q^{2}}\right)^{2}\int_{0}^{\infty}dp^{2}\,\tilde{\rho}(p^{2})\,e^{-p^{2}/q^{2}}\,.$
(94)

Using the generalized EBR, we can express the OPE of $\rho(q^{2})$ as the sum
of the hard and soft contributions:

$\displaystyle\big{[}\rho(q^{2})\big{]}_{\rm OPE}=\lim_{\delta\to
0}\left\\{\rho^{\pm}_{h}(q^{2})+\rho^{\pm}_{s}(q^{2})\right\\}\,,$ (95)
$\displaystyle\rho^{\pm}_{h}(q^{2})=\left(\frac{1}{q^{2}}\right)^{2-\delta}\int_{C_{\mp}}dp^{2}\,(p^{2})^{-\delta_{\pm}}\left[\tilde{\rho}(p^{2})\right]_{\rm
OPE}\,e^{-p^{2}/q^{2}}\,,$ (96)
$\displaystyle\rho^{\pm}_{s}(q^{2})=\left(\frac{1}{q^{2}}\right)^{2-\delta}\int_{C_{\mp}}dp^{2}\,(p^{2})^{-\delta_{\pm}}\tilde{\rho}(p^{2})\,\left\\{1-\frac{p^{2}}{q^{2}}+\frac{1}{2}\left(\frac{p^{2}}{q^{2}}\right)^{2}-\dots\right\\}\,,$
(97)

where the integrals are regularized in the same way as in the previous
example.

We can check that the Wilson coefficients of $[\tilde{\rho}(p^{2})]_{\rm OPE}$
do not have IR renormalons or IR divergences. In fact, the Borel transform of
$[\tilde{\rho}(p^{2})]_{\rm OPE}$ is given by

$\displaystyle B_{\tilde{\rho}}=$
$\displaystyle-\frac{p^{2}}{\pi}\left(\frac{\mu^{2}}{p^{2}}\right)^{u}\sum_{M=0}^{\infty}\left(\frac{m^{2}}{p^{2}}\right)^{M}\Gamma(M+u-1)$
$\displaystyle\times\sum_{n=0}^{M}\left(\frac{(n+u-1)_{M-n}}{(M-n)!}\right)^{2}\,(-1)^{n}\,X_{n}(u)\left\\{\sin(\pi
u)\,\tilde{f}_{M,n}(p^{2})+\cos(\pi u)\,\tilde{g}_{M}\right\\}\,,$ (98)

where

$\displaystyle\tilde{f}_{M,n}(p^{2})=1-2(M\\!+\\!u)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+(M\\!+\\!u)(M\\!+\\!u\\!-\\!1)\left\\{2\,\psi(M\\!-\\!n\\!+\\!1)-3\,\psi(M\\!+\\!u\\!-\\!1)+\log\\!\left(\frac{p^{2}}{m^{2}}\right)\right\\}\,,$
(99)
$\displaystyle\tilde{g}_{M}~{}~{}=~{}~{}-2\pi(M\\!+\\!u)(M\\!+\\!u\\!-\\!1)\,,$
(100)

and $\psi(z)=\Gamma^{\,\prime}(z)/\Gamma(z)$ denotes the digamma function.
(See App. D for the derivation.) $B_{\tilde{\rho}}$ has no poles on the
positive $u$ axis, hence $[\tilde{\rho}(p^{2})]_{\rm OPE}$ has no IR
renormalons. (The poles at $u=0,1$ of the gamma and digamma functions for
$M\leq 1$ are canceled by other factors in the summand.) On the other hand,
$B_{\tilde{\rho}}$ has poles on the negative $u$ axis, hence
$[\tilde{\rho}(p^{2})]_{\rm OPE}$ has UV renormalons. The closest pole is
located at $u=-1$. 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_{\tilde{\rho}}$ has no pole at $u=0$ and is independent of any
regularization parameter. Hence, $[\tilde{\rho}(p^{2})]_{\rm OPE}$ is also
free from IR divergences. This aspect is the same as that of the previous
example $\tilde{I}(p^{2})$.

We can calculate the Borel transform of $\rho_{h}^{\pm}(q^{2})$ by taking the
inverse dual transform of $B_{\tilde{\rho}}$.

$\displaystyle B_{\rho_{h}}=$
$\displaystyle-\frac{1}{\pi}\left(\frac{\mu^{2}}{q^{2}}\right)^{u}\sum_{M=0}^{\infty}\left(-\frac{m^{2}}{q^{2}}\right)^{M}\Gamma(M+u-1)$
$\displaystyle\times\sum_{n=0}^{M}\frac{(n+u-1)_{M-n}}{\\{(M-n)!\\}^{2}}\,X_{n}(u)\left\\{\sin(\pi
u)\,{f}_{M,n}(q^{2})+\cos(\pi u)\,{g}_{M}\right\\}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\times\Gamma(2-n-u-\delta)\,,$ (101)

where

$\displaystyle{f}_{M,n}(q^{2})=\tilde{f}_{M,n}(q^{2})+(M\\!+\\!u)(M\\!+\\!u\\!-\\!1)\,\psi(M\\!+\\!u\\!-\\!1)\,,$
(102)
$\displaystyle{g}_{M}~{}~{}~{}=~{}~{}~{}-\pi(M\\!+\\!u)(M\\!+\\!u\\!-\\!1)\,.$
(103)

IR renormalons are generated by the inverse dual transform. All of them stem
from the factor $\Gamma(2\\!-\\!n\\!-\\!u\\!-\\!\delta)$, just as in the
previous example, cf., Eq. (68). The generated IR renormalons have only
integer powers of $m^{2}/q^{2}$ and no logarithmic corrections are included.
This is because, the coefficient of $\log(q^{2}/m^{2})$ term in Eq. (102) has
no poles of $u$ due to $\sin(\pi u)$ factor. Hence, (a posteori) the simple
dual transform Eq. (94) is sufficient to remove all the renormalons. The same
factor $\Gamma(2\\!-\\!n\\!-\\!u\\!-\\!\delta)$ also generates IR divergences
$\sim\log\delta$ for $M\geq 2$ by Borel integral. The above formulas can be
used to obtain the perturbative series for the Wilson coefficients of
$[\tilde{\rho}(p^{2})]_{\rm OPE}$ and $[{\rho}(q^{2})]_{\rm OPE}$.

Since $[\tilde{\rho}(p^{2})]_{\rm OPE}$ does not have IR renormalons, the dual
space integral expression for $\rho^{\pm}_{h}(q^{2})$ converges well to its
exact value, as we include more terms of the series expansion in $m^{2}/p^{2}$
and $\bar{g}(p^{2})$ (after resummation of the UV renormalons). Using the dual
space expressions for $\rho^{\pm}_{h}(q^{2})$ and $\rho^{\pm}_{s}(q^{2})$, the
cancellation of the IR and UV divergences, as well as of the imaginary part,
can be shown in a similar manner to the previous example. This is because, for
instance, the UV behavior of the integrand of the soft contribution
$\rho^{\pm}_{s}(q^{2})$ is determined by $[\tilde{\rho}(p^{2})]_{\rm OPE}$,
and hence, the general term necessary for the proof is of the form
$(p^{2})^{n-\delta_{\pm}}\bar{g}(p^{2})^{k}=(p^{2})^{n-\delta_{\pm}}/\log(p^{2}/m^{2})^{k}$
and the identity Eq. (52) can be used. Furthermore, it is possible to perform
a simulation study and confirm good convergence of the OPE to the exact
$\rho(q^{2})$ after renormalon subtraction by the DSRS method. We perform a
similar fit as in Sec. 4.3 and obtain values of the nonperturbative parameters
which agree with the exact values with good accuracies. Thus, apart from the
resummation of the UV renormalons, an analysis similar to the previous example
can be made.

Finally we compare our study with the preceding analyses on the OPEs in the
$O(N)$ nonlinear $\sigma$ model. Refs. Bardeen:1976zh ; David:1982qv ;
David:1985xj  studied renormalon cancellation in the OPE of
$\langle\delta\alpha^{2}\rangle$ regularized by a cutoff, which mainly focused
on the LO Wilson coefficient and the ${\cal O}(m^{4})$ matrix element. Ref.
Beneke:1998eq  studied renormalon cancellation in the OPE of the renormalized
$\sigma$ self-energy. It showed the renormalon cancellation as well as
cancellation between IR and UV divergences in the entire OPE in the Borel
space. Ref. Lee:2001ws  investigated an application of the cancellation of the
imaginary part calculated in ref. Beneke:1998eq  (up to ${\cal O}(m^{2})$).
Novel aspects of the present paper are as follows. First, we propose the new
concept “dual space approach,” 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. (52) and (56)], 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.

##  5 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 $\sigma$ 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 $\Lambda_{\rm QCD}$. (B) How to incorporate the effects of
the beta function beyond one loop order ($b_{1},\,b_{2},\,\cdots\neq 0$).

(A) Renormalons are generally not in the simple power form of $\Lambda_{\rm
QCD}/Q$ but include corrections in the form of the anomalous dimension and
expansion in $\alpha_{s}(Q^{2})$, as given by Eq. (11). A method for
incorporating these corrections is given in App. B of ref. Hayashi:2021vdq  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},\cdots\neq 0$ in the dual space approach is
straightforward, and how this is done is already explained in Sec. 3 in
principle. Note that $\Lambda_{\rm\overline{MS}}$ in e.g. Eq. (11) 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_{\pm}$.
(ii) IR divergences occur from the integral over $p^{2}$ from the IR region
$p^{2}\sim 0$, beyond a certain order in expansion in $1/p$. By regularizing
the $p^{2}$ integral similarly to the nonlinear $\sigma$ model, the IR
divergences appear as $\sim\log\delta$ . This is the case, even including NLL
and beyond, since the coupling constant behaves as

$\alpha_{s}(p^{2a})=\frac{1}{b_{0}\log(p^{2a}/\Lambda_{\rm\overline{MS}}^{2})}-\frac{b_{1}}{b_{0}^{3}}\frac{\log\log(p^{2a}/\Lambda_{\rm\overline{MS}}^{2})}{\log(p^{2a}/\Lambda_{\rm\overline{MS}}^{2})^{2}}+{\cal
O}\left(\log(p^{2a}/\Lambda_{\rm\overline{MS}}^{2})^{-2}\right)\,,$ (104)

as $p^{2}\to 0$, and only the leading term causes the divergence. For instance
we can subtract the $\log\delta$ 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
$e^{-tp^{2}}=1-tp^{2}+\cdots$, while the hard expansion generates the OPE in
the dual space, $[\tilde{S}(p^{2})]_{\rm OPE}$, which is the double expansion
in $1/p^{2}$ and $\alpha_{s}(p^{2a})$. See Eq. (56). In the case
$b_{1},\,b_{2},\,\cdots\neq 0$, the generalized EBR method is guaranteed by
the equality

$\int_{C_{\pm}}dp^{2}\,\frac{(p^{2})^{a_{\pm}}}{\\{\log(p^{2}/m^{2})\\}^{n}}\,[\log\log(p^{2}/m^{2})]^{k}=0,~{}~{}~{}~{}~{}(a_{\pm}\in\mathbb{C},~{}n,k\in\mathbb{Z}_{>0})$
(105)

since the high energy asymptotic expansion of $\alpha_{s}(p^{2a})$ generates
such a term.

(2) In computing Wilson coefficients by perturbative QCD, it is customary to
use dimensional regularization. Hence, instead of the pseudo-dimensional
regularization we used in Sec. 4 for simplicity, we can use the full
dimensional regularization. This would help defining the nonperturbative
matrix elements universally, i.e., independently of observables, in a fully
consistent manner.

(3) In the dual space approach, we justify reconstruction of an observable by
the OPE using the EBR relation Eq. (105). It means that the IR divergence
should be treated in the $\log\delta$ (resummed) form rather than the
$1/\delta^{n}$ form which originates from the expansion of $1/(u+\delta)$ in
$u$ (in the Borel integral expression). In particular only in the $\log\delta$
form can we show the cancellation of the imaginary part. If we subtract IR
divergences in the expanded form, we are concerned about the ${\cal
O}(\delta^{0})$ (finite) term generated by the product of $1/\delta^{n}$ and
$\delta^{n}$. Such terms remain and can be dependent on $Q^{2}$, which would
give contributions different from the resummed prescription.

Figure 6: Example of a diagram which includes two $\delta\alpha$ propagators.

Finally we discuss an unsolved problem. In (A) above, we explained how to
incorporate corrections to the simple power form $(\Lambda_{\rm QCD}/Q)^{n}$
of renormalons into the DSRS method. We anticipate that the sum of these
corrections given as series expansion in $\alpha_{s}(Q^{2})$ possibly
diverges. For example, this type of corrections are expected to be included in
the diagram of Fig. 6 in the nonlinear $\sigma$ model, which contributes at
${\cal O}(1/N^{2})$. This diagram includes two $\delta\alpha$ propagators,
say, $D_{\alpha}(k_{1}^{2})$ and $D_{\alpha}(k_{2}^{2})$. From the region
where $k_{1}$ is hard and $k_{2}$ is soft, $k_{1}\gg k_{2}\sim m$, this
diagram generates the product of the matrix element
$\langle(\delta\alpha)^{2}\rangle$, which originates from
$D_{\alpha}(k_{2}^{2})$, and Wilson coefficients $C_{i}$ corresponding to the
contribution from $[D_{\alpha}(k_{1}^{2})]_{\rm OPE}$. The Wilson coefficients
$C_{i}$ include IR renormalons. Thus, the LO Wilson coefficient, corresponding
to the (hard)-(hard) region of the diagram, is expected to include a
renormalon of the form $C_{i}\langle(\delta\alpha)^{2}\rangle$, where $C_{i}$
is given by a divergent series. Although this problem is not specific to the
dual space approach, up to now we do not know how to remove a renormalon of
this type using the DSRS method. Stepwise inclusion of the corrections, first
through the anomalous dimension and then order by order in $\alpha_{s}(Q^{2})$
expansion, is only an asymptotic approximation in this case.

##  6 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
$\tilde{S}(p^{2})$ such that the OPE of $\tilde{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
$\alpha_{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 $\Lambda_{\rm 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.

The above properties of the DSRS method are derived under the following
assumptions. We assumed that (i) the OPE of $S(Q^{2})$ or $\tilde{S}(p^{2})$
coincides with its high-energy expansion, (ii) for each Wilson coefficient, we
can neglect contributions from singularities other than renormalon
singularities in the right-half Borel plane, and (iii) in the case of QCD, the
beta function beyond five loops does not alter the analytic structure of the
beta function known up to five loops Hayashi:2021vdq .

In particular, the EBR method is applied to the $p^{2}$ integral (inverse dual
transform). We emphasize that it constitutes the basis of our analysis of the
OPE.

  * (1)

We applied the EBR method to the Laplace ($p^{2}$) integral representation of
$S(Q^{2})$, which generates the OPE of $S(Q^{2})$. The one-parameter integral
form allows a simple proof of the EBR method.

  * (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 $\sigma$ 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. (105), 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 Hayashi:2021vdq . 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 $(\Lambda_{\rm QCD}/Q)^{n}$ of each
IR renormalon is incorporated by anomalous dimension and expansion in
$\alpha_{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.

## Acknowledgments

Y.H. acknowledges support from GP-PU at Tohoku University. This work was
supported by JSPS KAKENHI Grant Numbers JP21J10226, JP20J00328, JP20K03923,
JP19K14711, and JP18H05542.

##  Appendix A Relation to FTRS method

In this appendix we show that after resummation of artificial UV renormalons
in the FTRS method, it agrees with the DSRS method.

Let us consider a dimensionless and RG invariant observable $S(Q^{2})$. The
starting point of the FTRS method is to define the Fourier transform of
$S(Q^{2})$ as

$\tilde{S}^{\prime}(\tau)=\int
d^{3}\vec{x}\,e^{-i\vec{\tau}\cdot\vec{x}}r^{2au^{\prime}}S(r^{-2a})\,,~{}~{}~{}~{}~{}Q=r^{-a}\,,~{}~{}~{}|\vec{x}|=r\,,~{}~{}~{}|\vec{\tau}|=\tau\,.$
(106)

This is applied to the perturbation expansion
$S(Q^{2})=\sum_{n=0}^{\infty}c_{n}(L_{Q})\alpha_{s}(\mu^{2})^{n+1}$. Using

$S(r^{-2a})=\big{(}\mu^{2}r^{2a}\big{)}^{\hat{H}}\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\mu^{2})^{n+1},$
(107)

as the formal solution to the RG equation, with the operator $\hat{H}$ defined
in Eq. (19), we derive

$\displaystyle\tilde{S}^{\prime}(\tau)$ $\displaystyle=$ $\displaystyle\int
d^{3}\vec{x}\,e^{-i\vec{\tau}\cdot\vec{x}}r^{2au^{\prime}}\big{(}\mu^{2}r^{2a}\big{)}^{\hat{H}}\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\mu^{2})^{n+1}$
$\displaystyle=$
$\displaystyle-\frac{4\pi}{\tau^{3+2au^{\prime}}}\big{(}\mu^{2}/\tau^{2a}\big{)}^{\hat{H}}\sin(\pi
a(\hat{H}+u^{\prime}))\Gamma(2a(\hat{H}+u^{\prime})+2)\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\mu^{2})^{n+1}.$

This expression shows that in the Fourier transformed series, renormalons at
$u=-u^{\prime}+\frac{n}{a}$ are suppressed (one can see this by a similar
computation to Eqs. (15) and (16)), while renormalons at
$u=-u^{\prime}-\frac{2n+1}{2a}$ are newly created, where
$n=0,\,1,\,2,\,\cdots$. With a natural parameter choice, $a>0,\,u^{\prime}<0$,
IR renormalons are suppressed and UV renormalons are generated. The latter
renormalons are called “artificial UV renormalons.”

First we note the identity to resum these aritificial UV renormalons:

$\displaystyle\sin(\pi x)\Gamma(2x+2)$ $\displaystyle=$
$\displaystyle-\frac{\pi\Gamma(2x+2)}{\Gamma(-x)\Gamma(1+x)}=-\frac{\sqrt{\pi}}{\Gamma(-x)}2^{1+2x}\Gamma(x+3/2)$
(109) $\displaystyle=$
$\displaystyle-\frac{\sqrt{\pi}}{\Gamma(-x)}2^{1+2x}\int_{0}^{\infty}dy\,y^{x+1/2}\,e^{-y}\,.$

Using this identity we rewrite Eq. (LABEL:AppA-FT1) as

$\displaystyle\tilde{S}^{\prime}(\tau)$ $\displaystyle=$
$\displaystyle\frac{4\pi^{3/2}}{\tau^{3+2au^{\prime}}}\big{(}\mu^{2}/\tau^{2a}\big{)}^{\hat{H}}\frac{2^{1+2a(\hat{H}+u^{\prime})}}{\Gamma(-a(\hat{H}+u^{\prime}))}\int_{0}^{\infty}\\!\\!\\!dy\,y^{a(\hat{H}+u^{\prime})+1/2}\,e^{-y}\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\mu^{2})^{n+1}$
(110) $\displaystyle=$
$\displaystyle\frac{4\pi^{3/2}}{\tau^{3+2au^{\prime}}}\big{(}\mu^{2}/\tau^{2a}\big{)}^{\hat{H}}\frac{2^{1+2au^{\prime}}}{\Gamma(-a(\hat{H}+u^{\prime}))}\int_{0}^{\infty}\\!\\!\\!dy\,y^{au^{\prime}+1/2}\,e^{-y}\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\mu^{2}/(4y)^{a})^{n+1}$
$\displaystyle=$
$\displaystyle\frac{4\pi^{3/2}}{\tau^{3+2au^{\prime}}}\int_{0}^{\infty}\\!\\!\\!dy\,y^{au^{\prime}+1/2}\,e^{-y}\frac{2^{1+2au^{\prime}}}{\Gamma(-a(\hat{H}+u^{\prime}))}\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\tau^{2a}/(4y)^{a})^{n+1}\,.$

Here, we used the RG invariance of $S$.

The regularized $S(Q^{2})$ is obtained by the inverse Fourier transform of
$\tilde{S}^{\prime}(\tau)$ with contour deformation. Hence, the principal
value integral that eliminates renormalons are given by

$\displaystyle\big{[}S(Q^{2})\big{]}_{\rm PV}$ $\displaystyle=$
$\displaystyle\frac{r^{-2au^{\prime}-1}}{2\pi^{2}}\int_{0,\rm
PV}^{\infty}d\tau\,\tau\sin(\tau r)\tilde{S}^{\prime}(\tau)$ (111)
$\displaystyle=$
$\displaystyle\frac{r^{-2au^{\prime}-1}}{\sqrt{\pi}}\int_{0,\rm
PV}^{\infty}\frac{d\tau}{\tau^{2+2au^{\prime}}}\,\sin(\tau
r)\int_{0}^{\infty}dy\,y^{au^{\prime}+1/2}\,e^{-y}\frac{2^{2+2au^{\prime}}}{\Gamma(-a(\hat{H}+u^{\prime}))}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\tau^{2a}/(4y)^{a})^{n+1}$
$\displaystyle=$
$\displaystyle\frac{2r^{-2au^{\prime}-1}}{\sqrt{\pi}}\int_{0}^{\infty}dy\,\,e^{-y}\int_{0,\rm
PV}^{\infty}\frac{d\tau}{\tau^{2+2au^{\prime}}}\,\frac{\sin(2\sqrt{y}\tau
r)}{\Gamma(-a(\hat{H}+u^{\prime}))}\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\tau^{2a})^{n+1}$
$\displaystyle=$ $\displaystyle{2r^{-2au^{\prime}}}\int_{0,\rm
PV}^{\infty}\frac{d\tau}{\tau^{1+2au^{\prime}}}\,e^{-\tau^{2}r^{2}}\frac{1}{\Gamma(-a(\hat{H}+u^{\prime}))}\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(\tau^{2a})^{n+1}$
$\displaystyle=$ $\displaystyle t^{-au^{\prime}}\int_{0,{\rm
PV}}^{\infty}dp^{2}\,e^{-tp^{2}}\frac{1}{(p^{2})^{1+au^{\prime}}}\frac{1}{\Gamma(-a(\hat{H}+u^{\prime}))}\sum_{n=0}^{\infty}c_{n}(0)\alpha_{s}(p^{2a})^{n+1}\,.$

In the third equality the variable of integration is changed, and in the
fourth equality we perform the $y$ integration. Finally, we define
$\tau=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. (21) and (24).

##  Appendix 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 $\bar{u}$ to control the resummation.
Let us take the LO Wilson coefficient in the dual space $\tilde{C}_{0}(p^{2})$
and explain how to resum renormalons. $\tilde{C}_{0}(p^{2})$ can be expressed
by the following one-parameter integral form,

$\displaystyle\tilde{C}_{0}(p^{2})$ $\displaystyle=$
$\displaystyle\frac{1}{(p^{2})^{1+au^{\prime}}}\sum_{n=0}^{\infty}\alpha_{s}(\mu^{2})^{n+1}\tilde{c}_{n}(L_{p})$
(112) $\displaystyle=$
$\displaystyle\frac{1}{(p^{2})^{1+au^{\prime}}}\frac{\Gamma(a(\bar{u}+\hat{H}))}{\Gamma(a(\bar{u}+\hat{H}))}\sum_{n=0}^{\infty}\alpha_{s}(\mu^{2})^{n+1}\tilde{c}_{n}(L_{p})$
$\displaystyle=$
$\displaystyle\frac{1}{(p^{2})^{1+au^{\prime}}}\int_{0}^{\infty}d\zeta\,\zeta^{a(\bar{u}+\hat{H})-1}e^{-\zeta}\,\sum_{n=0}^{\infty}\alpha_{s}(\mu^{2})^{n+1}\tilde{c}^{\prime}_{n}(L_{p})$
$\displaystyle=$
$\displaystyle\frac{1}{(p^{2})^{1+au^{\prime}}}\int_{0}^{\infty}d\zeta\,\zeta^{a\bar{u}-1}e^{-\zeta}\,\sum_{n=0}^{\infty}\alpha_{s}(p^{2a}\zeta^{-a})^{n+1}\tilde{c}^{\prime}_{n}(0)\,.$

Here, $\tilde{c}^{\prime}_{n}$’s can be determined from $\tilde{c}_{n}$ or
from the perturbative series in the original space, Eq. (8), by the following
relation

$\displaystyle\sum_{n=0}^{\infty}\alpha_{s}(\mu^{2})^{n+1}\tilde{c}^{\prime}_{n}(0)$
$\displaystyle=\frac{1}{\Gamma(a(\bar{u}+\hat{H}))}\sum_{n=0}^{\infty}\alpha_{s}(\mu^{2})^{n+1}\tilde{c}_{n}(0)$
$\displaystyle=\frac{1}{\Gamma(a(\bar{u}+\hat{H}))\Gamma(-a(u^{\prime}+\hat{H}))}\sum_{n=0}^{\infty}\alpha_{s}(\mu^{2})^{n+1}c_{n}(0)\,.$
(113)

See Eqs. (18)–(20). In the last line of Eq. (112), after applying
$\zeta^{a\hat{H}}$ to $\alpha_{s}(\mu^{2})^{n+1}$, we have set
$\mu^{2}=p^{2}$. Note the RG running formula

$(\mu_{1}^{2}/\mu_{2}^{2})^{\hat{H}}\alpha_{s}(\mu_{1}^{2})^{n+1}=\alpha_{s}(\mu_{2}^{2})^{n+1}\,.$
(114)

It can be seen that the UV renormalons of
$\sum_{n=0}^{\infty}\alpha_{s}(\mu^{2})^{n+1}\tilde{c}_{n}(0)$ corresponding
to $u=-\bar{u}$, $-\bar{u}-1/a$, $-\bar{u}-2/a$, $\cdots$ are suppressed by
the $1/\Gamma(a(\bar{u}+\hat{H}))$ factor.

Then $\big{[}C_{0}\big{]}_{\pm}$ with resummation of UV renormalons is
calculated by

$\displaystyle\big{[}C_{0}(Q^{2})\big{]}_{\pm}$ $\displaystyle=$
$\displaystyle
Q^{2u^{\prime}}\int_{C_{\mp}}dp^{2}\,e^{-tp^{2}}\tilde{C}_{0}(p^{2})$ (115)
$\displaystyle=$ $\displaystyle
Q^{2u^{\prime}}\int_{C_{\mp}}dp^{2}\,e^{-tp^{2}}\frac{1}{(p^{2})^{1+au^{\prime}}}\int_{0}^{\infty}d\zeta\,\zeta^{a\bar{u}-1}e^{-\zeta}\,\sum_{n=0}^{\infty}\alpha_{s}(p^{2a}\zeta^{-a})^{n+1}\tilde{c}^{\prime}_{n}(0)$
$\displaystyle=$ $\displaystyle
Q^{2u^{\prime}}\int_{0}^{\infty}d\zeta\,\zeta^{a(\bar{u}-u^{\prime})-1}e^{-\zeta}\int_{C_{\mp}}d\bar{p}^{2}\,e^{-\zeta
t\bar{p}^{2}}\frac{1}{(\bar{p}^{2})^{1+au^{\prime}}}\,\sum_{n=0}^{\infty}\alpha_{s}(\bar{p}^{2a})^{n+1}\tilde{c}^{\prime}_{n}(0)$
$\displaystyle=$ $\displaystyle
Q^{2u^{\prime}}\int_{C_{\mp}}\frac{d\bar{p}^{2}}{(\bar{p}^{2})^{1+au^{\prime}}}\,\frac{\Gamma(a(\bar{u}-u^{\prime}))}{(1+t\bar{p}^{2})^{a(\bar{u}-u^{\prime})}}\sum_{n=0}^{\infty}\alpha_{s}(\bar{p}^{2a})^{n+1}\tilde{c}^{\prime}_{n}(0).$

Recall that $t=Q^{-2/a}$. In the third line we have changed the variable of
integration from $p^{2}$ to $\bar{p}^{2}=p^{2}\zeta^{-1}$.

##  Appendix 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. Jantzen:2011nz
.

Let us take the observable $I(t)$ studied in Sec. 4.2 as an example:

$\displaystyle I(t)$ $\displaystyle=$ $\displaystyle
t^{2}\int_{0}^{\infty}dp^{2}\,\tilde{I}(p^{2})\,e^{-tp^{2}}$ (116)
$\displaystyle=$ $\displaystyle\lim_{\delta\to
0}~{}t^{2-\delta}\int_{0}^{\infty}dp^{2}\,p^{-2\delta_{\pm}}\,\tilde{I}(p^{2})_{\pm}\,e^{-tp^{2}}\,.$

For convenience of notation, instead of deforming the contour to $C_{\pm}$, we
have taken the integral path on the positive real axis and instead shifted the
Landau pole infinitesimally above or below the axis, $m^{2}_{\pm}\equiv
m^{2}\pm i0$, and write
$\tilde{I}(p^{2})_{\pm}\equiv\tilde{I}(p^{2})|_{m^{2}\to m^{2}_{\pm}}$. We
note that $\tilde{I}(p^{2})$ is regular at ${\rm Re}\,p^{2}>0$.

We introduce the factorization scale $\mu_{F}$ such that $\Lambda_{\rm
QCD}^{2}\ll\mu_{F}^{2}\ll t^{-1}$ and separate the integral region to above
and below $\mu_{F}^{2}$:

$\displaystyle I_{\rm reg}$ $\displaystyle\equiv$ $\displaystyle
t^{2-\delta}\int_{0}^{\infty}dp^{2}\,p^{-2\delta_{\pm}}\,\tilde{I}(p^{2})_{\pm}\,e^{-tp^{2}}$
$\displaystyle=$ $\displaystyle
t^{2-\delta}\bigg{(}\int_{0}^{\mu_{F}^{2}}+\int_{\mu_{F}^{2}}^{\infty}\bigg{)}dp^{2}\,p^{-2\delta_{\pm}}\tilde{I}(p^{2})_{\pm}\,e^{-tp^{2}}$
$\displaystyle=$ $\displaystyle
t^{2-\delta}\int_{0}^{\mu_{F}^{2}}dp^{2}p^{-2\delta_{\pm}}\tilde{I}(p^{2})_{\pm}T_{tp^{2}}\big{[}e^{-tp^{2}}\big{]}+t^{2-\delta}\int_{\mu_{F}^{2}}^{\infty}dp^{2}p^{-2\delta_{\pm}}\big{[}\tilde{I}(p^{2})_{\pm}\big{]}_{\rm
OPE}e^{-tp^{2}}\,,$

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. ($[\tilde{I}(p^{2})_{\pm}]_{\rm OPE}$ represents double
expansion in $1/p^{2}$ and $1/\log(p^{2}/m_{\pm}^{2})$.)

We further rewrite the integral regions:

$\displaystyle I_{\rm reg}$ $\displaystyle=$ $\displaystyle
t^{2-\delta}\bigg{(}\int_{0}^{\infty}-\int_{\mu_{F}^{2}}^{\infty}\bigg{)}dp^{2}\,p^{-2\delta_{\pm}}\,\tilde{I}(p^{2})_{\pm}\,T_{tp^{2}}\big{[}e^{-tp^{2}}\big{]}$
(118)
$\displaystyle\quad\quad+t^{2-\delta}\bigg{(}\int_{0}^{\infty}-\int_{0}^{\mu_{F}^{2}}\bigg{)}dp^{2}\,p^{-2\delta_{\pm}}\,\big{[}\tilde{I}(p^{2})_{\pm}\big{]}_{\rm
OPE}\,e^{-tp^{2}}$ $\displaystyle=$ $\displaystyle
t^{2-\delta}\int_{0}^{\infty}dp^{2}p^{-2\delta_{\pm}}\tilde{I}(p^{2})_{\pm}T_{tp^{2}}\big{[}e^{-tp^{2}}\big{]}+t^{2-\delta}\int_{0}^{\infty}dp^{2}p^{-2\delta_{\pm}}\big{[}\tilde{I}(p^{2})_{\pm}\big{]}_{\rm
OPE}e^{-tp^{2}}$
$\displaystyle\quad\quad-t^{2-\delta}\bigg{(}\int_{0}^{\mu_{F}^{2}}+\int_{\mu_{F}^{2}}^{\infty}\bigg{)}dp^{2}p^{-2\delta_{\pm}}\big{[}\tilde{I}(p^{2})_{\pm}\big{]}_{\rm
OPE}T_{tp^{2}}\big{[}e^{-tp^{2}}\big{]}.$

In the last expression, the first two terms are equal to $I_{s}^{\pm}(t)$ and
$I_{h}^{\pm}(t)$, respectively. The last term is the sum of terms in the form
of Eq. (52) or (105), which are zero by analytic continuation. This completes
the proof.

##  Appendix D Calculations for OPEs in nonlinear $\sigma$ model

In this appendix we present details of the calculation for the OPEs of the
observables in the nonlinear $\sigma$ 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 $\sigma$ model.

Consider an observable given as an integral form

$S(q^{2})=\int\prod_{i}dx_{i}\int d^{D}p\,\,W(x_{i};q,p)\,G(p^{2})\,,$ (119)

where $q,p$ are vectors in $D$ dimension. $G(p^{2})$ is given as a
perturbative series

$G(p^{2})=(p^{2})^{-\nu}\sum_{k}c_{k}(L_{p})\,\bar{g}(\mu^{2})^{k+1}\,,~{}~{}~{}~{}~{}L_{p}=\log\left(\frac{\mu^{2}}{p^{2}}\right)\,.$
(120)

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
$p^{2}$:

$F(a)\equiv\int\prod_{i}dx_{i}\int
d^{D}p\,\,W(x_{i};q,p)\,\left(\frac{\mu^{2}}{p^{2}}\right)^{a}\,.$ (121)

Then, the Borel transform of $S(q^{2})$ is given in a factorized form as

$B_{S}(u)=F(\nu+u)\,B_{G}(u;p^{2}=\mu^{2})\,.$ (122)

This formula is useful in the case that one can calculate $B_{G}(u)$ and
$F(a)$, while direct evaluation of $B_{S}(u)$ via Eq. (119) is difficult.
(Proof)
By setting $\mu^{2}=p^{2}$ and using the running formula Eq. (114), we can
write

$\displaystyle
G(p^{2})=(p^{2})^{-\nu}\sum_{k}c_{k}(0)\,\bar{g}(p^{2})^{k+1}=\left(\frac{\mu^{2}}{p^{2}}\right)^{\hat{H}+\nu}\sum_{k}c_{k}(0)\,\bar{g}(\mu^{2})^{k+1}(\mu^{2})^{-\nu}\,,$
(123)
$\displaystyle\hat{H}=\bar{g}^{2}\frac{\partial}{\partial\bar{g}}\,,~{}~{}~{}~{}~{}\text{in
particular,}~{}~{}~{}\hat{H}\,\bar{g}^{k+1}=(k+1)\,\bar{g}^{k+2}\,.$ (124)

Inserting this into Eq. (119) and using Eq. (121), we obtain

$\displaystyle S(q^{2})$
$\displaystyle=F(\hat{H}+\nu)\sum_{k}c_{k}(0)\,\bar{g}(\mu^{2})^{k+1}\,(\mu^{2})^{-\nu}$
$\displaystyle=\sum_{j}\frac{F^{(j)}(\nu)}{j!}\,\hat{H}^{j}\sum_{k}c_{k}(0)\,\bar{g}(\mu^{2})^{k+1}\,(\mu^{2})^{-\nu}$
$\displaystyle=\sum_{j}\frac{F^{(j)}(\nu)}{j!}\,\sum_{k}c_{k}(0)\,\frac{(k+j)!}{k!}\,\bar{g}(\mu^{2})^{k+j+1}\,(\mu^{2})^{-\nu}\,.$
(125)

We take the Borel transform and obtain

$\displaystyle B_{S}(u)$
$\displaystyle=\sum_{j}\frac{F^{(j)}(\nu)}{j!}\,u^{j}\sum_{k}\,\frac{c_{k}(0)}{k!}\,u^{k}\,(\mu^{2})^{-\nu}$
$\displaystyle=F(\nu+u)\,B_{G}(u;p^{2}=\mu^{2})\,.$ (126)

###  D.2 Calculations of $B_{\tilde{I}}$, $X_{n}$, $B_{I_{h}}$

We can calculate the Borel transforms $B_{\tilde{I}}$ and $X_{n}$ as follows.
For the moment we fix $\mu^{2}=p^{2}$. Define

$\displaystyle f(x,\bar{g})$ $\displaystyle\equiv$
$\displaystyle\frac{\tilde{I}(p^{2})}{p^{2}}=\frac{\sqrt{1+4x}}{\log\Big{(}\frac{\sqrt{1+4x}+1}{\sqrt{1+4x}-1}\Big{)}}=\sum_{n=0}^{\infty}x^{n}\sum_{k=0}^{\infty}d_{n,k}(0)\,\bar{g}(p^{2})^{k+1}\,,$
(127)

where

$\displaystyle
x=\frac{m^{2}}{p^{2}},~{}~{}~{}\,\bar{g}(p^{2})=\frac{1}{\log(p^{2}/m^{2})}\,.$
(128)

The Taylor expansion of $f(x,\bar{g})$ in $\bar{g}(p^{2})$ is given by

$\displaystyle f(x,\bar{g})$ $\displaystyle=$
$\displaystyle\frac{\sqrt{1+4x}}{\log\Big{(}\frac{(\sqrt{1+4x}+1)^{2}}{4x}\Big{)}}=\frac{\bar{g}\sqrt{1+4x}}{1+2\bar{g}\log\Big{(}\frac{\sqrt{1+4x}+1}{2}\Big{)}}$
(129) $\displaystyle=$
$\displaystyle\sqrt{1+4x}\sum_{k=0}^{\infty}\bigg{[}-2\log\Big{(}\frac{\sqrt{1+4x}+1}{2}\Big{)}\bigg{]}^{k}\,\bar{g}^{k+1}.$

The Borel transform of $f$ is readily obtained as

$B_{f}(u,x)=\sqrt{1+4x}\sum_{k=0}^{\infty}\frac{1}{k!}\bigg{[}-2\log\Big{(}\frac{\sqrt{1+4x}+1}{2}\Big{)}\bigg{]}^{k}\,u^{k}=\sqrt{1+4x}\bigg{(}\frac{\sqrt{1+4x}+1}{2}\bigg{)}^{-2u}\,,$
(130)

and $B_{\tilde{I}}=p^{2}B_{f}$.

$X_{n}$ is included in $B_{f}$ as
$B_{f}(u,x)=\sum_{n=0}^{\infty}X_{n}(u)\,x^{n}$. Taylor expansion of $B_{f}$
in $x$ is calculated as follows:

$\displaystyle B_{f}(u,x)$ $\displaystyle=$
$\displaystyle\sum_{m=0}^{\infty}\frac{4^{u}\Gamma(1-2u)}{\Gamma(m+1)\Gamma(1-m-2u)}(1+4x)^{\frac{m+1}{2}}$
(131) $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\frac{4^{u+n}\Gamma(1-2u)\Gamma(3/2+m/2)}{\Gamma(m+1)\Gamma(1-m-2u)\Gamma(3/2+m/2-n)\Gamma(n+1)}x^{n}$
$\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\bigg{[}4^{n+u}\frac{\Gamma(3/2)\,_{3}F_{2}\big{(}\\{3/2,\,u,\,1/2+u\\},\\{1/2,\,3/2-n\\},1\big{)}}{\Gamma(3/2-n)\Gamma(1+n)}$
$\displaystyle\quad\quad+4^{n+u}\frac{\Gamma(1-2u)\,_{3}F_{2}\big{(}\\{1+n,\,-1/2+n+u,\,n+u\\},\\{n,\,1/2+n\\},1\big{)}}{\Gamma(2n)\Gamma(2-2n-2u)}\bigg{]}x^{n}$
$\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\bigg{[}\frac{\Gamma(-n-2u)}{\Gamma(n+1)\Gamma(2-2n-2u)}(4u^{2}-2u-2n)\bigg{]}x^{n},$

where we used an identity

${}_{3}F_{2}(\\{a+1,\,b,\,c\\},\\{a,\,d\\},1)=\frac{\Gamma(d)\Gamma(d-b-c-1)}{\Gamma(d-b)\Gamma(d-c)}\Big{(}d-b-c-1+\frac{bc}{a}\Big{)}.$
(132)

Let us now calculate the Borel transform $B_{f}(u,x)$ with respect to
$\bar{g}(\mu^{2})$, namely in the case $\mu$ is left as an independent
parameter. In this case, using the running formula Eq. (114), we have

$\displaystyle
f(x,\bar{g})=\sum_{n=0}^{\infty}x^{n}\sum_{k=0}^{\infty}d_{n,k}(0)\,\bar{g}(p^{2})^{k+1}=\left(\frac{\mu^{2}}{p^{2}}\right)^{\hat{H}}\sum_{n=0}^{\infty}x^{n}\sum_{k=0}^{\infty}d_{n,k}(0)\,\bar{g}(\mu^{2})^{k+1}\,.$
(133)

This has the same form as the first line of Eq. (125) with $\nu=0$ and
$F(a)=(\mu^{2}/p^{2})^{a}$. Hence,

$B_{f}(u,x;p^{2};\mu^{2})=\biggl{(}\frac{\mu^{2}}{p^{2}}\biggr{)}^{\\!u}\,B_{f}(u,x;\mu^{2};\mu^{2})=\sqrt{1+4x}\bigg{(}\frac{\sqrt{1+4x}+1}{2}\bigg{)}^{\\!-2u}\biggl{(}\frac{\mu^{2}}{p^{2}}\biggr{)}^{\\!u}\,.$
(134)

To calculate $B_{I_{h}}(u)$, we can apply the factorization formula Eq. (122)
in the case where $\nu=n-1$,
$G(p^{2})=p^{2}(m^{2}/p^{2})^{n}\tilde{C}_{n}(p^{2})$ [cf., Eqs. (61), (62)],
and

$F(a)=t^{2-\delta}\\!\int_{C_{\mp}}dp^{2}\,(p^{2})^{-\delta_{\pm}}\,e^{-tp^{2}}\left(\frac{\mu^{2}}{p^{2}}\right)^{a}=(\mu^{2}t)^{a}\,t\,\Gamma(1-\delta-a)\,.$
(135)

Then, we readily obtain Eq. (68).

###  D.3 Calculations of $B_{\rho_{h}}$, $B_{\tilde{\rho}}$

There are three regions when the EBR method is applied to $\rho(q^{2})$
defined by Eq. (93): (i) $k\sim q\gg m$, (ii) $(q+k)\sim m\ll q,k$, and (iii)
$k\sim m\ll q$. In order to regularize each integral in the expansion, we
employ dimensional regularization except for $D_{\alpha}(k^{2})$ whose
expression for $D=2$, Eq. (31), is used:151515 This treatment is for the sake
of simplicity, since keeping $D_{\alpha}(k^{2})$ also in dimensional
regularization makes the calculations tedious. In the end, the dual space
quantities $\tilde{\rho}$ and $B_{\tilde{\rho}}$ are obtained to be finite as
$D\to 2$, so this treatment should be justified.

$\rho(q^{2})_{\rm reg}\equiv q^{2}\left(\frac{\partial}{\partial
q^{2}}\right)^{2}\int\frac{d^{D}k}{(2\pi)^{D}}\frac{1}{(q+k)^{2}+m^{2}}\,N\left[D_{\alpha}(k^{2})\right]_{D=2}\,.$
(136)

To compute the dual space observable $[\tilde{\rho}(p^{2})]_{\rm 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.

$\displaystyle\rho(q^{2})_{\rm reg}^{\rm(i,ii)}$ $\displaystyle=$
$\displaystyle q^{2}\left(\frac{\partial}{\partial
q^{2}}\right)^{2}\int\frac{d^{D}k}{(2\pi)^{D}}\frac{1}{(q+k)^{2}+m^{2}}\,N\left[D_{\alpha}(k^{2})\right]_{\rm
OPE}$ (137) $\displaystyle=$ $\displaystyle
4\pi\mu^{2}\sum_{n,k}F(n+\hat{H}-1)\left(\frac{m^{2}}{\mu^{2}}\right)^{n}d_{n,k}(0)\,\bar{g}(\mu^{2})^{k+1}\,,$

where

$\displaystyle F(a)$ $\displaystyle=q^{2}\left(\frac{\partial}{\partial
q^{2}}\right)^{2}\int\frac{d^{D}k}{(2\pi)^{D}}\frac{1}{(q+k)^{2}+m^{2}}\,\left(\frac{\mu^{2}}{k^{2}}\right)^{a}$
$\displaystyle=\frac{m^{-2\delta}}{(4\pi)^{1-\delta}}\left(\frac{\mu^{2}}{m^{2}}\right)^{\\!a}\,\frac{\Gamma(1-a-\delta)\Gamma(a+\delta)}{\Gamma(1-\delta)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times
q^{2}\left(\frac{\partial}{\partial
q^{2}}\right)^{2}\,_{2}F_{1}\left(a,a+\delta;1-\delta;-\frac{q^{2}}{m^{2}}\right)\,.$
(138)

We set $D=2-2\delta$.

Hence, using the factorization formula, we obtain the Borel transform of
$\rho_{h}^{\pm}(q^{2})$ as

$B_{\rho_{h}}(u)=4\pi\mu^{2}\sum_{n=0}^{\infty}F(n+u-1)\,X_{n}(u)\,\left(\frac{m^{2}}{\mu^{2}}\right)^{n}\,.$
(139)

We can rewrite ${}_{2}F_{1}(a,b;c;-\frac{q^{2}}{m^{2}})$ in terms of
${}_{2}F_{1}(a^{\prime},b^{\prime};c^{\prime};-\frac{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 $\rho_{h}^{\pm}(q^{2})$ originate from the factor
$\Gamma(2-n-u-\delta)$ in $F(n+u-1)$. We can take the limit $\delta\to 0$ for
the other part, which gives the formula Eq. (101).

We can take the dual transform of Eq. (101) using Eq. (96), order by order in
expansion in $1/q^{2}$. The dual transform essentially generates
$1/\Gamma(2\\!-\\!M\\!-\\!u\\!-\\!\delta)$. Since
$\Gamma(2\\!-\\!n\\!-\\!u\\!-\\!\delta)/\Gamma(2\\!-\\!M\\!-\\!u\\!-\\!\delta)$
for $n\leq 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 $\delta\to 0$. Thus, we
obtain $B_{\tilde{\rho}}$ given by Eq. (98). It is straightforward to obtain
$[\tilde{\rho}(p^{2})]_{\rm OPE}$ from $B_{\tilde{\rho}}$ at each order of the
expansion in $1/p^{2}$ and $\bar{g}(p^{2})$. We expect that the obtained
result for $B_{\tilde{\rho}}(u)$ or $[\tilde{\rho}(p^{2})]_{\rm 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.161616 The algorithm for deriving $[\tilde{\rho}(p^{2})]_{\rm
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 $\rho_{h}^{\pm}(q^{2})$, and afterwards we
remove the regularization. Once $[\tilde{\rho}(p^{2})]_{\rm OPE}$ is
determined, $\rho_{h}^{\pm}(q^{2})$ can be redetermined regularization
dependently. The regularization dependence cancels against that of
$\rho_{s}^{\pm}(q^{2})$.

## References

  * (1) D. J. Gross and A. Neveu, Phys. Rev. D 10 (1974), 3235.
  * (2) B. E. Lautrup, Phys. Lett. B 69 (1977), 109-111
  * (3) G. ‘t Hooft, lectures given at the ”Ettore Majorana” Int. School of Subnuclear Physics, Erice, July 1977.
  * (4) M. Beneke, Phys. Rept. 317, 1 (1999). [hep-ph/9807443].
  * (5) U. Aglietti and Z. Ligeti, Phys. Lett. B 364 (1995), 75 [arXiv:hep-ph/9503209 [hep-ph]].
  * (6) M. Beneke and V. M. Braun, Nucl. Phys. B 426 (1994), 301-343 [arXiv:hep-ph/9402364 [hep-ph]].
  * (7) I. I. Y. Bigi, M. A. Shifman, N. G. Uraltsev and A. I. Vainshtein, Phys. Rev. D 50, 2234-2246 (1994) [arXiv:hep-ph/9402360 [hep-ph]].
  * (8) A. Pineda, Ph.D. Thesis (1998).
  * (9) A. H. Hoang, M. C. Smith, T. Stelzer and S. Willenbrock, Phys. Rev. D 59, 114014 (1999) [arXiv:hep-ph/9804227 [hep-ph]].
  * (10) M. Beneke, Phys. Lett. B 434, 115-125 (1998) [arXiv:hep-ph/9804241 [hep-ph]].
  * (11) M. Neubert and C. T. Sachrajda, Nucl. Phys. B 438 (1995), 235-260 [arXiv:hep-ph/9407394 [hep-ph]].
  * (12) P. Ball, M. Beneke and V. M. Braun, Phys. Rev. D 52 (1995), 3929-3948 [arXiv:hep-ph/9503492 [hep-ph]].
  * (13) A. A. Penin and N. Zerf, JHEP 04 (2014), 120 [arXiv:1401.7035 [hep-ph]].
  * (14) Y. Kiyo, G. Mishima and Y. Sumino, Phys. Lett. B 752 (2016) 122 Erratum: [Phys. Lett. B 772 (2017) 878] [arXiv:1510.07072 [hep-ph]].
  * (15) M. Beneke, A. Maier, J. Piclum and T. Rauh, PoS RADCOR2015 (2016), 035 [arXiv:1601.02949 [hep-ph]].
  * (16) A. Bazavov et al. [Fermilab Lattice, MILC and TUMQCD], Phys. Rev. D 98 (2018) no.5, 054517 [arXiv:1802.04248 [hep-lat]].
  * (17) C. Peset, A. Pineda and J. Segovia, JHEP 09 (2018), 167 [arXiv:1806.05197 [hep-ph]].
  * (18) A. H. Hoang, Z. Ligeti and A. V. Manohar, Phys. Rev. D 59 (1999), 074017 [arXiv:hep-ph/9811239 [hep-ph]].
  * (19) A. H. Hoang, Z. Ligeti and A. V. Manohar, Phys. Rev. Lett. 82 (1999), 277-280 [arXiv:hep-ph/9809423 [hep-ph]].
  * (20) A. Alberti, P. Gambino, K. J. Healey and S. Nandi, Phys. Rev. Lett. 114 (2015) no. 6, 061802 [arXiv:1411.6560 [hep-ph]].
  * (21) A. Bazavov, N. Brambilla, X. Garcia Tormo, i, P. Petreczky, J. Soto and A. Vairo, Phys. Rev. D 86 (2012), 114031 [arXiv:1205.6155 [hep-ph]].
  * (22) K. G. Wilson, Phys. Rev. 179 (1969) 1499.
  * (23) A. H. Mueller, Nucl. Phys. B 250 (1985), 327-350.
  * (24) H. Takaura, T. Kaneko, Y. Kiyo and Y. Sumino, Phys. Lett. B 789, 598-602 (2019) [arXiv:1808.01632 [hep-ph]].
  * (25) H. Takaura, T. Kaneko, Y. Kiyo and Y. Sumino, JHEP 04, 155 (2019) [arXiv:1808.01643 [hep-ph]].
  * (26) C. Ayala, X. Lobregat and A. Pineda, [arXiv:2005.12301 [hep-ph]].
  * (27) C. Ayala, X. Lobregat and A. Pineda, JHEP 12 (2020), 093 [arXiv:2009.01285 [hep-ph]].
  * (28) C. Ayala, X. Lobregat and A. Pineda, Phys. Rev. D 99 (2019) no.7, 074019 [arXiv:1902.07736 [hep-th]].
  * (29) C. Ayala, X. Lobregat and A. Pineda, Phys. Rev. D 101 (2020) no.3, 034002 [arXiv:1909.01370 [hep-ph]].
  * (30) T. Lee, Phys. Rev. D 67, 014020 (2003) [arXiv:hep-ph/0210032 [hep-ph]].
  * (31) H. Takaura, JHEP 10, 039 (2020) [arXiv:2002.00428 [hep-ph]].
  * (32) Y. Hayashi, Y. Sumino and H. Takaura, Phys. Lett. B 819 (2021), 136414 [arXiv:2012.15670 [hep-ph]].
  * (33) Y. Hayashi, Y. Sumino and H. Takaura, JHEP 02 (2022), 016 [arXiv:2106.03687 [hep-ph]].
  * (34) Y. Hayashi, JHEP 06 (2022), 157 [arXiv:2112.14408 [hep-ph]].
  * (35) M. A. Benitez-Rathgeb, D. Boito, A. H. Hoang and M. Jamin, JHEP 07, 016 (2022) [arXiv:2202.10957 [hep-ph]].
  * (36) M. A. Benitez-Rathgeb, D. Boito, A. H. Hoang and M. Jamin, JHEP 09, 223 (2022) [arXiv:2207.01116 [hep-ph]].
  * (37) M. Beneke and V. A. Smirnov, Nucl. Phys. B 522 (1998), 321-344 [arXiv:hep-ph/9711391 [hep-ph]].