text:
# Weakly nonlinear theory of grain boundary motion in patterns with crystalline symmetry

## Abstract

We study the motion of a grain boundary separating two otherwise stationary domains of hexagonal symmetry. Starting from an order parameter equation appropriate for hexagonal patterns, a multiple scale analysis leads to an analytical equation of motion for the boundary that shares many properties with that of a crystalline solid. We find that defect motion is generically opposed by a pinning force that arises from non-adiabatic corrections to the standard amplitude equation. The magnitude of this force depends sharply on the mis-orientation angle between adjacent domains so that the most easily pinned grain boundaries are those with an angle $`4^{}\theta 8^{}`$. Although pinning effects may be small, they do not vanish asymptotically near the onset of this subcritical bifurcation, and can be orders of magnitude larger than those present in smectic phases that bifurcate supercritically.

The microstructure of a condensed phase and the distribution of topological defects largely determine its mechanical and thermodynamical response, as well as the temporal evolution of its non-equilibrium configurations. Continuum or hydrodynamic approaches to phases with broken symmetry are now well understood , including the long wavelength description of topological defects . Nevertheless, a quantitative theory of defect motion (e.g., dislocation glide) remains a difficult problem because it requires short scale phenomena that lie beyond a hydrodynamic theory, and therefore beyond the degree of universality that such a description entails.

Bifurcations to states with broken symmetry are also encountered in a variety of physical, chemical and biological systems driven outside of thermodynamic equilibrium. The characteristic length scales involved are much larger than atomic dimensions, therefore making the observation and study of defects substantially easier. Here, amplitude or phase equations that focus on slow modulations of a base periodic pattern play the role of the hydrodynamic description . It has been noted, for example, that at leading order hexagonal patterns are the dissipative analogues of a two dimensional, isotropic solid . We focus in this paper on grain boundary motion in hexagonal patterns, and find many qualitative and quantitative similarities with grain boundary motion in crystalline solids. In contrast with the latter case, our starting model allows a detailed analysis of the breakdown of the long wavelength description of defect motion through an explicit multiple scale analysis. We are then able to derive analytically several results that are known only qualitatively or empirically in crystalline solids.

Hexagonal patterns formed in a spatially extended system (like in Langmuir monolayers  or block copolymer melts ) are often fairly disordered. They consist of randomly oriented grains separated by grain boundaries (arrays of dislocations in the low angle case), and are equivalent by symmetry to a polycrystalline solid. This nonequilibrium microstructure usually evolves on a very slow, fluctuation-dependent time scale due to pinning forces to defect motion. Our focus here is on a coarse grained model of a hexagonal pattern, and, in particular, on the motion of a grain boundary separating two domains with arbitrary mis-orientation. We extend to this case a recent study of defect motion and pinning near a supercritical bifurcation involving stripe patterns (or smectic phases) . We show that the coupling between the slow variables of an envelope description of the defect and the underlying (fast) periodicity of the base pattern leads to an effective periodic potential, the close analogue of the Peierls barrier acting on defects in a crystalline solid . Furthermore, we find that the magnitude of the potential barrier does not vanish as the subcritical bifurcation point is approached, contrary to the case of smectic phases . Hence, not only crystalline phases emerging from a subcritical bifurcation are harder than smectics in terms of the forces acting on topological defects, but also an envelope description that ignores the internal degrees of freedom of defects does not appear to be valid even near the bifurcation point. The self-generated pinning effects discussed here are expected to be a general feature of modulated phases, and can explain, for example, many grain boundary conformations that have been observed experimentally and numerically in block copolymer melts .

We study here the Swift-Hohenberg model of Rayleigh-Bénard convection with an additional quadratic term to allow the formation of hexagonal patterns ,

$$\frac{\psi }{t}=ϵ\psi \frac{1}{k_0^4}(k_0^2+^2)^2\psi +g_2\psi ^2\psi ^3.$$
(1)
The order parameter $`\psi (\stackrel{}{x},t)`$ is related to the vertical velocity at the mid plane of the convective cell, $`ϵ`$ is the reduced Rayleigh number, and $`g_2`$ can be related to deviations from Boussinesq behavior in the working fluid. The uniform solution $`\psi =0`$ becomes unstable for $`ϵ>0`$ to a periodic pattern of wavelength $`\lambda _0=2\pi /k_0`$. Hexagonal patterns are stable for $`|ϵ_m(g_2)|<ϵ<ϵ_M(g_2)`$, and roll patterns for $`ϵ>ϵ_M`$. In this study, we chose $`ϵ[0,ϵ_M]`$ so that only hexagonal patterns are stable.

An approximate stationary solution for a configuration containing a planar grain boundary between two uniform and symmetric domains of hexagons that have a relative mis-orientation angle $`\theta `$ ($`0\theta \pi /3`$, see Fig. 1) can be found by assuming that $`\psi _0=_{n=1}^6A_ne^{i\stackrel{}{k}_n\stackrel{}{x}}+c.c.`$, where $`A_n(x)`$ are slowly varying envelopes, and $`x`$ is the coordinate normal to the boundary. As $`x\mathrm{}`$, $`A_{1,2,3}A_0=(g_2+\sqrt{g_2^2+15ϵ})/15`$ exponentially fast outside a boundary layer of thickness $`\xi `$ \[$`A_n(x)=f_n(x/\xi )`$\], whereas $`A_{4,5,6}0`$. As $`x\mathrm{}`$, $`A_{4,5,6}A_0`$ and $`A_{1,2,3}0`$. We chose $`|\stackrel{}{k}_n|=k_0`$, $`\stackrel{}{k}_{1,4}=k_0(\mathrm{cos}(\varphi )\widehat{x}\mathrm{sin}(\varphi )\widehat{y})`$, where $`\varphi =\pi /6\theta /2`$, and the other vectors $`\stackrel{}{k}_n`$ are obtained from these by rotations of $`\pm 2\pi /3`$ as indicated in Figure 1. In analogy with crystalline solids, a small angle grain boundary is well described as an array of penta-hepta defects (the dislocations cores of a hexagonal pattern), separated on average by a distance of the order $`\lambda _0/\theta `$. However, since the projections of $`\stackrel{}{k}_{1,2,3}`$ on the $`x`$ axis are usually not commensurate, the patterns observed along the boundary are not periodic.

We next focus on grain boundary motion and extend our earlier results for stripe patterns . The amplitude equation formalism has already been used to study defect dynamics in hexagonal patterns . Amplitude equations are obtained from a standard multiscale analysis of Eq. (1). For example, the equation for $`A_1`$ follows from the solvability condition

$$_x^{x+l_x}\frac{dx^{}}{l_x}\underset{l_y\mathrm{}}{lim}_0^{l_y}\frac{dy^{}}{l_y}[L(\psi _0)+g_2\psi _0^2\psi _0^3]e^{i\stackrel{}{k}_1\stackrel{}{x}^{}}=0,$$
(2)
where $`L`$ is a linear operator that follows from Eq. (1)  and $`l_x`$ a length of $`O(\lambda _0)`$ to be specified later. When both $`ϵ`$ and $`g_20`$, the length scale of variation of the $`A_n`$’s (i.e. the grain boundary thickness $`\xi `$) is much larger than $`\lambda _0`$, and only non oscillatory terms contribute to the integral (2). The solvability condition leads to the usual amplitude equations that show that defects are either immobile or move with constant velocity .

If, on the other hand, $`ϵ`$ and $`g_2`$ are small but finite, then the amplitudes are not strictly constant within a period $`l_x`$. Any oscillatory term in (2) of wavevector $`\stackrel{}{K}`$ parallel to the normal $`x`$-axis does not integrate to zero, and the equation for the slowly varying amplitudes cannot be decoupled from the phase of the defect (yielding “non-adiabatic corrections” ). Although these contributions are small (non-analytic in both $`ϵ`$ and $`g_2`$), they may have dramatic effects. In the case of Eq. (1), non-vanishing terms can arise from the cubic nonlinearity, and are proportional to $`A_1^2A_4`$, $`\overline{A}_1\overline{A}_4^2`$, $`A_1A_3A_5`$, $`A_1\overline{A}_3\overline{A}_5`$, $`A_1A_2A_6`$ and $`A_1\overline{A}_2\overline{A}_6`$. An oscillatory term proportional to $`\mathrm{cos}(Kx^{})`$ in Eq. (2) contributes to order $`\mathrm{exp}(|K\xi |)`$ to the law of grain boundary of motion (see Eqs. (9) and (12) below). This contribution is simplest when both $`ϵ1`$ and $`g_21`$ so that $`\xi \lambda _0`$ and a single mode with lowest $`K`$ dominates. In this limit, the slowest non-adiabatic corrections are those proportional to $`A_1A_2A_6`$ and $`A_1\overline{A}_2\overline{A}_6`$, which have $`K=|\stackrel{}{k}_2+\stackrel{}{k}_6|=2k_0\mathrm{sin}(\theta /2)`$ (see Fig. 1). The solvability condition (2) now reads

$`{\displaystyle \frac{A_1}{t}}`$ $`=`$ $`{\displaystyle \frac{F}{\overline{A}_1}}`$ (3)
$``$ $`{\displaystyle _{x_0}^{x_0+l_x}}{\displaystyle \frac{dx}{l_x}}6A_1(A_2A_6e^{2ik_0x\mathrm{sin}(\theta /2)}+\mathrm{c}.\mathrm{c}.),`$ (4)
where $`l_x=\lambda _0/[2\mathrm{sin}(\theta /2)]`$. The first two terms in Eq. (4) represent the standard amplitude equation with the Lyapunov functional

$`F`$ $`=`$ $`{\displaystyle }d\stackrel{}{r}[ϵ{\displaystyle \underset{n=1}{\overset{6}{}}}|A_n|^2+{\displaystyle \frac{4}{k_0^2}}{\displaystyle \underset{n=1}{\overset{6}{}}}|D_nA_n|^2`$ (5)
$``$ $`2g_2(\overline{A}_1\overline{A}_2\overline{A}_3+\overline{A}_4\overline{A}_5\overline{A}_6+A_1A_2A_3+A_4A_5A_6)`$ (6)
$`+`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \underset{n=1}{\overset{6}{}}}|A_n|^4+3{\displaystyle \underset{n<m}{}}|A_n|^2|A_m|^2],`$ (7)
with $`D_n=/x_n`$ the derivative along $`\stackrel{}{k}_n`$. Equations similar to (4) can be derived for the remaining amplitudes $`A_n`$ (not shown). The last term in the r.h.s. of Eq. (4) is new and represents the dominant non-adiabatic correction in the limits $`ϵ1`$ and $`g_21`$.

In order to derive an equation of motion for the grain boundary from the equations for the $`A_n`$, we first denote by $`a_n(x)`$ ($`1n6`$) the leading order amplitudes of the stationary grain boundary, solutions of the one-dimensional coupled Ginzburg-Landau equations $`F/\overline{a}_n=0`$. We then look for solutions of the form $`A_n(x,t)=a_n(xx_{gb}(t))`$, where $`x_{gb}(t)`$ is the now time-dependent position of the grain boundary . We find,

$$D\dot{x}_{gb}=p_{\mathrm{hex}}\mathrm{sin}[2k_0x_{gb}\mathrm{sin}(\theta /2)],$$
(8)
where $`D=_{\mathrm{}}^{\mathrm{}}𝑑x_{n=1}^6(_xa_n)^2`$ is a friction coefficient, and

$`p_{\mathrm{hex}}`$ $`=`$ $`\mathrm{Max}_\beta {\displaystyle _{\mathrm{}}^{\mathrm{}}}dx\mathrm{cos}[2k_0\mathrm{sin}(2\theta )x+\beta ]\{3[a_2_x(a_2^2a_6)`$ (9)
$`+`$ $`a_6_x(a_6^2a_2)]]+12[a_1_x(a_1a_2a_6)+a_3_x(a_3a_2a_6)`$ (10)
$`+`$ $`a_4_x(a_4a_2a_6)+a_5_x(a_5a_2a_6)]\}`$ (11)
is the (dimensionless) amplitude of a pinning force of wavelength $`\lambda _0/[2\mathrm{sin}(\theta /2)]`$. This periodic force explicitly arises from the coupling between the slowly varying wave amplitudes and the periodicity of the base state in the integral term of Eq. (4). Equations (8) and (9) show that a planar grain boundary initially located at an arbitrary position relaxes toward the nearest minimum of the periodic pinning potential. The stationary and stable positions of the boundary are thus discrete, separated from each other by a distance $`\mathrm{\Delta }x_{gb}=\lambda _0/[2\mathrm{sin}(\theta /2)]>\lambda _0`$. The two wavevectors with the smallest projection on the grain boundary normal ($`\stackrel{}{k}_2`$ and $`\stackrel{}{k}_6`$, see Figure 1) set the wavevector ($`\stackrel{}{k}_2+\stackrel{}{k}_6`$) of the periodic pinning potential. The usual amplitude equation formalism would instead predict that $`p_{\mathrm{hex}}=0`$, and $`x_{gb}`$ is arbitrary and decoupled from the phase of the pattern. Figure 1 shows three successive stable locations of the the grain boundary obtained by numerically solving Eq. (1). The values of $`\mathrm{\Delta }x_{gb}`$ determined numerically agree very well with the analytic result.

We expect equations of the form of Eq. (8) to be generic at a bifurcation, and not limited to hexagon-hexagon grain boundaries. In particular, we anticipate that non-adiabatic effects are important in block copolymer melts, and can explain the conformations of planar interfaces observed in these systems .

Remarkably, the present result for a pattern of hexagonal symmetry is analogous to the Peierls force acting on dislocations in crystalline solids . Peierls calculated the energy of a single dislocation by summing over the interactions between atoms within two neighboring layers, their displacements given by continuous elasticity as a first approximation. The energy of the dislocation oscillates as a function of its position so that it can glide only if a force of finite amplitude acts on it (the Peierls’ force). Here, we find that a similar force acts on assemblies of dislocations organized in arrays. Like in the simpler geometry studied by Peierls, the amplitude of the pinning force decays exponentially with the spatial thickness of the defect $`\xi `$. From Eq. (9), we find

$$p_{\mathrm{hex}}c^{}A_0^4e^{2k_0\mathrm{sin}(\theta /2)\xi a^{}},$$
(12)
with $`c^{}`$ and $`a^{}`$ dimensionless constants of order unity that can a priori depend on the mis-orientation $`\theta `$. At this point, we note an important distinction between defect motion in hexagonal patterns that emerge at a subcritical bifurcation, and in stripe patterns that bifurcate supercritically. The supercritical case was discussed in ref.  for a $`90^{}`$ boundary separating two domains of stripes given by Eq. (1) with $`g_2=0`$. The corresponding pinning force acting on the boundary, $`p_{\mathrm{stripe}}`$, satisfies a relation similar to Eq. (12), with $`\xi 1/\sqrt{ϵ}`$. Hence $`p_{\mathrm{stripe}}\mathrm{exp}(1/\sqrt{ϵ})0`$, as the control parameter $`ϵ0`$. In the hexagonal phase, however, $`\xi (ϵ=0)15\lambda _0/(8\sqrt{6}\pi g_2)`$ is finite. Therefore $`p_{\mathrm{hex}}`$, although small, does not vanishes when $`ϵ0`$ (nor when $`ϵ|ϵ_m(g_2)|`$).

Figure 2 shows typical variations of the pinning force as a function of $`ϵ`$ for different values of $`g_2`$ ($`\theta =\pi /9`$). The curves are given by Eq. (9), where the amplitudes $`a_n`$ have been obtained by numerically integrating the system of equations $`F/\overline{a}_n=0`$. For a value of $`g_2`$ as small as $`0.3`$, $`p_{\mathrm{hex}}`$ is many orders of magnitude larger than $`p_{\mathrm{stripe}}`$. Therefore, non-adiabatic effects in hexagonal, “crystalline-like”, patterns are difficult to avoid. Defects need to overcome much higher activation barriers, and will be much less easily un-pinned either upon the application of external stresses or random noise (that would be represented by additional terms in the r.h.s. of Eq. (8)).

Finally we study the dependence of the pinning force with grain boundary mis-orientation $`\theta `$. Figure 3 displays $`p_{\mathrm{hex}}(\theta )`$ for different values of the parameters $`g_2`$ and $`ϵ`$. We find that the least mobile grain boundaries, i.e. those for which the amplitude of the pinning force is maximal, have a low angle $`\theta _M`$, typically such that $`4^{}\theta _M8^{}`$. Both $`\theta _M`$ and the overall shape of $`p_{\mathrm{hex}}(\theta )`$ depend weakly on $`g_2`$ and $`ϵ`$. Figure 3 also shows the pinning force vs. $`\theta `$ obtained from a direct numerical solution of Eq. (1) with $`g_2=0.3`$ and $`ϵ=0.05`$. Here $`p_{\mathrm{hex}}`$ is estimated by fitting grain boundary relaxation trajectories to Eq. (8), with the friction coefficient $`D`$ assumed to be given by the analytic result. The numerical results compare reasonably well with the theory given the numerical difficulties in tracking the relaxation of the grain boundary in a finite sized system. Although the maximum value is lower and the curve flatter in the numerical case, the right order of magnitude is obtained, as well as a similar range for $`\theta _M`$. The discrepancy can also be attributed in part to non-adiabatic corrections of higher order to Eq. (9), that may not be negligible at $`g_20.3`$.

Our present analysis can qualitatively explain the properties of polycrystalline, partially ordered configurations that are typically observed at long times in many pattern forming systems with this symmetry . Although the evolution defined by Eq. (1) is driven by the minimization of a Lyapunov functional, we expect that grain boundaries (and other topological defects) will become pinned at long times as the driving force for microstructure coarsening decreases. Therefore, the system will eventually reach metastable, glassy-like configurations that can only order by slow activated processes (presumably logarithmic in time, as already observed in  with a random forcing term added to Eq. (1)). We additionally note that defects observed in cold metals are essentially low angle grain boundaries in the range $`5^{}`$ and $`10^{}`$ , a value that compares very well with the values obtained by the present theory.

In summary, we have analyzed the motion of grain boundaries in hexagonal patterns from an order parameter equation and extended the standard Ginzburg-Landau equation for the slowly varying amplitude to incorporate non-analytic corrections. Like in crystalline phases, defect motion is opposed by short range forces with periodicity and amplitude that strongly depend on the mis-orientation angle between domains. Although small, these pinning forces can not be neglected asymptotically at long times in a coarsening system, even near onset, and are orders of magnitude higher than those produced in patterns of smectic symmetry.

We are indebted to François Drolet for useful discussions. This research has been supported by the U.S. Department of Energy, contract No. DE-FG05-95ER14566.


format: text

text:
# Extraction of 𝛼_𝑠 from the Gross–Llewellyn Smith sum rule using Borel resummation

## I Introduction

Many low–energy Quantum Chromodynamic (QCD) observables, including the Gross–Llewellyn Smith (GLS) sum rule, with characteristic energy scale a few GeVs are analyzed in perturbative QCD within the framework of operator product expansion (OPE). In this scheme usually the most important is the perturbative contribution from the Wilson coefficient of the unit operator, and there are nonperturbative, power suppressed, higher–twist contributions. Generally, the higher order coefficients in the perturbative contribution grow rapidly due to the asymptotic nature of the perturbative expansion. The uncalculated higher order corrections are thus expected to be large, and this can cause a large uncertainty in data analysis that employs the unprocessed, finite order perturbative expansion. It is therefore important to properly handle the divergent perturbation expansion via resummation, since it can give a more accurate result with reduced theoretical uncertainty. Besides, resummation serves to give a well defined meaning to the higher–twist contributions. Without a proper resummation of the perturbative part, the higher–twist contributions are ambiguous .

An often used resummation technique is the Borel resummation. It has a sound theoretical basis since it is built on our understanding about the singularities in the Borel plane which cause the divergence of the perturbative expansion. Its use generally improves the quality of data analysis, as can be seen from the reduced dependence on the renormalization scheme and scale, and from the reduced dependence on the uncertainty of the uncalculated next higher order perturbative coefficient.

At moderate values of the strong coupling $`\alpha _s(Q)`$ at a few GeV, the Borel integral receives most of its value from the interval between the origin and the first infrared (IR) renormalon singularity, and just beyond it, in the Borel plane. Let us call this loosely defined interval, for convenience, the primary interval. In Borel resummation it is thus very important to describe the Borel transform, which determines the Borel integral, as accurately as possible in the primary interval using the calculated first terms of perturbation theory.

For this purpose two steps can be taken: (1) an explicit incorporation of the first IR renormalon singularity in the Borel transform, and (2) use of an optimal conformal mapping. With the usual power expansion of the Borel transform about the origin the information on the renormalon singularity is lost. To remedy this, one may explicitly incorporate the first renormalon singularity by writing a Borel transform $`\stackrel{~}{D}(b)`$, which behaves as $`1/(1b/b_0)^{1+\nu }`$ around the singularity at $`b=b_0`$, as

$$\stackrel{~}{D}(b)=\frac{R(b)}{(1b/b_0)^{1+\nu }}$$
(1)
with $`R(b)\stackrel{~}{D}(b)(1b/b_0)^{1+\nu }`$. The function $`R(b)`$***This function was first introduced in  in Borel resummation and independently in  in renormalon residue calculation. is by definition bounded and has softer singularity at the first IR renormalon. Hence we can expect that the Borel transform in the form of (1), with $`R(b)`$ perturbatively expanded about the origin, would give a better approximation than the direct expansion of $`\stackrel{~}{D}(b)`$.

The step (2) utilizes the information on the locations of the singularities in the Borel plane. Use of conformal mapping in Borel resummation has a long history , and its use in perturbative QCD was particularly emphasized in . On the Borel plane there are IR renormalon singularities on the positive real axis and ultraviolet (UV) renormalons on the negative axis. There are also instanton-caused singularities, which can be safely ignored in our case. By pushing the singularities away from the primary interval one can obtain a smoother $`R`$ in the new primary interval on the mapped plane. This would render the perturbation of $`R`$ in the mapped plane to converge better. Among the several mappings considered in the literature we find that the one proposed in  is particularly suited when used in combination with the Borel transform in the form of Eq. (1). This mapping moves the first IR renormalon singularity that defines the primary interval to a point within the unit circle and all other singularities to the unit circle. Since the effect of the first IR renormalon is softened by the step (1), we expect this mapping in our case to be better suited than, for example, a mapping that moves all singularities to the unit circle .

These techniques were applied to the hadronic tau decay width  and to the hadronic contribution to the muon anomalous magnetic moment . In this work we apply them to the GLS sum rule. The CCFR analysis  of the GLS sum rule was based on evaluation of the truncated perturbation series (TPS) in $`\overline{\mathrm{MS}}`$ scheme. Aside from the inherent ambiguity of the higher–twist contributions, this method gives predictions which are not stable under the inclusion of an additional term ($`\alpha _s^4`$) in the TPS. As we shall see, these problems can be avoided with the use of Borel resummation.

A crucial new element of our analysis comes with the calculation of the nonperturbative contribution. Aside from the perturbative part, an important contribution to the GLS sum rule comes from the nonperturbative, hadronic matrix element of the twist–four operator. Being nonperturbative, this contribution is usually fitted using the QCD sum rule calculation. Recently it was proposed by one of us , motivated by an observation that the nonperturbative amplitudes in lower dimensional solvable models have a simple analyticity in complex coupling plane, that these higher–twist contributions can in principle be calculated from the Borel resummation of the perturbation series. The proposal was based on the conjecture that the higher–twist contributions have cuts only along the positive real axis in the complex coupling plane, which allows to relate the real part of the nonperturbative amplitude to its perturbatively calculable imaginary part. This scheme was shown to work well in model field theories. When applied to some of the solvable lower dimensional theories, it allowed the associated nonperturbative amplitudes to be accurately calculated from the first terms of the perturbation theory in the respective theories.

From our analysis we obtain for the strong coupling parameter the central value $`\alpha _s(M_Z)0.117`$. Compared to the corresponding CCFR central value $`\alpha _s(M_Z)=0.114`$, our value is closer to the world average $`\alpha _s(M_Z)0.118`$. The main improvement comes from the correct incorporation of the renormalon singularity on the Borel amplitude and the calculation of the nonperturbative contribution.

The paper is organized in the following way. In Sec. II we describe the resummation method, incorporating in it the known structure of the leading IR renormalon and the nonperturbative part, as well as the conformal mapping. Section III contains the numerical application of the method to the GLS sum rule, leading to predictions for $`\alpha _s(M_Z)`$. In Sec. IV we compare our predictions with those of other methods. In Sec. V we discuss some general features of the Borel resummation and the OPE approaches to understand precisely where the two methods deviate from each other. In Sec. VI we summarize our results and present conclusions.

## II The method

In this section we give an overview of our method used for the QCD analysis of the GLS sum rule. Its implementation in detail will be given in the following section. The QCD correction $`\mathrm{\Delta }(Q^2)`$ to the GLS sum rule is defined by

$$_0^1𝑑xF_3^{\nu N}(x,Q^2)=3(1\mathrm{\Delta }(Q^2)),$$
(2)
where $`F_3^{\nu N}`$ is the non-singlet deep inelastic scattering (DIS) structure function in $`\nu N`$ scattering. Here we shall ignore the target mass correction since it is irrelevant for our present discussion, but it will be included in the numerical analysis presented in the next section.

We first begin by reviewing the old but important problem  with the conventional QCD formulation of $`\mathrm{\Delta }(Q^2)`$ in OPE framework, which is widely used in data analysis. This problem is not confined to the GLS sum rule, but generic to any perturbative OPE formulations. $`\mathrm{\Delta }(Q^2)`$ in OPE up to twist–four operator is given by

$$\mathrm{\Delta }(Q^2)=W_0(\alpha _s(Q))+W_1(\alpha _s(Q))\frac{O}{Q^2},$$
(3)
where $`\alpha _s(Q)`$ is the strong coupling constant, and $`O`$ is the reduced nucleonic matrix element of the twist–four operator that was first derived in Ref.

$`O_\mu `$ $`=`$ $`\overline{u}\stackrel{~}{G}_{\mu \nu }\gamma ^\nu \gamma _5u+\overline{d}\stackrel{~}{G}_{\mu \nu }\gamma ^\nu \gamma _5d,`$ (4)
$`\stackrel{~}{G}^{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{2}}ϵ^{\mu \nu \alpha \beta }G_{\alpha \beta }^a{\displaystyle \frac{\lambda ^a}{2}},P|O_\mu |P_{\mathrm{spin}\mathrm{averaged}}=2p_\mu O.`$ (5)
Here, $`\lambda ^a`$ are the usual Gell-Mann matrices and we used the notations of Ref. . Throughout the article we shall consider only the twist–four contribution as the nonperturbative effect, and ignore higher twist contributions since they are believed to be small.

In conventional QCD analysis the Wilson coefficients $`W_i,(i=0,1)`$, in (3) are taken from the finite order, perturbative QCD calculation in an usual renormalization scheme, say $`\overline{\mathrm{MS}}`$ scheme, and the reduced matrix element $`O`$ from data fitting or QCD sum rule calculations, etc. However, this scheme is, in principle, fundamentally flawed, since perturbatively the Wilson coefficients are not well–defined. In perturbative calculation of the Wilson coefficients the quantum fluctuations of all energy scale contribute, and in particular at large orders the contribution from the far infrared regimes, where perturbative QCD should fail, is large and gives rise to a same sign, factorially growing large order behavior . Thus, without some kind of resummation of the divergent perturbation series, the Wilson coefficients are not well–defined, and therefore neither is the OPE (3). As a consequence, no well–defined meaning that is independent of the definition of the Wilson coefficient $`W_0`$, can be assigned to the matrix element $`O`$ .

One could in principle introduce an infrared cutoff $`\mu ^2(Q^2)`$ in the perturbative calculation of $`W_0`$, and regard the twist–four term to contain all the low momentum $`k^2K^2<\mu ^2`$ contributions . This would remove the infrared renormalon and make the perturbation expansion for $`W_0`$ convergent, and the twist–four term to be $`\mu ^2`$ dependent . However, to do this in practice is impossible, because it is impossible to compute arbitrarily complicated Feynman amplitudes of arbitrarily high power, and in particular their small momentum contributions, as stressed by the authors of Ref. .

The problem discussed so far is not of academic nature only, but has important practical implications. One might still think the OPE (3), with finite order perturbative Wilson coefficients, is a good approximation scheme, since at any rate the perturbative Wilson coefficients can be regarded as a good approximation at reasonably small values of the strong coupling constant. Actually, this would be the case, provided that the nonperturbative, higher–twist effect is far larger than the ignored next higher perturbative term in $`W_0`$. In practice, however, this condition is not supported by data analysis. For instance, if this were the case, we would expect little variation in the fitted values for the twist–four contribution over the order of perturbation in $`W_0`$. But the variation is not small at all. In the QCD sum rule calculation using the next–to–leading order (NLO) $`W_0`$ , the twist–four contribution was found not to be small, roughly equal to the perturbative correction at $`Q=1\mathrm{GeV}`$. But, it was observed in Refs.  that the twist–four contribution virtually vanishes when fitted with the next–to–next–to–leading order (NNLO) $`W_0`$ against the parton distribution functions extracted from experiments, and most of the twist–four contribution extracted in the NLO fit can be accounted for by the perturbative NNLO contributions. This can be interpreted as a clear manifestation of the inherent ambiguity of the OPE approach. Moreover, this tendency of strong dependence of higher–twist contribution on the order of perturbation appears not to be special to the GLS sum rule, but generic. The recent new estimate of the gluon condensate , from fitting the vectorial spectral function of hadronic tau decay using NNLO Adler function, gives a small central value of only $`1/3`$ that of the original QCD sum rule estimate  which uses the leading order perturbation. These examples strongly indicate that the inherent ambiguity of the OPE can have important consequences in practical applications, and demands a careful treatment.

Borel resummation resolves this problem, which proceeds as follows . The perturbation series for $`W_0`$

$$W_0=\underset{0}{\overset{\mathrm{}}{}}w_na(Q)^{n+1}\left[a(Q)\alpha _s(Q)/\pi \right],$$
(6)
which is, being of same sign at large ordersHere we ignore, for the moment, the UV renormalons, which give rise to sign-alternating large order behaviors. Being Borel summable, the UV renormalons can in principle be treated separately, and do not affect our discussion in any essential way., non-Borel resummable at physical, positive coupling $`a(Q)`$. So it is first Borel resummed at negative $`a(Q)`$ which yields a Borel resummed amplitude $`\mathrm{\Delta }_\mathrm{P}(a(Q))`$. Then to obtain a Borel resummed physical amplitude one may analytically continue $`\mathrm{\Delta }_\mathrm{P}`$ to positive $`a(Q)`$ in the complex coupling plane. This gives for $`a(Q)>0`$

$$\mathrm{\Delta }_\mathrm{P}(a(Q)\pm \mathrm{i}ϵ)=\frac{1}{\beta _0}_{0\pm \mathrm{i}ϵ}^{\mathrm{}\pm \mathrm{i}ϵ}𝑑be^{b/\beta _0a(Q)}\stackrel{~}{W}_0(b)$$
(7)
with

$$\stackrel{~}{W}_0(b)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{w_n}{n!}(b/\beta _0)^n.$$
(8)
Inserted for normalization convenience, $`\beta _0`$ is the one–loop coefficient of the QCD $`\beta `$-function. Eq. (7) shows explicitly that in case of existing singularities on the line of integration, i.e., IR renormalons, the real part of the Borel integral $`\mathrm{\Delta }_\mathrm{P}`$ is the (generalization of the) Cauchy Principal Value. The Borel transform (8) is believed to be convergent on the unit disk $`|b|<1`$, and is known to have a branch cut along the positive real axis beginning at $`b=1`$. Near the branch cut it behaves as

$$\stackrel{~}{W}_0(b)=\frac{C}{\mathrm{\Gamma }(\nu )}\beta _0^{1+\nu }(1b)^{1\nu }(1+O(1b))+\left(\text{Analytic part}\right)$$
(9)
with

$$\nu =(\beta _1/\beta _0\gamma _2)/\beta _0,$$
(10)
where $`\beta _1`$ and $`\gamma _2`$ are respectively the two–loop coefficient of the QCD $`\beta `$–function and the one–loop coefficient of the anomalous dimension of the twist–four operator appearing in the OPE (3). In our case $`\nu `$ is positive; for instance, $`\nu =32/81`$ when three quark flavors ($`n_f=3`$) are active .§§§ According to Ref. : $`\gamma _2=(N_c1/N_c)/3=8/9`$. Our convention for parameters $`\beta _j`$ (and $`c_j\beta _j/\beta _0`$) is specified by Eq. (24) in the next Section. For $`n_f=3`$, we have $`\beta _0=9/4`$, $`\beta _1=4`$. Therefore, $`\nu =32/81`$. Note that the Borel transform beyond the convergence disk of the series (8) can be obtained by analytic continuation. The analytic part, which is analytic around the singularity, is not calculable, but the residue $`C`$, a real number, is calculable perturbatively .

Because of the branch cut the Borel integrals in (7) develop imaginary parts beginning at $`b=1`$. Since the QCD correction $`\mathrm{\Delta }(Q^2)`$ must be real, clearly $`\mathrm{\Delta }_\mathrm{P}`$ alone cannot reproduce the true amplitude. There must be something else. Precisely at this point, the nonperturbative amplitude, denoted $`\mathrm{\Delta }_{\mathrm{NP}}(a(Q))`$, comes to the rescue, which cancels the imaginary parts of $`\mathrm{\Delta }_\mathrm{P}`$, rendering the sum of the two to be real. This $`\mathrm{\Delta }_{\mathrm{NP}}`$ can be shown to be directly related to the twist–four contribution in the OPE (3) , with its overall form governed by the associated RG equation, and may be regarded as a refinement of the latter, now free from the inherent ambiguity of the OPE. The QCD correction may now be written as

$`\mathrm{\Delta }(Q^2)`$ $`=`$ $`\mathrm{\Delta }_\mathrm{P}(a\pm \mathrm{i}ϵ)+\mathrm{\Delta }_{\mathrm{NP}}(a\pm \mathrm{i}ϵ)`$ (11)
$`=`$ $`\mathrm{Re}[\mathrm{\Delta }_\mathrm{P}(a\pm \mathrm{i}ϵ)]+\mathrm{Re}[\mathrm{\Delta }_{\mathrm{NP}}(a\pm \mathrm{i}ϵ)].`$ (12)
Either sign can be taken because the result is independent of the sign chosen.

What can we tell about the nonperturbative amplitude $`\mathrm{\Delta }_{\mathrm{NP}}`$? Its imaginary part at a positive coupling is certainly calculable from the perturbation theory because it is essentially the imaginary part of $`\mathrm{\Delta }_\mathrm{P}`$, albeit with opposite sign, which is calculable in principle from the perturbation theory. However, as we see in (12), what we need is the real part, which is certainly not directly calculable.

There is, however, an intriguing possibility of perturbative calculation of the real part. From the perspective of Borel resummation the sole reason for the introduction of the nonperturbative amplitude above was to cancel the imaginary parts arising from the analytic continuation of $`\mathrm{\Delta }_\mathrm{P}`$ to physical coupling in the complex coupling plane. With the Borel integral (7) and the singular Borel transform (9), we can see that $`\mathrm{\Delta }_\mathrm{P}(a(Q))`$, which by definition can have a singularity only along the positive real axis in the complex $`a(Q)`$ plane, has a branch cut of the form (see Appendix A)

$$C(a(Q))^\nu e^{1/\beta _0a(Q)}[\mathrm{\hspace{0.17em}1}+O(a)].$$
(13)
The nonperturbative amplitude $`\mathrm{\Delta }_{\mathrm{NP}}`$ should cancel the imaginary part arising from this at positive coupling. The simplest functional form for $`\mathrm{\Delta }_{\mathrm{NP}}`$ that can achieve this purpose can be obtained by postulating that, as has $`\mathrm{\Delta }_\mathrm{P}`$, the nonperturbative amplitude have a branch cut only along the positive real axis in the coupling plane. This then leads to the following conjecture for $`\mathrm{\Delta }_{\mathrm{NP}}`$ :

$$\mathrm{\Delta }_{\mathrm{NP}}(a(Q))=C(a(Q))^\nu e^{1/\beta _0a(Q)}[\mathrm{\hspace{0.17em}1}+O(a)]$$
(14)
This has a very important implication because it allows to relate the real part of the nonperturbative amplitude to the calculable imaginary part. In this paper we will adopt this conjecture, which appears very plausible, at least to us, and as will be mentioned shortly it is supported by some lower dimensional solvable models. Now with (12) and (14) we can write the QCD correction with only the calculable $`\mathrm{\Delta }_\mathrm{P}`$ as

$$\mathrm{\Delta }_\mathrm{P}(Q^2)+\mathrm{\Delta }_{\mathrm{NP}}(Q^2)=\left\{\mathrm{Re}\mathrm{cot}(\nu \pi )\mathrm{Im}\right\}[\mathrm{\Delta }_\mathrm{P}(a(Q)\pm \mathrm{i}ϵ)].$$
(15)
This equation is the basis of our numerical analysis in the following section. In Appendix A we present, for reference, some explicit formulas leading to (13)–(15).

The argument that led to the determination of the nonperturbative amplitude above did not rely on any special property of OPE, but only on the general property of Borel resummation of a same sign divergent perturbation theory. Thus it can be applied to any perturbation theory with a same sign large order behavior. Application of this scheme to the lower dimensional, solvable models such as the double well potential and the two-dimensional nonlinear $`\sigma `$ model in large $`N`$ limit, allowed an accurate calculation of the associated nonperturbative amplitude using only the first terms of the perturbation expansion in the respective models . Moreover, interestingly, the numerical estimate of the gluon condensate from this scheme applied to the Adler function gives a value virtually identical to the new estimate .

Given the QCD correction in the form (15) with $`\mathrm{\Delta }_\mathrm{P}(a\pm \mathrm{i}ϵ)`$ given by (7), our aim is now to describe the Borel transform $`\stackrel{~}{W}_0`$ as accurately as possible in the primary interval using the known first terms of the perturbation series. For this, following the outline in Introduction, we write the Borel integral (7) as

$$\mathrm{\Delta }_\mathrm{P}(a(Q)\pm \mathrm{i}ϵ)=\frac{1}{\beta _0}_{0\pm \mathrm{i}ϵ}^{\mathrm{}\pm \mathrm{i}ϵ}𝑑be^{b/\beta _0a(Q)}\frac{R(b)}{(1b)^{1+\nu }}$$
(16)
with $`R(b)`$ now defined by

$$R(b)(1b)^{1+\nu }\stackrel{~}{W}_0(b)$$
(17)
and $`\nu `$ given in (10) ($`\nu =32/81`$). This step is expected to greatly improve the perturbative description of the Borel transform in the primary interval because it implements the renormalon singularity correctly, and renders us to deal with a much softer singularity. The singularity of $`R(b)`$ at $`b=1`$ is a branch cut and thus softer than that of $`\stackrel{~}{W}_0(b)`$.

We can obtain further improvement by use of a conformal mapping that exploits the known locations of the singularities of $`\stackrel{~}{W}_0(b)`$. The latter is known to have renormalon singularities at non-zero integer values of $`b`$ on the real axis . To speed up the convergence of the perturbative expansion of $`R(b)`$, we may push the singularities, save the unavoidable first IR renormalon, as far away from the origin as possible. This way we can reduce the influence of the renormalon singularities and make $`R(b(w))`$ smoother around the primary interval. One such a mapping we consider is

$$w(b)=\frac{\sqrt{1+b}\sqrt{1b/2}}{\sqrt{1+b}+\sqrt{1b/2}}b(w)=\frac{8w}{(3w^22w+3)}.$$
(18)
This maps the first IR renormalon to $`w=1/3`$ and all other singularities to the unit circle $`|w|=1`$. We expect this mapping combined with the implementation (17) to provide an optimized environment for the Borel integral. In the mapped plane the Borel integral now assumes the form

$$\mathrm{\Delta }_\mathrm{P}(a(Q)\pm \mathrm{i}ϵ)=\frac{1}{\beta _0}_{𝒞_\pm }𝑑we^{b(w)/\beta _0a(Q)}\frac{db(w)}{dw}\frac{R(b(w))}{(1b(w))^{1+\nu }},$$
(19)
where one of the integration contours $`𝒞_{}`$ is shown in Fig. 1. Again, since the answer is independent of the sign chosen, either contour can be taken. In the next section we perform a numerical analysis of our implementation of the QCD correction: Eqs. (15) and (19).

## III Numerical analysis

In this Section we will apply the method described in the previous Section to the Gross–Llewellyn Smith (GLS) sum rule, deducing values of the QCD coupling parameter $`\alpha _s^{\overline{\mathrm{MS}}}(M_Z^2)`$ from the GLS values extracted from experiments. In the first Subsection, we will present the resummed expression for the contributions of the three massless quark flavors. As a byproduct of the obtained expression, we will obtain an estimate of the next–to–next–to–next–to–leading ($`\mathrm{N}^3\mathrm{L}`$) coefficient $`w_3`$ of the perturbative expansion. In the subsequent Subsection, we will include the effects of the massive fourth quark flavor ($`c`$–quark) to the GLS observable. In the last Subsection, the available measured GLS values will be confronted with our resummed expression and values of $`\alpha _s^{\overline{\mathrm{MS}}}(M_Z^2)`$ will be extracted.

### A GLS – the massless $`n_f=3`$ part

The GLS quantity $`M_3(Q^2)`$ is the following integral (first moment) of the charged–current non–polarized DIS structure function $`F_3(F_3^{\nu p}+F_3^{\nu n})/2`$ over the Bjorken parameter $`x`$:

$$M_3(Q^2)=\frac{1}{3}_0^1𝑑xF_3(x;Q^2)\frac{\zeta ^2}{x^2}\left[1+2\left(1+\frac{4m_N^2x^2}{Q^2}\right)^{1/2}\right],$$
(20)
where $`\zeta 2x/(1+\sqrt{1+4m_N^2x^2/Q^2})`$ is the Nachtmann variable , and $`m_N`$ is the nucleon mass. The quantity $`M_3`$ is the first Nachtmann moment of $`F_3`$ which absorbs all the kinematical power corrections (target mass corrections: TMC) $`(m_N^2/Q^2)^n`$. The quantity $`Q^2=q^2`$ is the virtuality of the exchanged gauge boson; $`Q^2`$ characterizes the typical process momenta of the observable $`M_3(Q^2)`$. When expanding the integrand in Eq. (20) in powers of $`m_N^2/Q^2`$, we can rewrite

$$M_3(Q^2)=_0^1𝑑xF_3(x;Q^2)\left[1\frac{2}{3}\frac{m_N^2x^2}{Q^2}+\left(\frac{m_N^2x^2}{Q^2}\right)^2+𝒪\left(\frac{m_N^6}{Q^6}\right)\right].$$
(21)
The above GLS moment can be written in the form

$$M_3(Q^2)3\left(1\mathrm{\Delta }(Q^2)\right),$$
(22)
where the “canonical” quantity $`\mathrm{\Delta }(Q^2)`$ has the following power expansion $`W_0`$ \[cf. Eq. (3)\]:

$$\mathrm{\Delta }(Q^2)W_0(Q^2)=a(1+w_1a+w_2a^2+w_3a^3+\mathrm{}).$$
(23)
Here, $`a\alpha _s(\mu ^2;c_2,c_3,\mathrm{})/\pi `$ is the QCD coupling parameter with a given choice of the renormalization scale (RScl) $`\mu ^2`$ and the renormalization scheme (RSch) parameters $`c_j`$ ($`j2`$). The evolution of $`a`$ with the RScl is governed by the renormalization group equation (RGE)

$$\frac{a(\mu ^2;c_2,c_3,\mathrm{})}{\mathrm{ln}\mu ^2}=\beta _0a^2(1+c_1a+c_2a^2+c_3a^3+\mathrm{}),$$
(24)
where the parameters $`c_j\beta _j/\beta _0`$ ($`j2`$) characterize the choice of the RSch, and $`\beta _0`$ and $`c_1\beta _1/\beta _0`$ are universal constants. $`\beta _0=(112n_f/3)/4`$, $`c_1=(10238n_f/3)/(16\beta _0)`$, where $`n_f`$ is the number of active quark flavors. The evolution of parameter $`a`$ with $`c_j`$’s ($`j2`$) is governed by analogous differential equations, which follow from the Stevenson equation (see Appendix A of Ref. ). The next–to–leading (NL) coefficient $`w_1`$ has been calculated in Ref. , and the NNL coefficient $`w_2`$ in Ref. . At the specific RScl $`\mu ^2=Q^2`$ and $`\overline{\mathrm{MS}}`$ RSch, and when the number of the active quark flavors is $`n_f=3`$, these coefficients have the values The superscript ‘(0)’ in $`w_j^{(0)}`$ denotes the special RScl–choice $`\mu ^2=Q^2`$ and $`\overline{\mathrm{MS}}`$ RSch.

$`w_1^{(0)}w_1(\mu ^2=Q^2)`$ $`=`$ $`3.58333,`$ (25)
$`w_2^{(0)}w_2(\mu ^2=Q^2;c_2^{\overline{\mathrm{MS}}})`$ $`=`$ $`20.2153+w_2^{(0)}(l.l.),`$ (26)
$`w_2^{(0)}(l.l.)`$ $`=`$ $`1.23954.`$ (27)
In the NNL coefficient, the small “light–by–light” part was separated off. The “light–by–light” contribution should not be included in resummations of $`\mathrm{\Delta }(Q^2;\mathrm{P})`$, as will be argued at the end of this Subsection.

The Borel integral (19), which will be together with (15) the basis for our resummation, is independent of the choice of the RScl $`\mu ^2`$ and the RSch ($`c_2,c_3,\mathrm{}`$) used in the integrand.\**\**\** The location of the renormalon pole $`b=1`$ and the power $`\nu =32/81`$ (10) are independent of the choice of the RScl and RSch. Thus, we can rewrite it as (see also Fig. 1)

$`\mathrm{\Delta }_\mathrm{P}(Q^2\pm iϵ)`$ $`=`$ $`e^{\pm i\varphi }{\displaystyle \frac{1}{\beta _0}}{\displaystyle _0^1}𝑑x{\displaystyle \frac{db(w)}{dw}}\mathrm{exp}\left[{\displaystyle \frac{b(w)}{\beta _0a(\mu ^2)}}\right]{\displaystyle \frac{R(b(w);\mu ^2/Q^2)}{(1b(w))^{113/81}}}|_{w=xe^{\pm i\varphi }},`$ (28)
At this stage we can, as a byproduct of formula (17), obtain an estimate of the yet unknown $`\mathrm{N}^3\mathrm{L}`$ coefficient $`w_3`$ appearing in the perturbative expansion (23) . If working with the RScl $`\mu ^2=Q^2`$ and in the $`\overline{\mathrm{MS}}`$ scheme, the expansion of $`R(b(w))`$ in powers of $`w`$ gives us:

$`R(b(w))`$ $`=`$ $`1+0.526749w+0.709369w^2`$ (30)
$`+(43.2574+0.277464w_3^{(0)})w^3+𝒪(w^4),`$
where we excluded the “light–by–light” part (27) of the NNL coefficient $`w_2^{(0)}`$ (26). Looking at the coefficients appearing in (30), it is reasonable to expect that the coefficient $`R_3`$ at $`w^3`$ is $`|R_3|1`$. If we assume $`R_3=1\pm 1`$, we obtain a rather stringent estimate $`w_3^{(0)}=159.5\pm 3.6`$. If we adopt a more cautious assumption $`|R_3|\stackrel{<}{}10^1`$, we obtain $`w_3^{(0)}160\pm 30`$. If we apply Padé approximant (PA) $`[1/1]_R(w)`$ to the expansion (30) and re–expand it back in powers of $`w`$ up to $`w^3`$, we obtain an estimate $`w_3^{(0)}=159.3`$. On the other hand, if applying the $`[1/2]_R(w)`$ and demanding $`w_{\mathrm{pole}}=1`$ (i.e., $`b=+2`$, $`\mathrm{IR}_2`$), the prediction is $`w_3^{(0)}=158.5`$; if demanding $`w_{\mathrm{pole}}=1`$ (i.e., $`b=1`$, $`\mathrm{UV}_1`$), the prediction is $`w_3^{(0)}=157.0`$. Very similar estimates are obtained if we do not apply the conformal transformation $`b(w)`$ (18). Therefore, we will adopt the following estimate for $`w_3^{(0)}`$

$$w_3^{(0)}\left[w_3(\mu ^2=Q^2;\overline{\mathrm{MS}})\right]=158\pm 30.$$
(31)
We emphasize that this estimate excludes the “light–by–light” contributions, which are assumed to be suppressed at the $`\mathrm{N}^3\mathrm{L}`$ order. The exclusion of the “light–by–light” contributions reduces the (NNL) perturbative expansion of the (non-polarized) GLS sum rule to that of the Bjorken polarized sum rule (BjPSR) . It is interesting that the effective charge method (ECH)  and the (TPS) principle of minimal sensitivity (PMS)  predict $`w_3^{(0)}130`$ , based on the assumption that $`c_3^{\mathrm{ECH}}c_3^{\mathrm{PMS}}c_3^{\overline{\mathrm{MS}}}`$. Further, an RScl– and RSch–invariant method  that is somewhat related to the PA and PMS approaches, also predicts $`w_3^{(0)}130`$. This is at the lower end of our new estimate (31). We thus conclude that the explicit (and exact) structure of the leading IR renormalon of the GLS–observable (BjPSR–observable) in the Borel plane, as given in Eq. (17), is responsible for the somewhat higher estimate of $`w_3^{(0)}`$ in comparison with the ECH, PMS, and PA–related methods of resummation.

We add here a few remarks on the question of the “light–by–light” contributions. The power expansion of the massless part of the GLS sum rule to $`\mathrm{N}^3\mathrm{L}`$ order $`a^3`$ (see Refs. ), in the $`\overline{\mathrm{MS}}`$ scheme and at the RScl $`\mu ^2=Q^2`$, is given by

$`({\displaystyle \frac{3}{4}}C_F)W_0(Q^2)`$ $`=`$ $`a({\displaystyle \frac{3}{4}}C_F)+a^2C_F({\displaystyle \frac{21}{32}}C_F{\displaystyle \frac{23}{16}}C_A+{\displaystyle \frac{1}{4}}n_f)`$ (35)
$`+a^3[{\displaystyle \frac{3}{128}}C_F^3+({\displaystyle \frac{1241}{576}}{\displaystyle \frac{11}{12}}\zeta _3)C_F^2C_A+({\displaystyle \frac{5437}{864}}+{\displaystyle \frac{55}{24}}\zeta _5)C_FC_A^2`$
$`+({\displaystyle \frac{133}{1152}}{\displaystyle \frac{5}{24}}\zeta _3)C_F^2n_f+({\displaystyle \frac{3535}{1728}}+{\displaystyle \frac{3}{8}}\zeta _3{\displaystyle \frac{5}{12}}\zeta _5)C_FC_An_f{\displaystyle \frac{115}{864}}C_Fn_f^2`$
$`+n_f{\displaystyle \frac{d^{abc}d^{abc}}{N_c}}({\displaystyle \frac{11}{192}}+{\displaystyle \frac{1}{8}}\zeta _3)]+𝒪(a^4),`$
where the Casimir coefficients for QCD ($`N_c=3`$) are $`C_F=4/3`$, $`C_A=3`$, and the group–theoretical factor appearing in the last term in (35) is

$$d^{abc}d^{abc}=\frac{40}{3}.$$
(36)
This group–theoretical factor is not present in the calculation of the GLS sum rule up to the two–loop ($`a^2`$) order, and it appears for the first time at the three–loop ($`a^3`$) order. This term is called “light–by–light,” it corresponds to diagrams with a new topology involving exchange of three gluons. In our calculation we will add this “light–by–light” contribution $`\mathrm{\Delta }_{l.l.}(Q^2)`$, given by the last term in Eq. (35), as a separate term not included in our resummation approach. The reason for this is the following: Resummation approaches cannot be expected to predict (and resum) those higher order terms which are characterized by new higher order group–theoretical factors, when we have only one such term ($`a^3`$) explicitly available. We assume that such resummed “light–by–light” contributions are small, comparable to the quite small $`a^3`$ “light–by–light” contribution in (35). Similar considerations can be found in Refs. , in cases of various observables and beta functions.

### B Inclusion of the massive quark flavor ($`c`$) contribution and nuclear corrections

In the approximation of massless quarks, the calculation of the $`\mathrm{N}^3\mathrm{L}`$ ($`a^3`$) QCD correction to the GLS sum rule has been carried out in Ref. . The inclusion of the heavy ($`c`$–)quark contributions is important at the precision level at which we are working. In addition, it is important in order to estimate the scales $`Q^2`$ where the ($`c`$–)quarks can be treated as massless or massive quarks. This point is important because it indicates which number of light flavors $`n_f`$ should be used in the (resummed) “massless” part of the the perturbation series. The calculation of the heavy flavor contribution to the GLS sum rule up to the second order ($`a^2`$) was performed in Ref.  and discussed in . According to this approach: (a) the quarks $`u,d,s`$ are massless and result in the dominant $`n_f=3`$ massless QCD contribution to the GLS sum rule $`\mathrm{\Delta }(Q^2)`$; (b) for $`Q^22`$$`4\mathrm{GeV}^2`$, the massive flavor is the $`c`$–quark and it contributes a relatively small correction to the aforementioned $`n_f=3`$ massless contribution. The other heavy quark flavors are ignored due to the strong suppression by the mixing angles in the Cabibbo–Kobayashi–Maskawa matrix and by the small values of $`Q^2`$.

The heavy ($`c`$–)flavor correction contributions to $`\mathrm{\Delta }(Q^2)`$ are

$`\mathrm{\Delta }_c(Q^2)`$ $`=`$ $`\mathrm{\Delta }_c^{(1)}+\mathrm{\Delta }_c^{(2)},`$ (37)
$`\mathrm{\Delta }_c(Q^2)^{(1)}`$ $`=`$ $`\left[{\displaystyle \frac{1}{3(1+\xi )}}a(\mu ^2){\displaystyle \frac{C_F}{4}}\left\{{\displaystyle \frac{1}{1+\xi }}+2{\displaystyle \frac{\mathrm{ln}(1+\xi )}{1+\xi }}\right\}\right]\mathrm{sin}^2\theta _c,`$ (38)
$`\mathrm{\Delta }_c(Q^2)^{(2)}`$ $`=`$ $`a(\mu ^2)^2{\displaystyle \frac{C_FT_F}{16}}[({\displaystyle \frac{1}{105}}\xi ^2+{\displaystyle \frac{16}{45}}\xi )\mathrm{ln}\xi `$ (42)
$`+{\displaystyle \frac{1}{\lambda ^4}}\left({\displaystyle \frac{2}{105}}\xi +{\displaystyle \frac{2783}{315}}+{\displaystyle \frac{6740}{63}}{\displaystyle \frac{1}{\xi }}+{\displaystyle \frac{137552}{315}}{\displaystyle \frac{1}{\xi ^2}}+{\displaystyle \frac{62528}{105}}{\displaystyle \frac{1}{\xi ^3}}\right)`$
$`{\displaystyle \frac{1}{\lambda ^5}}({\displaystyle \frac{1}{105}}\xi ^2+{\displaystyle \frac{142}{315}}\xi +{\displaystyle \frac{494}{63}}+{\displaystyle \frac{1516}{21}}{\displaystyle \frac{1}{\xi }}+{\displaystyle \frac{23024}{63}}{\displaystyle \frac{1}{\xi ^2}}+{\displaystyle \frac{298432}{315}}{\displaystyle \frac{1}{\xi ^3}}`$
$`+{\displaystyle \frac{102656}{105}}{\displaystyle \frac{1}{\xi ^4}})\mathrm{ln}({\displaystyle \frac{\lambda +1}{\lambda 1}}){\displaystyle \frac{20}{3}}{\displaystyle \frac{1}{\xi ^2}}\mathrm{ln}^2({\displaystyle \frac{\lambda +1}{\lambda 1}})],`$
where $`\xi =Q^2/m_c^2`$, $`\lambda =\sqrt{1+4m_c^2/Q^2}`$, $`C_F=4/3`$, $`T_F=1/2`$. The heavy flavor corrections (38)–(42) will not be included in the resummation procedure for $`\mathrm{\Delta }`$, but will be added separately.

In addition to the discussed heavy flavor contributions, there are also nuclear corrections to the GLS sum rule, due to the nuclear effects in the $`F_3`$ structure functions. These effects were calculated in Ref.  and were found for the iron target (used by the CCFR Collaboration) to be small:

$$\mathrm{\Delta }_{\mathrm{Fe}}(Q^2)\frac{410^3\mathrm{GeV}^2}{Q^2}.$$
(43)
This contribution, at $`Q^2=2`$$`3\mathrm{GeV}^2`$, turns out to be in its magnitude smaller than the heavy flavor effects by a factor of 3–4, but larger than the “light–by–light” contribution by a factor of about 2.

### C Extraction of $`\alpha _s`$ values

The values of the GLS sum rule $`M_3(Q^2)`$, for various specific $`Q^2`$, have been obtained from the experimentally extracted values of the structure function $`F_3(x;Q^2)`$, by the Fermilab CCFR Collaboration . They have large systematic experimental uncertainties, primarily because of the uncertainties in the normalization of $`xF_3`$ and in the integration in the regions $`x1`$ and $`x>0.5`$. The systematic uncertainties are somewhat smaller, comparable to the statistical experimental uncertainties, only when the values of the exchanged boson virtuality are $`Q^2=2.00,3.16\mathrm{GeV}^2`$ (see their Table III). These values, including the target mass correction terms of (20)–(21), as well as the nuclear correction term (43), are†††††† In the entries of Table III of Ref. , the target mass corrections are included, but the nuclear corrections neglected.

$`Q^2=2.00\mathrm{GeV}^2:M_3(Q^2)`$ $`=`$ $`2.49\pm 0.13\mathrm{\Delta }(Q^2)=0.168\pm 0.043,`$ (44)
$`Q^2=3.16\mathrm{GeV}^2:M_3(Q^2)`$ $`=`$ $`2.55\pm 0.12\mathrm{\Delta }(Q^2)=0.149\pm 0.039.`$ (45)
Here we added in quadrature the statistical and the systematic experimental uncertainties. In addition, the values $`\mathrm{\Delta }(Q^2)`$ in (44)–(45) have the aforementioned nuclear correction contribution (43) subtracted out; however, $`\mathrm{\Delta }(Q^2)`$ still includes the massless (perturbative and nonperturbative) contribution (15), the light–by–light contribution $`w_2(l.l.)a^3`$ \[cf. (27)\] and the $`c`$–quark contribution (37)–(42).

For a given value of $`\alpha _s(Q^2,\overline{\mathrm{MS}})`$, our expression for $`\mathrm{\Delta }(Q^2)=\mathrm{\Delta }_\mathrm{P}(Q^2)+\mathrm{\Delta }_{\mathrm{NP}}(Q^2)+\mathrm{\Delta }_{l.l.}(Q^2)+\mathrm{\Delta }_c(Q^2)`$, given in the previous Subsections, still has the freedom (ambiguity) of the choice of the RScl $`\mu ^2`$ and the RSch.

In the following analysis, we fix first the RSch to be $`\overline{\mathrm{MS}}`$, with the beta function being the PA $`\beta (x)=[2/3](x)`$ based on the $`\mathrm{N}^3\mathrm{L}`$ TPS of the $`\beta _{\overline{\mathrm{MS}}}(x)`$. This PA choice of the beta function has reasonable behavior in the region of large $`x\alpha _s/\pi `$ and has been used in Refs.  in the analyses of low energy QCD observables. We will comment later on how the results change when we use $`\mathrm{N}^3\mathrm{L}`$ TPS of the $`\beta _{\overline{\mathrm{MS}}}(x)`$, and when we change the RSch even more drastically.

The RScl–dependence is yet another source of the theoretical ambiguity. In Figs. 2 and 3 we show the dependence of the predicted values of $`\mathrm{\Delta }(Q^2)`$ on the RScl $`\mu ^2`$ at the NNL and $`\mathrm{N}^3\mathrm{L}`$ level, at $`Q^2=2.00,3.16\mathrm{GeV}^2`$, respectively, for given fixed representative values of $`\alpha _s(Q^2,\overline{\mathrm{MS}})`$. The NNL level means that we take for $`R(b(w))`$ in the Borel integrations in (28) and (15) the NNL TPS (quadratic polynomial in $`w`$), i.e., only the knowledge of $`w_1`$ and $`w_2`$ coefficients is used. The $`\mathrm{N}^3\mathrm{L}`$ level means that $`R(b(w))`$ is the $`\mathrm{N}^3\mathrm{L}`$ TPS (cubic polynomial in $`w`$), i.e., we use, in addition, $`w_3^{(0)}=158`$ according to the estimate (31). From Figs. 2,3 we see that the $`\mathrm{N}^3\mathrm{L}`$ expressions drastically improve the stability of the predictions under the change of the RScl. In the case $`Q^2=2.00\mathrm{GeV}^2`$, the $`\mathrm{N}^3\mathrm{L}`$ values of the massless quantity $`\mathrm{\Delta }_{\mathrm{P}+\mathrm{NP}}(Q^2)`$ achieve minimal $`\mathrm{ln}\mu ^2`$–sensitivity at $`\mu ^2/Q^23.3`$ : $`d\mathrm{\Delta }_{\mathrm{P}+\mathrm{NP}}(Q^2)/d\mathrm{ln}\mu ^21.4510^4`$. At $`Q^2=3.16\mathrm{GeV}^2`$, there is unfortunately no point of minimal RScl–sensitivity; the slope is negative and getting weaker when $`\mu ^2`$ increases. However, the RScl–sensitivity is very weak when $`\mu ^2/Q^2>2.5`$, and almost stabilizes when $`\mu ^2/Q^23.3`$ : $`d\mathrm{\Delta }_{\mathrm{P}+\mathrm{NP}}(Q^2)/d\mathrm{ln}\mu ^21.5510^3`$. Therefore, we will take, for definiteness, the RScl choice $`\mu ^2=3.3Q^2`$ in the case of both $`Q^2`$ values. Later on, we will comment on the ambiguity of our results when the RScl is varied.

Having fixed the RSch ($`\overline{\mathrm{MS}}`$, with $`[2/3]_\beta (x)`$ PA) and the RScl ($`\mu ^2=3.3Q^2`$), our expressions for $`\mathrm{\Delta }(Q^2)=\mathrm{\Delta }(Q^2,\mathrm{P}+\mathrm{NP}+l.l.+c)`$ become unambiguous functions of the input value of $`\alpha _s(Q^2;\overline{\mathrm{MS}})`$. Adjusting simply the latter value so that the experimental values (44)–(45) are achieved, we obtain:

$`Q^2=2.00\mathrm{GeV}^2:\alpha _s(Q^2;\overline{\mathrm{MS}})`$ $`=`$ $`0.348_{0.093}^{+\mathrm{?}}(\mathrm{exp}.)`$ (46)
$`Q^2=3.16\mathrm{GeV}^2:\alpha _s(Q^2;\overline{\mathrm{MS}})`$ $`=`$ $`0.305_{0.074}^{+0.132}(\mathrm{exp}.)`$ (47)
The central, upper and lower values here correspond to the pertaining experimental values of $`\mathrm{\Delta }(Q^2)`$ in (44)–(45). In the case $`Q^2=2.00\mathrm{GeV}^2`$, the upper bound for the coupling parameter cannot be obtained from the experimental data because the applied resummation method predicts values of $`\mathrm{\Delta }(Q^2)`$ which are always lower than the presently allowed experimental upper bound: $`\mathrm{\Delta }(Q^2)_{\mathrm{meth}.}0.196`$ \[$`<\mathrm{\Delta }(Q^2)_{\mathrm{max}\mathrm{exp}.}=0.168+0.043`$\]. The situation in the case $`Q^2=3.16\mathrm{GeV}^2`$ is somewhat similar, the experimental upper bound being now slightly below the maximal value allowed by the method. This situation is presented graphically in Figs. 4 and 5 which show $`\mathrm{\Delta }(Q^2)`$ as function of $`\alpha _s(Q^2)`$ as predicted by the applied method, for $`Q^2=2.00,3.16\mathrm{GeV}^2`$ cases, respectively. The central, upper and lower experimental bounds are included as straight dotted lines.

We now return to the question of the uncertainties of our predictions under the variation of the RScl and RSch. If we vary the RScl parameter $`\xi ^2\mu ^2/Q^2`$ around 3.3 across the interval $`1.5\xi ^25`$, the predictions for $`\alpha _s(Q^2;\overline{\mathrm{MS}})`$ vary by at most $`\pm 0.002`$. On the other hand, variation of the RSch leads to larger ambiguities. For example, if we repeat the calculation in the ‘t Hooft RSch ($`c_2=c_3=0`$, we take $`[2/3]_\beta (x)`$ for definiteness),‡‡‡‡‡‡ For comparison, $`c_2^{\overline{\mathrm{MS}}}=4.471`$ $`c_3^{\overline{\mathrm{MS}}}=20.99`$, for $`n_f=3`$. keeping $`\xi ^2=3.3`$, the predictions for $`\alpha _s(Q^2;\overline{\mathrm{MS}})`$ change by $`0.009`$ ($`Q^2=2.00\mathrm{GeV}^2`$) and $`0.005`$ ($`Q^2=3.16\mathrm{GeV}^2`$). We will regard these values as characteristic values for the RSch–uncertainties of our results. The replacement of the $`\overline{\mathrm{MS}}`$ $`[2/3]_\beta (x)`$ by the $`\mathrm{N}^3\mathrm{L}`$ TPS $`\overline{\mathrm{MS}}`$ beta function changes our predictions for $`\alpha _s(Q^2;\overline{\mathrm{MS}})`$ by only about $`0.001`$.

Yet another source of the theoretical uncertainty in our predictions may be the truncation in the ($`\mathrm{N}^3\mathrm{L}`$) TPS $`R(b(w))`$. We regard the uncertainty $`\pm 30`$ in the estimated value of $`w_3^{(0)}`$, Eq. (31), as the major source of the truncation uncertainty. This changes the prediction for $`\alpha _s(Q^2;\overline{\mathrm{MS}})`$ by only $`\pm 0.001`$.

The other sources of theoretical uncertainty come from the massive ($`c`$–)quark contributions presented in Subsec. III B, due to the uncertainties $`m_c=1.25\pm 0.10\mathrm{GeV}`$ and $`\mathrm{sin}(\theta _{\mathrm{Cabibbo}})=|V_{cd}|=0.224\pm 0.016`$ . The resulting uncertainties of the predictions for $`\alpha _s(Q^2;\overline{\mathrm{MS}})`$ are $`\pm 0.002`$ and $`\pm 0.002`$, respectively, when $`Q^2=2.00\mathrm{GeV}^2`$, and $`\pm 0.001`$ and $`\pm 0.001`$ when $`Q^2=3.16\mathrm{GeV}^2`$.

The final predictions for $`\alpha _s`$ are presented in Table I.

In the Table, we presented the central predictions for $`\alpha _s(Q^2;\overline{\mathrm{MS}})`$ when $`Q^2=2.00`$ and $`3.16\mathrm{GeV}^2`$, and the uncertainties of the predictions due to various sources. Further, we RGE–evolved these predictions to the canonical scale $`M_Z^2`$ and included the results in Table I. This RGE evolution was carried out by using the $`[2/3]`$ Padé approximant of the four–loop $`\overline{\mathrm{MS}}`$ TPS beta function, using the values of the four–loop coefficient $`c_3(n_f)`$  and the corresponding three–loop matching conditions  for the flavor thresholds.\**\**\** For details on the corresponding evolution uncertainties, we refer to Ref. . They include the variation when the \[2/3\] Padé form of the beta function is replaced by the TPS form. When adding in quadrature the various theoretical uncertainties, the predictions of Table I can be summarized as

$`Q^2=2.00\mathrm{GeV}^2:\alpha _s(Q^2;\overline{\mathrm{MS}})`$ $`=`$ $`0.348_{0.093}^{+\mathrm{?}}(\mathrm{exp}.)\pm 0.010(\mathrm{th}.),`$ (48)
$`\alpha _s(M_Z^2;\overline{\mathrm{MS}})`$ $`=`$ $`0.1166_{0.0115}^{+\mathrm{?}}(\mathrm{exp}.)\pm 0.0010(\mathrm{th}.);`$ (49)
and

$`Q^2=3.16\mathrm{GeV}^2:\alpha _s(Q^2;\overline{\mathrm{MS}})`$ $`=`$ $`0.305_{0.074}^{+0.132}(\mathrm{exp}.)\pm 0.006(\mathrm{th}.),`$ (50)
$`\alpha _s(M_Z^2;\overline{\mathrm{MS}})`$ $`=`$ $`0.1167_{0.0118}^{+0.0128}(\mathrm{exp}.)\pm 0.0008(\mathrm{th}.).`$ (51)
We see that the predictions of the applied method suggest that the experimental data on the GLS should be refined significantly in order to increase the predictive power for the QCD coupling parameter.

Another observation, evident from Figs. 25, is that the nonperturbative massless contributions (NP) to $`\mathrm{\Delta }(Q^2)`$ are very significant, and negative. They have their origin, as explained in the previous Section, in the first IR renormalon singularity at $`b=1`$, or equivalently, they correspond approximately to the $`d=2`$ massless power corrections ($`1/Q^2`$). For $`Q^2=2.00,3.16\mathrm{GeV}^2`$, they lead to about $`11\%`$, $`8\%`$, decrease of the value of $`\mathrm{\Delta }(Q^2)`$, respectively. This is to be contrasted with the heavy ($`c`$–)quark contribution\*†\*†\*† The latter is, to a large degree, a $`d=2`$ massive power correction $`m_c^2/Q^2`$, see Subsec. III B. which is positive and leads to only about $`3.6\%`$, $`3.1\%`$ increase of $`\mathrm{\Delta }(Q^2)`$, respectively. If the $`c`$–quark contribution were not included, the central predicted values in (48)–(51) would change to $`\alpha _s(Q^2)=0.367`$, $`0.316`$ \[$`\alpha _s(M_Z^2)=0.1183`$, $`0.1180`$\] for $`Q^2=2.00`$, $`3.16\mathrm{GeV}^2`$, respectively.

The small negative “light–by–light” contribution was separated from our resummation and then added as the term $`\mathrm{\Delta }_{l.l.}(Q^2)w_2^{(0)}(l.l.)a^3(\mu ^2;c_2,\mathrm{})`$. The “light–by–light” part of $`\mathrm{\Delta }_{l.l.}(Q^2)`$, by the special topology of the Feynman diagrams representing it, is a quasiobservable in the sense that it is RScl– and RSch–invariant. Thus, the coefficient $`w_2^{(0)}(l.l.)`$ (27) is the leading coefficient of that quasiobservable and is therefore unchanged under the changes of the RScl and RSch. The “light–by–light” part decreases $`\mathrm{\Delta }(Q^2)`$ by only about $`0.4\%`$ and $`0.3\%`$, for $`Q^2=2.00,3.16\mathrm{GeV}^2`$, respectively.

## IV Comparison with other approaches

One may ask how crucial is the introduction of the conformal transformation (18) for obtaining the numerical predictions (48)–(51). If we repeat the same analysis, but this time without the conformal transformation, and keeping $`\xi ^2=3.3`$, we obtain

$`Q^2=2.00\mathrm{GeV}^2:\alpha _s(Q^2;\overline{\mathrm{MS}})`$ $`=`$ $`0.346_{0.092}^{+\mathrm{?}}(\mathrm{exp}.),\alpha _\mathrm{s}(\mathrm{M}_\mathrm{Z}^2;\overline{\mathrm{MS}})=0.1163_{0.0113}^{+\mathrm{?}}(\mathrm{exp}.);`$ (52)
$`Q^2=3.16\mathrm{GeV}^2:\alpha _s(Q^2;\overline{\mathrm{MS}})`$ $`=`$ $`0.304_{0.073}^{+0.125}(\mathrm{exp}.)),\alpha _\mathrm{s}(\mathrm{M}_\mathrm{Z}^2;\overline{\mathrm{MS}})=0.1165^{+0.0124}_{0.0117}(\mathrm{exp}.).`$ (53)
These results are very close to the results (48)–(51). Thus, we see that the introduction of the conformal transformation (18), which had the task of reducing the influence of the UV and the nonleading IR renormalons, does not influence significantly the predictions. Therefore, we can conclude that these renormalon singularities (at $`b=1,2,\mathrm{}`$ and $`b=2,3,\mathrm{}`$) are in GLS numerically much less important than the leading IR renormalon singularity (at $`b=1`$), even when no conformal transformation is introduced.

We can ask how our predictions compare with those of other, alternative, OPE based methods which, in contrast with the method applied here, do not take into account explicitly the structure of the leading IR renormalon singularity in the Borel plane.

One such an alternative method is the (TPS) PMS optimization of the perturbative contribution, which fixes the RScl and RSch in the TPS in a judicious manner . Resummations of the GLS sum rule based on this method were theoretically and numerically investigated in 1992 by the authors of Ref. . They were confronting the TPS results with the measured values, paying particular attention to the RScl– and RSch–dependence of the predicted values of $`\alpha _s(M_Z^2;\overline{\mathrm{MS}})`$. For the nonperturbative massless (twist–four) $`d=2`$ contribution, they employed the positive value as obtained in Ref.  \[$`\mathrm{\Delta }_{\mathrm{NP}}(Q^2)0.1\mathrm{GeV}^2/Q^2`$\]. Further, the authors of Ref.  accounted for the quark mass threshold effects ($`n_f^{\mathrm{eff}.}`$) by introducing a judiciously chosen weighted average of $`n_f=3,4,5`$.\*‡\*‡\*‡ This is different from our approach, where we separately added the contributions (37)–(42) of the heavy ($`c`$–)quark as corrections to the massless $`n_f=3`$ GLS sum rule, as recently suggested in Ref. . Furthermore, they used the GLS measured values available at that time \[$`\mathrm{\Delta }(Q^2=3\mathrm{G}\mathrm{e}\mathrm{V}^2)=0.167\pm 0.027`$\] which differ significantly from the presently available values (45). Their central value prediction was $`\alpha _s(M_Z^2)=0.115`$, which is lower than our central value predictions (49),(51).

The CCFR group  carried out a numerical analysis similar to that of the authors of Ref. , but with the newer, lower, experimental data (44)–(45) for $`\mathrm{\Delta }(Q^2)`$, and using in the TPS–part $`\overline{\mathrm{MS}}`$ RSch (and RScl $`\mu ^2=Q^2`$). Their central value is $`\alpha _s(3\mathrm{GeV}^2)=0.278`$  The values of $`\alpha _s(3\mathrm{GeV}^2)`$ obtained from our central values of Eqs. (48) and (50) are $`0.309`$ and $`0.310`$, respectively. and $`\alpha _s(M_Z^2)=0.114`$, thus slightly lower than that of Ref. , and significantly lower than our central value predictions (49) and (51). The principal reason for this difference shall be discussed in the following section. Further, if they included in their method the $`\mathrm{N}^3\mathrm{LO}`$ term in the TPS, with $`w_3^{(0)}`$ as estimated in Eq. (31), the predicted value of $`\alpha _s(M_Z^2)`$ would decrease by about $`0.002`$.

Furthermore, the CCFR group mentioned that their central value increases to $`\alpha _s(3\mathrm{GeV}^2)0.305`$  and $`\alpha _s(M_Z^2)=0.118`$  We note that the RGE evolution of $`\alpha _s`$ from $`Q^2`$ to $`M_Z^2`$ gives in our approach different results: $`\alpha _s(3\mathrm{GeV}^2)=0.278(0.305)`$ gives $`\alpha _s(M_Z^2)=0.1124(0.1162)`$ when using the three–loop or four–loop TPS $`\beta `$–function, $`0.1123(0.1161)`$ when using the (four–loop) $`[2/3]`$ PA $`\beta `$–function – we use the corresponding two–loop and three–loop matching conditions for the flavor thresholds ; other details given in Ref. . The CCFR Collaboration apparently uses an approximate three–loop RGE evolution \[a truncated expansion in inverse powers of $`\mathrm{ln}(\mu ^2/\mathrm{\Lambda }_{\overline{\mathrm{MS}}})`$\] and different matching conditions for the flavor thresholds, possibly of Ref. , giving them the values $`0.114(0.118)`$.  when they set the twist–four ($`d=2`$) contribution approximately equal to zero. Such higher–twist values for the GLS sum rule are suggested by the calculation by the authors of Ref.  based on an IR renormalon model with dispersive approach of Ref. , and also by the calculation by the authors of Ref.  based on the bag model. It looks reasonable that the result with the aforementioned IR renormalon model method gives prediction rather close to our prediction, since our method also accounts for the IR renormalon contribution, although in a different manner. However, the calculations in Refs.  apparently do not give us a clear handle on how to treat the perturbative part, i.e., whether to take it as a LO, NLO, or NNLO TPS, or in any other form. This is in contrast with our method, where the perturbative and nonperturbative parts are clearly connected with each other. We discuss this aspect in more detail in the next section.

At this point, we would like to point out that the PMS method has a signal casting doubts on its applicability in the discussed GLS cases– namely, the PMS RScl is very low in this case: $`\mu _{\mathrm{PMS}}^20.203Q^2`$. For $`Q^2=2.00,3.16\mathrm{GeV}^2`$, this implies the scales $`\mu _{\mathrm{PMS}}0.64,0.80\mathrm{GeV}`$, respectively, which may be too low for the application of perturbative approaches such as PMS. The same problem appears when applying the effective charge (ECH) method  to these GLS cases.

The present world average for the QCD coupling parameter is $`\alpha _s^{\overline{\mathrm{MS}}}(M_\text{z}^2)=0.1173\pm 0.0020`$ by Ref.  and $`0.1184\pm 0.0031`$ by Ref. . Comparing this with our predictions (49) and (51), we see that the method applied in the present paper gives us the central values which agree well with the present world average. We wish to point out that this agreement suggests that the method applied in the present paper for the nonperturbative massless correction to $`\mathrm{\Delta }(Q^2)`$ is at least consistent with the experiments. If this correction were zero, or had the opposite sign, the obtained central prediction for $`\alpha _s(M_Z^2)`$ would be at the lower edge or even outside the interval of the present world average. These considerations do not necessarily imply, but indicate, that the applied method gives the correct nonperturbative contributions. For a more definite statement in this respect, the experimental uncertainties in the GLS sum rule would have to be reduced significantly.

## V Borel resummation versus OPE approach

In the discussions so far, we considered the amplitudes of Borel resummation or OPE approach only at fixed values of $`Q^2`$. Here we mean, for convenience, by the OPE approach the usual perturbative expansion plus a power suppressed term representing the twist–four contribution. In this section we consider them over a continuous range of $`Q^2`$. This slight change of view will reveal the characteristic features of the two approaches, and enable us to better understand the cause of the significant difference in the extracted strong coupling constants seen in the previous section.

We first note the remarkable stability of the Borel resummed amplitudes over the order of perturbation involved in their calculation. In Fig. 6 (a)-(b), we plot, over the interval $`1<Q^2<10,`$ in $`\mathrm{GeV}^2`$, the real part of the Borel resummed amplitudes for $`\mathrm{\Delta }_\mathrm{P}(Q^2)`$ using NLO, NNLO, and N3LO perturbations and the corresponding amplitudes of the ordinary TPS’s in the OPE approach. In all of the plots in Fig. 6 the N3LO QCD $`\beta `$-function was used in the running of the strong coupling. We also take the RG scale at $`\mu ^2=Q^2`$, since the RG scale dependence is sufficiently small for our present discussion. The aforementioned stability of the Borel resummed amplitudes becomes clear when the two figures are compared. Note the variation in the Borel resummed amplitudes is very small, whereas the TPS amplitudes have significant order dependence. While this stability is not completely unexpected, because the leading renormalon singularity is effectively softened by the use of the function $`R(b)`$ in the Borel integration, the degree of the stability is still remarkable. This suggests that the renormalon–induced asymptotic behavior of the perturbative coefficients sets in quite early in perturbation, and that the use of $`R(b)`$ and conformal mapping in Borel resummation is very efficient in handling the renormalon singularity. We also note in passing that the stability of the Borel resummed amplitudes for $`\mathrm{\Delta }_\mathrm{P}(Q^2)+\mathrm{\Delta }_{\mathrm{NP}}(Q^2)`$, which are not shown in the figure, is comparable to that of $`\mathrm{\Delta }_\mathrm{P}(Q^2)`$.

In the previous section we have seen there is a significant difference between the Borel resummation and the OPE approach in the prediction of the strong coupling constant. While there are obvious differences in the two approaches, it was not clear what aspect of the Borel resummation is primarily responsible for the difference. Is it because of the perturbative part $`\mathrm{\Delta }_\mathrm{P}(Q^2)`$ or because of our specific implementation of the nonperturbative part $`\mathrm{\Delta }_{\mathrm{NP}}(Q^2)`$, or both?

To answer this question we plot in Fig. 6 (c) the Borel resummed $`\mathrm{\Delta }_\mathrm{P}(Q^2)+\mathrm{\Delta }_{\mathrm{NP}}(Q^2)`$ against the OPE amplitudes $`(\mathrm{NLO}\mathrm{TPS})+0.1/Q^2`$, whose power term representing the twist–four contribution is from the sum rule calculation ,\*∥\*∥\*∥We note, however, there are some variations in the estimate of the twist–four contribution. The sum rule calculation of Ref.  predicts the ($`d=2`$) power term $`0.16\mathrm{GeV}^2/Q^2`$. and $`(\mathrm{NNLO}\mathrm{TPS})+0.02/Q^2`$. The small power term in the latter was chosen for the amplitude to match the NLO OPE amplitude at large values of $`Q^2`$ in the plots. In the figure we first notice that the NLO and NNLO OPE amplitudes with a large difference in twist–four contribution agree reasonably well over the whole range of $`Q^2`$ considered. This implies that the higher–twist term in the NLO amplitude can be largely accounted for by the NNLO perturbative term, which is in qualitative agreement with the observation in Refs. . This also shows that the use of the NLO sum rule calculation of the higher–twist contribution with a TPS of different order, which is not an uncommon practice, can be dangerous. Higher–twist contributions calculated at a given order of the leading perturbative contribution should never be used with a TPS of different order.

On the other hand, the Borel resummed amplitude is in a reasonably good agreement with the OPE amplitudes at large $`Q^2(>4\mathrm{GeV}^2`$), but deviates significantly at small momenta. Obviously, this deviation at small momenta explains the difference in the prediction of the strong coupling. Before we answer the origin of this deviation, we note that the good agreement of the two approaches at large momenta is a nontrivial result. Even though the amplitudes from the two approaches should agree at very high momenta (or at small couplings), since they have the same low order perturbations, the degree of agreement seen here is unlikely to be a random consequence. For instance, if our $`\mathrm{\Delta }_{\mathrm{NP}}(Q^2)`$ had the wrong sign or were zero, then we would see a significant difference at high momenta (see the dotted plot for the wrong sign case). Thus this good agreement of our Borel resummed amplitude with the OPE amplitudes at high momenta may be regarded as a partial support for our prescription of the nonperturbative part.

Now back to the question of what causes the deviation at low momenta. It is not difficult to see $`\mathrm{\Delta }_\mathrm{P}(Q^2)`$ must be responsible for the deviation. The reason is as follows. Since the nonperturbative part $`\mathrm{\Delta }_{\mathrm{NP}}(Q^2)`$ essentially behaves like a power suppressed term,\***\***\*** $`\mathrm{\Delta }_{\mathrm{NP}}(Q^2)\alpha _s(Q^2)^{\gamma _2/\beta _0}/Q^2=\alpha _s(Q^2)^{32/81}/Q^2`$, in accordance with Eqs. (14)–(15) and the OPE calculation of Ref. . and $`\mathrm{\Delta }_\mathrm{P}(Q^2)+\mathrm{\Delta }_{\mathrm{NP}}(Q^2)`$ is in agreement with the OPE amplitudes at large momenta, so should it be at low momenta, too, were $`\mathrm{\Delta }_\mathrm{P}(Q^2)`$ to behave like an OPE amplitude. Thus the primary cause of the difference in the predicted strong coupling constant must be that the Borel resummed $`\mathrm{\Delta }_\mathrm{P}(Q^2)`$ at low momenta cannot be parametrized in the form of an OPE amplitude. We can see this explicitly by looking at the plots in Fig. 6 (d) where, as an example, the NNLO $`\mathrm{Re}[\mathrm{\Delta }_\mathrm{P}(Q^2)]`$ is plotted against an OPE amplitude $`(\mathrm{NNLO}\mathrm{TPS})+0.16/Q^2`$. The power suppressed term in the latter was fixed so that the two amplitudes match at high momenta in the plots. Clearly, they deviate significantly at low momenta, with the Borel resummed growing more slowly than the OPE amplitude as the momentum is decreased.

This suggests the OPE amplitude tends to overestimate at low momenta. We can easily see that this tendency arises from the bad functional form of the polynomial Borel transform (TPS of (8)) around the renormalon singularity at $`b=1`$. In Fig. 7 (a) the N3LO Borel transforms in the two approaches are plotted in the variable $`b`$. For $`b>1`$ the Borel transform in our approach defined through $`R(b)`$ in (17) is complex, and its real part is plotted. It is obvious that the OPE Borel transform is badly broken around and beyond the renormalon singularity. When the coupling is small this is not a serious problem because the dominant contribution to the Borel integral (7) comes from the region close to the origin. However, as the coupling becomes larger the relevant integration region extends to the renormalon singularity, and beyond, and as we see in the plots, the OPE Borel transform can grossly overestimate at large couplings. The amplitudes obtained from these Borel transforms are plotted in Fig. 7 (b). As expected, at small momenta (large coupling) the OPE amplitude is larger than the Borel resummed. Note, on the other hand, at high momenta it is smaller than the latter. This is because the OPE Borel transform in the region $`0<b<1`$, from which the dominant contribution comes at small couplings, is smaller than the other one, which is a characteristic feature rendered automatically by the correct implementation of the renormalon singularity in the latter. This difference between the Borel resummed and the OPE amplitude at small couplings may be regarded as the resummation of the unaccounted higher order terms in the same sign asymptotic series.

That the Borel resummed $`\mathrm{\Delta }_\mathrm{P}(Q^2)`$ cannot be approximated by an OPE amplitude of a TPS plus a power suppressed term representing the renormalon effect may appear contradictory to the common opinion which states otherwise. The latter opinion, which is based on the factorially growing large order behavior and the running coupling, would be true in a sense, provided the strong coupling were sufficiently small, and a TPS of large order ($`1/\alpha _s`$) was used. In reality, however, the strong couplings at the low momenta we consider are not so small, and there is no guarantee that the Borel resummed $`\mathrm{\Delta }_\mathrm{P}(Q^2)`$ at those momenta can be parameterized as an OPE amplitude. The example here clearly shows that cannot be generally true. A Borel resummed amplitude can have a much more complex functional behavior than the sum of a TPS and a power term intended to account for the renormalon effect. This consideration suggests that several existing analyses of low energy QCD observables based on the OPE approach should be reexamined, since the issues raised here are likely to be relevant there, too.

To sum up this section, we have made two observations concerning the Borel resummation and the OPE approach. First, our method of calculation of the perturbative plus nonperturbative contribution in Borel resummation is consistent at larger $`Q^2>4\mathrm{GeV}^2`$ with the OPE approach using QCD sum rule calculation, and secondly, at low energies the OPE amplitude tends to overestimate, and the Borel resummed amplitude with a proper incorporation of the leading renormalon cannot be approximated by an OPE amplitude of a TPS with a power suppressed term. It is the second observation that directly accounts for the differences in the extracted strong coupling constants from the two approaches.

## VI Conclusions

We performed a resummation of the Gross–Llewellyn Smith (GLS) sum rule by fully accounting for the correct known form of the leading infrared (IR) renormalon singularity at $`b=1`$ in the Borel plane. As one direct consequence of this singularity, the resummed “perturbative” part $`\mathrm{\Delta }_\mathrm{P}(\alpha _s(Q))`$ of the GLS has a branch cut of the form $`(\alpha _s(Q^2))^\nu \mathrm{exp}(\pi /\beta _0\alpha _s(Q))[1+𝒪(\alpha _s)]`$, i.e., a twist–four term $`(1)^\nu (1/Q^2)\alpha _s^{\gamma _2/\beta _0}[1+𝒪(\alpha _s)]`$ with the branch cut discontinuity factor $`\mathrm{exp}(\pm i\pi \nu )`$. Here, $`\gamma _2=8/9`$ is the known one–loop coefficient of the anomalous dimension of the corresponding twist–four $`d=2`$ operator appearing in the OPE. A crucial element of the analysis was the fixing of the “nonperturbative” part as the negative of the aforementioned branch cut term, thus making the resummed (“perturbative” plus “nonperturbative”) GLS sum rule manifestly real. This procedure is free from the known ambiguity of separation of the “perturbative” and “nonperturbative” parts. In the Borel resummation of the “perturbative” part, we further employed a conformal transformation to minimize the numerical influence of other renormalon singularities. All this allowed us to perform the resummation of the massless part of the GLS sum rule, i.e., of the contributions of the three light quark flavors. The contributions of the heavy ($`c`$–)quark were added separately, as were the target nuclear correction contributions and the light–by–light contributions. These three types of contributions turned out to be small, in contrast to the “nonperturbative” contributions which turned out to be significant. The calculations were performed in the $`\overline{\mathrm{MS}}`$ renormalization scheme, and the renormalization scale $`\mu ^2`$ ($`Q^2`$) was taken in the region of the smallest $`\mu ^2`$–sensitivity of our results.

We then confronted the resummed expressions with the Fermilab CCFR Collaboration data  for the GLS (at $`Q^2=2,3.16\mathrm{GeV}^2`$) which already include the target mass corrections. Our central value prediction for the QCD coupling parameter, corresponding to the central GLS values of the CCFR, is $`\alpha _s(M_Z)0.117`$ \[see Eqs. (48) and (51)\], in good agreement with the present world average. This is different from the central value predictions of previous analyses of the GLS sum rule by the CCFR Collaboration  ($`\alpha _s(M_Z)0.114`$) and by the authors of  ($`\alpha _s(M_Z)0.115`$) which are below or at the lower edge of the world average. We have seen that our approach to the calculation of the nonperturbative contribution is consistent with the OPE approach, and the main reason for the difference between the two approaches is that at small $`Q^2<4\mathrm{GeV}^2`$, the OPE approach tends to overestimate and the Borel resummed perturbative contribution cannot be approximated by an OPE amplitude.

The GLS sum rule, at the low gauge boson transfer momenta $`Q^2=2`$$`4\mathrm{GeV}^2`$, is a very important quantity to measure, because it has apparently a strong nonperturbative component, stronger than in some other low–energy QCD observables such as the semihadronic tau decay rate. The more precise experimental values of the GLS sum rule would help determining the higher–twist contributions more accurately.

###### Acknowledgements.

The work of C.C. was supported by DGIP (UTFSM). The work of G.C. was supported by the FONDECYT (Chile) Grant No. 1010094. The work of T.L. and K.S.J. was supported in part by the BK21 Core Program.

## Appendix A Branch cut singularity and the nonperturbative part

In this Appendix, we show explicitly formulas leading to Eqs. (13)–(15), using the identity (7) and the leading IR renormalon singularity structure (9). The latter structure around $`b=1`$ can be rewritten more explicitly as

$`\stackrel{~}{W}_0(b)`$ $`=`$ $`{\displaystyle \frac{C}{\mathrm{\Gamma }(\nu )}}\beta _0^{1+\nu }(1b)^{1\nu }\left[1+\kappa _1(1b)+\kappa _2(1b)^2+\mathrm{}\right]`$ (A.2)
$`+(\mathrm{Analytic}\mathrm{part}).`$
The leading part is known ($`\nu `$), while the coefficients $`\kappa _j`$ ($`j1`$) of the subleading parts are not yet known. Inserting the above expansion into the Borel integration formula (7) and performing the change of the integration variable $`b=1+\beta _0at`$, we obtain

$$\mathrm{Im}\mathrm{\Delta }_\mathrm{P}(a(Q)\pm \mathrm{i}ϵ)=\frac{C}{\mathrm{\Gamma }(\nu )}e^{1/\beta _0a(Q)}a(Q)^\nu f_\pm (a(Q)),$$
(A.3)
where

$`f_\pm (a)`$ $`=`$ $`\mathrm{Im}{\displaystyle _{0\pm \mathrm{i}ϵ}^{\mathrm{}\pm \mathrm{i}ϵ}}𝑑te^t(t)^{1\nu }\left[1+\kappa _1(\beta _0a)(t)+\kappa _2(\beta _0a)^2(t)^2+\mathrm{}\right]`$ (A.4)
$`=`$ $`\mathrm{sin}(\pi \nu ){\displaystyle _0^{\mathrm{}}}𝑑te^tt^{1\nu }\left[1+\kappa _1(\beta _0a)(t)+\kappa _2(\beta _0a)^2(t)^2+\mathrm{}\right]`$ (A.5)
$`=`$ $`\mathrm{\Gamma }(\nu )\mathrm{sin}(\pi \nu )\left[1+\kappa _1(\beta _0a)\nu +\kappa _2(\beta _0a)^2\nu (\nu 1)+\mathrm{}\right]`$ (A.6)
Now requiring the imaginary parts of $`\mathrm{\Delta }_\mathrm{P}`$ at positive $`a(Q)`$ to match those in (A.3) $`\mathrm{\Delta }_\mathrm{P}`$ for general complex $`a(Q)`$ can be written as

$`\mathrm{\Delta }_\mathrm{P}(a(Q))`$ $`=`$ $`Ce^{1/\beta _0a(Q)}\{(a(Q))^\nu +\kappa _1(\nu )\beta _0(a(Q))^{\nu +1}`$ (A.9)
$`+\kappa _2(\nu )(\nu +1)\beta _0^2(a(Q))^{\nu +2}+\mathrm{}\}`$
$`+\stackrel{~}{\mathrm{\Delta }}_\mathrm{P}(a(Q)),`$
where $`\stackrel{~}{\mathrm{\Delta }}_\mathrm{P}(a)`$ is the part with no singularities (no cuts) for $`a>0`$. \*††\*††\*†† The absolute value of the singular term in (A.9) can be rewritten as $`(\mathrm{\Lambda }^2/Q^2)a(Q)^{+\gamma _2/\beta _0}[1+𝒪(a)]`$, where $`\mathrm{\Lambda }`$ is a $`Q`$–independent scale, and $`\gamma _2`$ is the one–loop coefficient of the anomalous dimension of the twist–four ($`d=2`$) operator $`O(Q)\mathrm{\Lambda }^2a(Q)^{+\gamma _2/\beta _0}`$ appearing in the OPE (3). According to Ref. , the nonperturbative part $`\mathrm{\Delta }_{\mathrm{NP}}(a(Q)\pm \mathrm{i}ϵ)`$ must cancel the imaginary part of $`\mathrm{\Delta }_\mathrm{P}(a(Q)\pm \mathrm{i}ϵ)`$, and $`\mathrm{\Delta }_{\mathrm{NP}}`$ was chosen to have the simplest, presumably the most natural, form – i.e., just the negative of the branch cut term of (A.9)

$`\mathrm{\Delta }_{\mathrm{NP}}(a(Q))`$ $`=`$ $`+Ce^{1/\beta _0a(Q)}\{(a(Q))^\nu +\kappa _1(\nu )\beta _0(a(Q))^{\nu +1}`$ (A.11)
$`+\kappa _2(\nu )(\nu +1)\beta _0^2(a(Q))^{\nu +2}+\mathrm{}\}.`$
Further, since $`(a\mathrm{i}ϵ)^{\nu +n}=a^{\nu +n}\mathrm{exp}[\pm \mathrm{i}(\nu n)\pi ]`$, expressions (A.3) and (A.11) immediately relate $`\mathrm{Re}\mathrm{\Delta }_{\mathrm{NP}}(a\pm \mathrm{i}ϵ)`$ with the (calculable) $`\mathrm{Im}\mathrm{\Delta }_\mathrm{P}(a\pm \mathrm{i}ϵ)`$

$$\mathrm{Re}\mathrm{\Delta }_{\mathrm{NP}}(a\pm \mathrm{i}ϵ)=\mathrm{cot}(\nu \pi )\mathrm{Im}\mathrm{\Delta }_\mathrm{P}(a\pm \mathrm{i}ϵ),$$
(A.12)
giving the result (15).