format: text

text:
# essai

## 1 Introduction

The standard collinear factorization in deep inelastic scattering (DIS) at twist two is based on separating hard physics associated with scales of order of the momentum transfer $`Q^2`$ into coefficient functions, and softer physics associated with the structure of the target into operator matrix elements which have the interpretation of parton distribution functions.

Considering the kinematic limit $`x1`$ (where Bjorken $`x`$ is defined by $`xq^2/2pq=Q^2/2pq`$), an additional, well-separated scale emerges, $`W^2=Q^2(1x)/x`$ which represents the invariant mass of the hadronic system. This is the characteristic scale for dynamics of the jet in the final state. When $`W^2`$ is much smaller than $`Q^2`$, perturbative corrections (Sudakov double logs) related to soft and collinear radiation which form the jet dominate the coefficient functions. These corrections become so large that the perturbative expansion breaks down. Consequently, they must be resummed to all orders. In standard analysis of structure function data the difficult large-$`x`$ region is usually avoided simply by putting a lower cut on $`W^2`$. For example, in the MRST analysis  a cut is put such that $`W^2>12.5`$ GeV2.

To describe the structure functions at large Bjorken $`x`$, it becomes natural to further factorize the process and treat the jet separately. Such factorization can indeed be accomplished , facilitating the resummation of the logarithmic corrections. A schematic picture of large-$`x`$ factorization is shown in Fig. 1.

It is well known that higher-twist corrections are also enhanced at large $`x`$: they appear as powers of $`\mathrm{\Lambda }^2/W^2`$ rather than $`\mathrm{\Lambda }^2/Q^2`$. At sufficiently small $`W^2`$ the Operator Product Expansion (OPE) itself breaks down and needs to be resummed. Given the complexity of this expansion , the only hope this could ever be achieved is if a small subset of matrix elements dominates at each twist, simplifying greatly the expansion. There are two ways in which the leading twist simplifies at large $`x`$. Firstly, the valence quarks dominate over the gluons and the sea. Secondly, as explained above (Fig. 1) the dynamics of the outgoing jet, associated with the scale $`W^2`$, decouples from the rest of the process. The question then arises whether any of this holds beyond the leading twist.

Recently, a first theoretical analysis of the large-$`x`$ limit of twist four has been performed , suggesting that indeed, this kind of simplification occurs to all orders in the twist expansion, and thus its resummation at large $`x`$ may be possible.

Ref.  first establishes, in the concrete case of $`F_2`$, the cancellation of the infrared renormalon ambiguity at twist two with the quadratically divergent contribution of twist-four operators, associated with their mixing with twist two. Renormalons reflect the ambiguity in separating contributions of different twists. Independently of how this separation is implemented, within twist four there is a twist-two like ingredient, which is proportional to the twist two matrix elements. The question that arises is what is the significance of this twist-two like component compared with other “genuine” higher-twist contributions. Renormalon-based models for power corrections  assume its dominance whereas other approaches may neglect it altogether.

The conjecture of  is that at large $`x`$ the twist-two like ingredient within the higher twist indeed dominates. Analysis of the leading-order coefficient functions of twist four for $`F_2`$ and $`F_L`$ provides a physical picture which is consistent with this conjecture. The final state which dominates the twist-four contribution at large $`x`$ is the same as the twist-two final state: in the case of $`F_2`$ this state contains a single energetic quark which initiates the jet. Moreover, the analysis of  also reveals that the multi-parton initial states at twist four, which distinguish it from twist two, are dominated at large $`x`$ by particular configurations which are twist-two like: a single quark in the initial state carries most of the momentum, whereas additional gluons carry small momentum fractions. Consequently, the corresponding twist-four matrix elements effectively depend only on the light-cone separation between the quark fields, similarly to twist two.

The conclusion is that factorization, similar to the one described by Fig. 1, may hold beyond the perturbative level. The essential difference is that beyond twist two, gluons of virtuality of the order $`W^2`$ are exchanged between the jet and the remnants of the target. Their effect, to all orders in the OPE, can be taken into account through a single “shape function” of $`N/Q^2`$, which multiplies the twist-two contribution in moment space. This way, the conjecture of ref.  translates into a concrete understanding of how higher-twist effects should be parametrized at large $`x`$.

The need to resum power-corrections to all orders as a consequence of kinematic thresholds exists in a wide range of hard processes, including, for example, Drell-Yan pair production and event-shape distributions. In spite of the different nature of these processes, the analogy is quite useful. In particular, the appearance of a new non-perturbative distribution , the shape function, near the kinematic threshold, is common to all these processes, and so is the relation between the resummed perturbative distribution and the properties of this non-perturbative function .

It should be stressed that the proposed non-perturbative factorization formula for DIS structure functions has not been derived from the OPE. It requires that the Sudakov resummed coefficient function of higher-twist contributions to $`F_2`$ would be the same as that of twist two. Moreover, the large-$`N`$ behaviour of higher-twist anomalous dimensions should coincide, asymptotically, with that of twist two. This highly non-trivial structure is yet to be verified by explicit calculations.

Having assumed the dominance of these higher-twist contributions which mix under renormalization with the leading twist (and are therefore proportional to the latter), renormalon resummation becomes absolutely essential for any higher-twist analysis. Indeed, the separation between the leading and higher twist is ambiguous. But, more importantly, prior to introducing any parametrization of power suppressed contributions one must be sure that any large perturbative corrections have already been taken into account . This concerns, in particular, running-coupling (or renormalon-related) perturbative corrections which are always parametrically larger than the corresponding ambiguity, and thus, by the previous assumption larger than the higher-twist contribution itself. Failing to take these perturbative corrections into account, the parametrization of the higher twist becomes meaningless. For example, the values of the extracted “higher twist” parameters will strongly depend on the order of truncation of the perturbative expansion, the renormalization scale used, etc.

The unresolved issue of higher twist in DIS structure functions dates far back. It was understood very early on that higher-twist effects are especially important at large $`x`$. The limitations have always been both insufficient data and lack of theoretical understanding of this region. In the recent years there has been much activity in this field, primarily owing to progress in perturbative calculations and resummation. Several groups have performed fits to data , incorporating (or not) soft gluon resummation and parametrizing the higher twist in various ways. What lacks however, is a physical picture that stands behind these parametrizations. In particular, one is bound to ask eventually how the higher-twist parameters are related to operator matrix elements within the OPE, namely which multi-parton correlations are being measured. Ref.  suggests an answer111A previous suggestion with some similarity is that of ref. .. More theoretical and experimental work is required to establish it. Here we make a further step in this direction. We concentrate mainly on perturbative aspects, constraining further the structure of the Sudakov exponent and examining where we stand on the long track to power accuracy.

As mentioned above, large-$`x`$ kinematics singles out very specific radiative corrections which dominate the coefficient function: these are the logarithmically-enhanced terms associated with the evolution of the jet. The jet function can be written as an exponential in moment space, where the exponentiation kernel at a given order can be systematically deduced from the standard loop expansion of the same order. To a first approximation (leading-log accuracy) the kernel simply corresponds to a single gluon emission calculation, whereas multiple emission is accounted for by the exponentiation (or solving the evolution equation with that kernel). Refining the calculation of the kernel to next-to-leading logarithmic accuracy amounts to taking into account the leading contribution associated with the running of the coupling. A natural procedure to further improve the approximation to the kernel is a renormalon calculation with a single dressed gluon . In this way some running-coupling effects are taken into account to any logarithmic accuracy, making the exponent renormalization-scale invariant. Consequently the related power-corrections can be systematically parametrized, in accordance with .

The essential difference between dressed gluon exponentiation (DGE), and other approaches to Sudakov resummation, is that it takes into account some all-order information on the exponentiation kernel itself. However, since the renormalon calculation is restricted to the large-$`N_f`$ (or, equivalently, large $`\beta _0`$) limit, it is necessary to combine it with what is known about the first few orders in the kernel from fixed-order calculations. One of the goals of this paper is to combine the large-$`\beta _0`$ result  for the Sudakov exponent of $`F_2`$ with the state-of-the-art knowledge of the $`\overline{\mathrm{MS}}`$ anomalous dimension and coefficient function. The latter practically222The non-Abelian contribution to the $`\overline{\mathrm{MS}}`$ anomalous dimension is still not available analytically as a function of $`N`$. However, there is a reliable estimate  of the leading $`\mathrm{ln}N`$ contribution based on the computation of specific moments . allows a complete  next-to-next-to-leading logarithmic (NNLL) accuracy calculation, so together with the all-order large-$`N_f`$ result, it seems that the Sudakov exponent as a whole is well constrained. Nevertheless, we will see that without making further assumptions concerning the structure of the Sudakov exponent, there still is a significant uncertainty concerning the magnitude of subleading logs. This has important consequences for phenomenology.

In this paper we apply, for the first time, DGE with the corresponding parametrization of higher-twist corrections to DIS data. We restrict our attention here to the large-$`x`$ region, by analysing moments $`N5`$. The moment space analysis is advantageous from a theoretical point of view in several respects. First of all, it facilitates the implementation of target-mass corrections through the use of Nachtmann weights. In addition, it simplifies significantly the resummation formulae for the coefficient function as well as the higher-twist corrections.

The paper is organized as follows. The theoretical part (section 2) begins by recalling the non-perturbative factorization formula of . We then return to discuss factorization on the perturbative level and its consequences for resummation (section 2.1). We emphasize in particular the fact that Sudakov logs in the physical evolution kernel are exclusively related to the jet function and not to the soft function. This translates into a more predictive formulation of Sudakov resummation than that commonly presented in the literature. This discussion is followed by an all-order analysis of twist-two at large $`N`$, first (section 2.2) for the physical kernel, and then (section 2.3) for the coefficient function. Sections 2.4 and 2.5 are devoted to combine the all-order large-$`N_f`$ information with the available coefficients in QCD and discuss in some detail the uncertainty involved in such a procedure. Section 3 is devoted to data analysis. Matching the DGE exponent into the known  next-to-next-to-leading order (NNLO) result we obtain an improved perturbative prediction for the scaling violation at large $`N`$, which we use as a baseline for the study of power corrections. The experimental moments $`N=4`$$`11`$ of $`F_2`$ for the proton are calculated based on SLAC and BCDMS data. Each of the moments is fitted separately, to extract the valence quark distribution, $`\alpha _s`$ and the magnitude of the higher twist. We use the stability of the extracted value of $`\alpha _s`$ as a function of $`N`$ as a measure of the quality of our theoretical description of $`F_2`$, and compare the results to the standard NNLO analysis. In section 4 we summarize our conclusions.

## 2 Factorization and the Sudakov exponent by DGE

### 2.1 Factorization

The non-perturbative factorization formula  for $`F_2`$ at large $`N`$ takes the form:

$`F_2^N(Q^2){\displaystyle _0^1}𝑑xx^{N2}F_2(x,Q^2)=H\left(Q^2\right)J_N(Q^2;\mu _F^2)q_N(\mu _F^2)J^{\text{NP}}\left(N\mathrm{\Lambda }^2/Q^2\right).`$ (1)
We will assume that it holds up to corrections of order $`1/N`$. Here, the first three factors correspond to the standard, perturbative large-$`x`$ factorization : $`H\left(Q^2\right)`$ is the hard part of the coefficient function depending on the momentum transfer $`Q^2`$, $`J_N(Q^2;\mu _F^2)`$ is the Sudakov-resumed jet function depending primarily on the invariant mass of the jet $`Q^2/N`$, and $`q_N(\mu _F^2)`$ is the twist-two quark matrix element. The last factor $`J^{\text{NP}}\left(N\mathrm{\Lambda }^2/Q^2\right)=1+\kappa _1\frac{N\mathrm{\Lambda }^2}{Q^2}+\mathrm{}`$, resums the dominant higher twist corrections to all orders .

Some important clarifications are due concerning the perturbative333The reader in referred to  for discussion on power corrections. factorization formula. First, the coefficients in the hard function $`H\left(Q^2\right)`$ are finite at large $`N`$, so this function does not play any rôle in the large-$`N`$ limit. The jet and the quark distribution are defined as follows:

$$J(Q^2(1\xi )/\xi ;\mu _F^2)=\frac{1}{2\pi }\mathrm{Im}d^4\stackrel{~}{y}\mathrm{e}^{i(q+zp)\stackrel{~}{y}}0\left|\mathrm{Tr}\left(\mathrm{\Psi }(\stackrel{~}{y})\overline{\mathrm{\Psi }}(0)p/\right)\right|0,$$
(2)
with $`\xi x/z`$, and

$$q(z;\mu _F^2)=_{\mathrm{}}^{\mathrm{}}\frac{d(py)}{py}\mathrm{e}^{izpy}p\left|\overline{\mathrm{\Psi }}(y)y/\mathrm{\Psi }(0)\right|p(y^2=0)$$
(3)
where the gauge in (2) and (3) is chosen as $`A_{}=0`$ and $`A_+=0`$, respectively. Here the lightcone ‘$`+`$’ direction is defined by the incoming hadron $`p`$, and the lightcone ‘$``$’ direction by the outgoing jet, $`(q+xp)`$. $`J`$ is simply a projection of the quark propagator in the axial gauge $`pA=0`$. At leading order $`\mathrm{\Psi }(\stackrel{~}{y})\overline{\mathrm{\Psi }}(0)`$ is the free quark propagator and $`J=\delta (zx)`$; at higher orders $`J`$ depends on $`(q+zp)^2=Q^2(1\xi )/\xi Q^2(1\xi )`$ but also, through the gauge fixing, on the hard scale $`pq`$. The latter dependence is unrelated to logarithmically enhanced terms and it can be neglected. Corrections on this scale appear through the hard function $`H(Q^2)`$.

In spite of the explicit dependence on the gauge, $`J(Q^2(1\xi )/\xi ;\mu _F^2)`$ and $`q(z;\mu _F^2)`$ are gauge invariant. In case of different gauge choices, an appropriate path-ordered exponential should be introduced. For $`J`$ the path goes along the incoming quark direction ($`p`$), connecting the point $`\stackrel{~}{y}`$ to infinity and then infinity to $`0`$, whereas for $`q`$ the path goes along the jet direction between $`y`$ and $`0`$. The different paths are shown by the dashed lines in figure 2. The functions in (1) are defined in moment space, $`J_N(Q^2;\mu _F^2)=_0^1𝑑xx^{N1}J((1x)Q^2;\mu _F^2)`$ and $`q_N(\mu _F^2)=_0^1𝑑xx^{N1}q(x;\mu _F^2)`$, where convolution in the momentum fraction simply becomes a product.

The $`\mu _F^2`$ dependence in (2) and (3) is there to remind us that the definition of these matrix elements requires a prescription (factorization scheme and scale) to deal with infrared and ultraviolet divergence, respectively. Since $`F_2`$ is a physical quantity, this dependence must cancel out in the product.

The standard formulation of large-$`x`$ factorization  includes an additional soft function $`V`$:

$$F_2(x,Q^2)=H\left(Q^2\right)J\left((1x)Q^2\right)V(x;\mu _F^2)q(x;\mu _F^2).$$
(4)
defined  by the path-ordered exponential shown in the lower picture in figure 2, which represents the DIS process in the Eikonal approximation. The proof of factorization  relies on separating a generic Feynman diagram into subprocesses, each depending on a single scale. In inclusive DIS there are two external scales: $`Q^2`$ and $`W^2=Q^2(1x)/xQ^2/N`$. The corresponding subprocesses are $`H(Q^2)`$ and $`J(Q^2(1x))`$. Momentum scales which are parametrically smaller than $`W^2`$, e.g. $`𝒪\left(Q^2(1x)^2\right)`$, are accounted for by the soft subprocess $`V(x,\mu _F^2)`$.

It is important to note that $`V(x,\mu _F^2)`$ does not depend444Dependence of this function on $`Q^2`$ or on $`Q^2/N^2`$ is sometimes introduced through the factorization scale. This, however, can be misleading, and we will avoid it in this formal discussion. on any external scale in the DIS process, and it is therefore possible and useful(!) to eliminate it from the factorization formula . Although not written explicitly in eq. (4), in order to define the soft function one actually performs two factorization procedures with two different scales $`\mu _{F_{1,2}}^2`$: one in the definition of $`J((1x)Q^2;\mu _{F_1}^2)`$ and the other in the definition of $`q(x;\mu _{F_2}^2)`$. This can be intuitively undersood from Fig. 1, where the soft blob is attached to jet and to the quark distribution. In any case some component of the soft blob is included in the jet and some in the quark distribution. Since $`V`$ represents the evolution of the quark distribution as well as that of the jet function with $`\mu _F^2`$, it can be eliminated altogether by choosing the same scale in both. This is demonstrated in figure 3.

The factorization formula is then

$`F_2^N(Q^2)`$ $`=`$ $`H\left(Q^2\right)J_N(Q^2;\mu _{F_1}^2)V_N(\mu _{F_1}^2,\mu _{F_2}^2)q_N(\mu _{F_2}^2)`$ (5)
$`=`$ $`H\left(Q^2\right)J_N(Q^2;\mu _{F_2}^2)q_N(\mu _{F_2}^2).`$
where

$$J_N(Q^2;\mu _{F_2}^2)J_N(Q^2;\mu _{F_1}^2)V_N(\mu _{F_1}^2,\mu _{F_2}^2).$$

In the factorization formula (1) the soft function appears only implicitly through the dependence of $`J_N(Q^2;\mu _F^2)`$ and $`q_N(\mu _F^2)`$ on $`\mu _F^2`$. This function has some remarkable properties which will become relevant in the following. As discussed in ref.  (see section 7 there) and in  its exponent contains just a single $`\mathrm{ln}N`$ at any order in the coupling. It thus follows from the factorization formula that apart from these $`\mu _F^2`$-related single logs, the only source of logarithmically enhanced contributions to the $`F_2`$ coefficient function is in the constraint on the invariant mass of the jet. This is true to all orders in perturbation theory. As we discuss further in the next section, this implies, in particular, that the jet function $`J((1x)Q^2;\mu _F^2)`$ alone determines the Sudakov exponent in the physical kernel $`d\mathrm{ln}F_2^N(Q^2)/d\mathrm{ln}Q^2`$.

We stress that the absence of additional sources of Sudakov logs is strictly related to the inclusive nature of structure functions. Note that, due to different kinematic circumstances, sensitivity to large-angle soft emission which is associated with the scale $`Q^2/N^2`$ does generate Sudakov double logs in other observables. Examples are provided by Drell-Yan pair production and event-shape distributions.

Finally, in the notation of , our statement that the only source of Sudakov logs in the physical kernel is in the jet function means that $`D_n^{\mathrm{DIS}}`$ vanish to all orders. Indeed, recently, it was shown , based on explicit calculations in $`\overline{\mathrm{MS}}`$ that $`D_2^{\mathrm{DIS}}=D_1^{\mathrm{DIS}}=0`$, although the conviction that $`D^{\mathrm{DIS}}`$ vanishes to all orders was not made.

### 2.2 The physical kernel

It is instructive to formulate first Sudakov resummation in the physical kernel defined by the logarithmic derivative $`d\mathrm{ln}F_2^N(Q^2)/d\mathrm{ln}Q^2`$. We claim that the latter can be written, to any order in perturbation theory, as

$`{\displaystyle \frac{d\mathrm{ln}F_2^N(Q^2)}{d\mathrm{ln}Q^2}}={\displaystyle \frac{d\mathrm{ln}J_N(Q^2;\mu _F^2)}{d\mathrm{ln}Q^2}}={\displaystyle \frac{C_F}{\beta _0}}{\displaystyle _0^1}{\displaystyle \frac{x^{N1}1}{1x}}𝒦\left((1x)Q^2\right)`$
$`={\displaystyle \frac{C_F}{\beta _0}}{\displaystyle _0^{\mathrm{}}}𝑑uT(u)\left(Q^2/\mathrm{\Lambda }^2\right)^u\left(N^u1\right)\mathrm{\Gamma }(u)B[𝒦](u),`$ (6)
where, as before, terms that are finite at large $`N`$ are neglected. Here we used the scheme invariant Borel representation ,

$$𝒦(\mu ^2)=_0^{\mathrm{}}𝑑uB[𝒦](u)T(u)\left(\mu ^2/\mathrm{\Lambda }^2\right)^u,$$
(7)
and $`T(u)`$ is the Laplace transform of the coupling555Eq. (2.2) is written in a renormalization scheme invariant way. Both $`T(u)`$ and the coefficients of $`B[𝒦](u)`$ would change depending on the renormalization group equation for the coupling, such that the product remains invariant.,

$$A(Q^2)=_0^{\mathrm{}}𝑑uT(u)\left(Q^2/\mathrm{\Lambda }^2\right)^u.$$
(8)
In the following we will use the ‘t Hooft coupling $`A\alpha _s(Q^2)\beta _0/\pi `$, defined by

$$\frac{dA}{d\mathrm{ln}Q^2}=A^2\left(1+\delta A\right);\delta \beta _1/\beta _{0}^{}{}_{}{}^{2},$$
(9)
so

$$T(u)=(u\delta )^{u\delta }\mathrm{exp}(u\delta )/\mathrm{\Gamma }(1+u\delta ).$$
(10)
The function

$$B[𝒦](u)=1+\stackrel{~}{k}_1\frac{u}{1!}+\stackrel{~}{k}_2\frac{u^2}{2!}+\stackrel{~}{k}_3\frac{u^3}{3!}+\mathrm{}$$

is the Borel representation of an $`N`$-independent effective charge, which is a function of the coupling having the expansion

$$𝒦=A+k_1A^2+k_2A^3+\mathrm{}.$$

Our notation here is adopted for the large-$`\beta _0`$ limit, where $`\delta =0`$, $`T(u)=1`$ and $`\stackrel{~}{k}_n=k_n`$. Beyond this limit the general relation can be written as an expansion in $`\delta `$,

$$\stackrel{~}{k}_n=k_n\delta k_{n1}n\left(\mathrm{\Psi }(n+1)+\gamma _E1\right)+𝒪(\delta ^2).$$
(11)
Like the physical kernel itself, $`𝒦`$ is does not depend on any arbitrary factorization procedure. Note that even if $`𝒦`$ is free of any infrared renormalons (as is the case in the large-$`\beta _0`$ limit, and perhaps also in the full theory – see below), these do appear in $`d\mathrm{ln}F_2^N(Q^2)/d\mathrm{ln}Q^2`$ through the factor $`\left(N^u1\right)\mathrm{\Gamma }(u)`$ which is associated with soft and collinear radiation under the large-$`x`$ phase-space constraint . This factor is understood as follows: gluon splitting off a quark gives rise to a $`1/(1x)_+`$ singularity. The phase-space constraint determining the maximal gluon virtuality is $`k^2<W^2=Q^2(1x)/x`$. Using the Borel representation of the dressed gluon propagator, $`k^2`$ is raised to a power $`u`$. Performing the phase-space integration, phase-space constraint promotes the $`x1`$ singularity to a cut: $`(1x)^{1u}`$. Going to moment space and neglecting $`𝒪(1/N)`$ contributions one obtains: $`(N^u1)\mathrm{\Gamma }(u)`$. As a result, at any positive integer value of $`u`$ where $`B[𝒦](u)`$ does not vanish, the physical kernel will have a renormalon singularity.

The all-order $`B[𝒦](u)`$ in the large-$`N_f`$ limit has been computed in , and it is given by

$$B[𝒦](u)|_{\mathrm{large}\beta _0}=\frac{1}{2}\frac{\mathrm{sin}\pi u}{\pi u}\mathrm{e}^{Cu}\left[\frac{1}{1u}+\frac{1}{1u/2}\right],$$
(12)
where $`C=5/3`$ is the constant from the $`\overline{\mathrm{MS}}`$ renormalization of the fermion loop666We assume that $`\mathrm{\Lambda }`$ in (2.2) is defined in this renormalization scheme, while the renormalization-group equation is (9).. Additionally, the first two coefficients $`k_1`$ and $`k_2`$ are known (see e.g. ), by combining results for the coefficient function and the anomalous dimension in the $`\overline{\mathrm{MS}}`$ factorization scheme. The relations will be given below.

Examining the (large $`N_f`$) renormalon structure of the physical kernel (2.2) with (12) we deduce that non-perturbative corrections at large $`N`$ appear as  a multiplicative factor in moment space, as summarized by eq. (1), with

$$J^{\text{NP}}(N\mathrm{\Lambda }^2/Q^2)=\mathrm{exp}\left\{\omega _1\frac{C_F}{\beta _0}\left(\frac{N\mathrm{\Lambda }^2}{Q^2}\right)\frac{1}{2}\omega _2\frac{C_F}{\beta _0}\left(\frac{N\mathrm{\Lambda }^2}{Q^2}\right)^2\right\}.$$
(13)
This power-correction model effectively resums the twist expansion, under the assumption that the higher-twist contributions which dominate at large $`N`$, are those which mix with the leading twist. Fixing the parameters $`\omega _1`$ and $`\omega _2`$ requires the knowledge of the dominant contribution to the matrix elements of the corresponding higher-twist operators . Clearly, the non-perturbative contribution $`J^{\text{NP}}(N\mathrm{\Lambda }^2/Q^2)`$ has a meaning only within a given regularization prescription for the renormalons in the perturbative jet function $`J_N(Q^2;\mu _F^2)`$.

In the large-$`\beta _0`$ limit $`𝒦`$ is free of any renormalon ambiguity. $`B[𝒦](u)`$ is an analytic function at any finite $`u`$, so renormalon singularities in the physical kernel are exclusively related to the end-point singularity at $`x1`$. If $`𝒦`$ is free of infrared renormalons also in the full theory, it follows777As usual , the nature of infrared renormalon singularities is related with the anomalous dimension of higher-twist operators. that at large $`N`$ the anomalous dimension of higher-twist operators (which mix with twist two), coincides with that of the twist two, as was already anticipated888Let us recall that the anomalous dimension of twist three is known to coincide asymptotically with that of twist two . in . Another possibility not yet excluded is that beyond the large-$`\beta _0`$ limit, $`B[𝒦](u)`$ would have infrared renormalons. For example, the poles at $`u=1`$ and $`u=2`$, which are exactly compensated by the sine factor in (12), might become branch points. Such a cut structure would imply that the large $`N`$ asymptotics of higher-twist anomalous dimensions differ from that of the leading twist, and a more complicated pattern of power corrections at large $`N`$ would emerge.

Finally, let us return to the last point we made on section 2.1. The absence an additional source of Sudakov logs of the form $`𝒟^{\mathrm{DIS}}((1x)^2Q^2)`$ in inclusive DIS (such a term was included, for example, in ) can also be understood from general consideration on power corrections. In DIS the external Lorentz invariants are $`Q^2`$ and $`W^2Q^2/N`$. On these grounds alone one expects power corrections to appear on either of these two scales and not on any other scale. For example, corrections such as $`\mathrm{\Lambda }^2N^2/Q^2`$ should not appear. On the other hand, had there been a function $`𝒟^{\mathrm{DIS}}((1x)^2Q^2)`$ in eq. (2.2), renormalons associated with the scale $`Q^2/N^2`$ would have appeared999The corresponding Borel integrand includes the factor $`(N^{2u}1)\mathrm{\Gamma }(2u)`$, so unless $`B[𝒟](u)`$ vanishes at all integer and half integer values of $`u`$, there would be power corrections on the scale $`Q^2/N^2`$. upon performing the integration over $`x`$, similarly to the situation in Drell-Yan production . As we argued in section 2.1, the perturbative factorization of structure functions does not require to introduce such a function and this problem is avoided altogether.

### 2.3 The jet function

Going from the physical kernel to $`\mathrm{ln}F_2^N`$ itself, the differential equation (2.2) should be supplemented by an initial condition. Within the OPE, the initial condition is defined as the matrix element of the twist-two quark operator with $`N`$ covariant derivatives, $`q_N(\mu _F^2)`$. The operator has a non-trivial renormalization, which takes the form

$$\frac{d\mathrm{ln}q_N(\mu _F^2)}{d\mathrm{ln}\mu _F^2}=\frac{C_F}{\beta _0}\mathrm{ln}N𝒜$$
(14)
up to corrections which are finite at large $`N`$. The effective charge

$$𝒜=A+a_2A^2+a_3A^3+\mathrm{}$$

is associated with the soft function $`V`$ in section 2.1. It corresponds to the $`1/(1x)_+`$ singular part in the Altarelli-Parisi splitting function, which can be computed using the Eikonal approximation, or the Wilson lines shown in the lower picture in Fig. 2. It has been noted that $`𝒜`$ is a rather universal object : it is the anomalous dimension of a Wilson line with a cusp. Upon integrating (2.2) and (14) and using the factorization formula (1) to write $`\mathrm{ln}J_N`$ as the difference between $`\mathrm{ln}F_2^N`$ and $`\mathrm{ln}q_N`$ (both regularized by the Borel variable $`u0`$), we obtain the following jet function:

$`\mathrm{ln}J_N(Q^2;\mu _F^2)`$ $`=`$ $`{\displaystyle \frac{C_F}{\beta _0}}{\displaystyle _0^{\mathrm{}}}𝑑uT(u)\left(Q^2/\mathrm{\Lambda }^2\right)^u`$ (15)
$`\left[\mathrm{\Gamma }(u)\left(N^u1\right){\displaystyle \frac{B[𝒦](u)}{u}}+\left(\mu _F^2/Q^2\right)^u{\displaystyle \frac{B[𝒜](u)}{u}}\mathrm{ln}N\right],`$
where, as before we use the scheme invariant Borel representation, where

$$B[𝒜](u)=1+\stackrel{~}{a}_2\frac{u}{1!}+\stackrel{~}{a}_3\frac{u^2}{2!}+\mathrm{}.$$
(16)
Note that the second term in the squared brackets cancels the $`u=0`$ singularity of the first, so in $`\mathrm{ln}J_N(Q^2;\mu _F^2)`$ the integral over $`u`$ exists. This is in contrast with $`\mathrm{ln}q_N`$ and $`\mathrm{ln}F_2^N`$, which are not perturbative entities. The dependence of $`\mathrm{ln}J_N(Q^2;\mu _F^2)`$ on the factorization scale $`\mu _F`$ and on the factorization scheme is guaranteed to cancel that of $`q_N(\mu _F^2)`$.

Eq. (15) summarises the perturbative jet function of $`F_2`$, namely the $`\mathrm{ln}N`$ enhanced contribution to the corresponding coefficient function. It is written in an all-order, scheme-invariant way. There are various ways to define the Borel integral in (15). In general, different regularizations of the renormalons differ by power terms of the form (13). Therefore, once power terms are included it does not matter which particular regularization is chosen : difference between different regularizations amounts to redefining $`\omega _i`$.

Here we choose to define (15) by truncating the corresponding perturbative series

$$\mathrm{ln}J_N(Q^2;\mu _F^2)=\frac{C_F}{\beta _0}\underset{k=0}{\overset{\mathrm{}}{}}\underset{l=1}{\overset{k+1}{}}C_{k,l}(\mathrm{ln}N)^lA^k,$$
(17)
with $`A\beta _0\alpha _s/\pi `$, at the minimal term, i.e. at the point where the terms start increasing in magnitude. Even this can be done is several ways, which differ by power corrections. One possibility would be to perform the sum over $`l`$ for a given $`k`$ and truncate the sum over $`k`$ at the minimal term. We will use another: let us rewrite the sum as in resummation with a fixed logarithmic accuracy,

$$\mathrm{ln}J_N(Q^2;\mu _F^2)=\frac{C_F}{\beta _0}\underset{m=0}{\overset{\mathrm{}}{}}g_m(\lambda )A^{m1},$$
(18)
where $`g_m(\lambda )`$ are functions of $`\lambda =A\mathrm{ln}N`$, and truncate the sum over $`m`$ at the minimal term.

It is straightforward to expand the Borel integrand in (15), and derive explicit expressions for $`g_m(\lambda )`$. Let us choose $`\mu _F^2=Q^2`$ and write $`g_m(\lambda )`$ as a series in $`\delta =\beta _1/\beta _0^2`$:

$$g_m(\lambda )g_m^{(1)}(\lambda )+\delta g_m^{(2)}(\lambda )+\delta ^2g_m^{(3)}(\lambda )+\mathrm{},$$
(19)
where first term appears in the large-$`\beta _0`$ limit (setting $`T(u)=1`$) and second and on are associated with integrating over the two-loop coupling using (10). The expansion in powers of $`\delta `$ is useful since we will eventually discard higher powers of $`\delta `$, which enter for the first time at rather high orders (note that $`\delta `$ itself is not small $`\delta (N_f=4)=0.7392`$).

The functions $`g_m^{(1)}(\lambda )`$ and $`g_m^{(2)}(\lambda )`$ are given by

$`g_m^{(1)}(\lambda )`$ $`=`$ $`r_m(\lambda )\stackrel{~}{p}_m/m!+\theta (m1)\lambda (\stackrel{~}{p}_m\stackrel{~}{a}_{m+1})/m`$ (20)
$`g_m^{(2)}(\lambda )`$ $`=`$ $`\theta (m1)\stackrel{~}{r}_m(\lambda )\stackrel{~}{p}_{m1}/(m1)!+\theta (m2)\lambda (\psi (m)+\gamma _E1)(\stackrel{~}{p}_{m1}\stackrel{~}{a}_m),`$
where the coefficients $`\stackrel{~}{p}_n`$ ($`𝒫=A+p_1A^2+p_2A^3+\mathrm{}`$) are defined through the expansion of

$$B[𝒫](u)\mathrm{\Gamma }(1u)B[𝒦](u)=1+\stackrel{~}{p}_1\frac{u}{1!}+\stackrel{~}{p}_2\frac{u^2}{2!}+\mathrm{}.$$
(21)
In (20), aside from a linear term associated with the factorization procedure, the $`\lambda `$-dependence appears only through the functions

$`r_m(\lambda )`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(m+n)!}{(n+2)!}}\lambda ^{n+2}=\mathrm{\Gamma }(m1)\left[(1\lambda )^{1m}\lambda (m1)1\right]`$
$`\stackrel{~}{r}_m(\lambda )`$ $`=`$ $`\left({\displaystyle \frac{d}{dm}}+\gamma _E1\right)r_m(\lambda ),`$ (22)
which are defined101010For $`m=0,1`$ the limit should be taken carefully, yielding logarithms of $`(1\lambda )`$. for $`m0`$. Note that both $`r_m(\lambda )`$ and $`\stackrel{~}{r}_m(\lambda )`$ are free of a linear term: their expansions start at order $`\lambda ^2`$. These functions have two generic properties  associated with their origin in an integral over the running coupling: they increase factorially with $`m`$ due to infrared renormalons and posses an increasing singularity at $`\lambda =1`$, corresponding to $`Q^2/N\mathrm{\Lambda }^2`$, where the perturbative treatment looses its validity.

Terms with $`\delta ^2`$ enter first at NNLL accuracy ($`m=2`$):

$$g_2^{(3)}(\lambda )=\frac{1}{2}\frac{(\lambda +\mathrm{ln}(1\lambda ))^2}{1\lambda }.$$
(23)
The lowest order contribution here is $`\mathrm{\Delta }\mathrm{ln}J_N(Q^2;\mu _F^2)=(C_F/\beta _0)\frac{1}{8}\delta ^2A^5\mathrm{ln}^4N`$, where we kept only the highest power of $`\mathrm{ln}N`$. Contrary to other terms of the same logarithmic accuracy (NNLL) this term has no large numerical coefficients. Normalizing by the magnitude of the leading term in $`\mathrm{ln}J_N(Q^2;\mu _F^2)=(C_F/\beta _0)\frac{1}{2}A\mathrm{ln}^2N`$, and substituting $`A\mathrm{ln}N=1`$ (for a worst-case estimate) the corrections is of relative size of $`(\delta A/2)^2`$, i.e. $`0.3`$% at $`Q^2=10`$ GeV2. In the following we will work with exact NNLL resummation, supplemented by estimates of subleading logs (N3LL and beyond, up to the minimal term in the series). As we will see below, estimates of subleading logs are anyway of limited accuracy. So there we will take into account just the first two terms in (19), neglecting contributions $`𝒪(\delta ^2)`$, which are tiny. Note that dealing with N3LL and beyond we will also be neglecting running-coupling effects associated with four-loop contributions to the $`\beta `$ function, as these corrections are neglected here in the relation between the ‘t Hooft coupling and $`\overline{\mathrm{MS}}`$.

It is useful to relate the Borel representation $`B[𝒦](u)`$ to the commonly used functions $`𝒜`$ and $``$. The relation is:

$$B[𝒦](u)=B[𝒜](u)uB[](u)$$
(24)
or

$$𝒦=𝒜+\frac{d}{d\mathrm{ln}Q^2}.$$
(25)
If we write $`=b_1A+b_2A^2+b_3A^3+\mathrm{}`$ with

$$B[](u)=\stackrel{~}{b}_1+\stackrel{~}{b}_2\frac{u}{1!}+\stackrel{~}{b}_3\frac{u^2}{2!}+\mathrm{}.$$
(26)
then $`\stackrel{~}{k}_n=\stackrel{~}{a}_{n+1}n\stackrel{~}{b}_n`$. Finally, it is straightforward to write the relation with $`a_n^{\overline{\text{MS}}}`$ and $`b_n^{\overline{\text{MS}}}`$ defined in the $`\overline{\mathrm{MS}}`$ renormalization scheme111111Note that our notation for the coefficients of $`𝒜`$ and $``$ is somewhat different from . The relation is $`A_n=C_F4^n\beta _0^{n1}a_n^{\overline{\text{MS}}}`$ and similarly for $`B_n`$.. For the first two orders (which are known in QCD) we have

$`\begin{array}{ccc}\stackrel{~}{a}_1=a_1=1\hfill & \stackrel{~}{b}_1=b_1=3/4\hfill & \stackrel{~}{k}_0=1\hfill \\ \stackrel{~}{a}_2=a_2=a_2^{\overline{\text{MS}}}\hfill & \stackrel{~}{b}_2=b_2=b_2^{\overline{\text{MS}}}\hfill & \stackrel{~}{k}_1=\stackrel{~}{a}_2\stackrel{~}{b}_1\hfill \\ \stackrel{~}{a}_3=a_3\delta a_2=a_3^{\overline{\text{MS}}}+\delta _2^{\overline{\text{MS}}}\delta a_2^{\overline{\text{MS}}}\hfill & & \stackrel{~}{k}_2=\stackrel{~}{a}_32\stackrel{~}{b}_2,\hfill \end{array}`$ (30)
where $`\delta _2^{\overline{\text{MS}}}=\beta _2^{\overline{\text{MS}}}/\beta _0^3`$. The $`\overline{\mathrm{MS}}`$ coefficients are given by ,

$`a_2^{\overline{\text{MS}}}(N_c,\beta _0)`$ $`=`$ $`{\displaystyle \frac{5}{3}}+\left({\displaystyle \frac{1}{3}}{\displaystyle \frac{\pi ^2}{12}}\right){\displaystyle \frac{C_A}{\beta _0}},`$ (31)
$`a_3^{\overline{\text{MS}}}(N_c=3,\beta _0)`$ $`=`$ $`{\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{\beta _0}}\left(\left({\displaystyle \frac{55}{16}}3\zeta _3\right)C_F+\left({\displaystyle \frac{253}{72}}{\displaystyle \frac{5}{18}}\pi ^2+{\displaystyle \frac{7}{2}}\zeta _3\right)C_A\right){\displaystyle \frac{19.5\pm 0.2}{\beta _0^2}}`$
$`b_2^{\overline{\text{MS}}}(N_c,\beta _0)`$ $`=`$ $`{\displaystyle \frac{247}{72}}+{\displaystyle \frac{1}{6}}\pi ^2+{\displaystyle \frac{1}{\beta _0}}\left(\left({\displaystyle \frac{3}{32}}{\displaystyle \frac{3}{2}}\zeta _3+{\displaystyle \frac{1}{8}}\pi ^2\right)C_F+\left({\displaystyle \frac{73}{144}}+{\displaystyle \frac{5}{2}}\zeta _3\right)C_A\right),`$
where the purely non-Abelian ingredient in $`a_3^{\overline{\text{MS}}}`$ is still not known exactly. However, a reliable estimate  was extracted from the first few moments  for $`N_c=3`$, allowing us to write the above.

### 2.4 Naive non-Abelianization

As explained in the introduction, our task is to take into account all the perturbative information available on $`\mathrm{ln}J`$, namely the large-$`N_f`$ limit and the first few coefficients known in QCD. The first question is, of course, to what extent does “naive non-Abelianization” provide an estimate of the coefficients. The answer, as reflected in the numerical values of the first few coefficients for $`N_c=3`$,

$`a_2^{\overline{\text{MS}}}`$ $`=`$ $`1.667{\displaystyle \frac{1.467}{\beta _0}},`$ (32)
$`a_3^{\overline{\text{MS}}}`$ $`=`$ $`0.333+{\displaystyle \frac{14.714}{\beta _0}}{\displaystyle \frac{19.5\pm 0.2}{\beta _0^2}}`$
$`b_2^{\overline{\text{MS}}}`$ $`=`$ $`1.786+{\displaystyle \frac{6.6104}{\beta _0}},`$
is quite clear: the terms that are leading in the large-$`\beta _0`$ limit, are not at all dominant for the realistic value $`\beta _0(N_f=4)2`$.

This is also reflected in Table 1 (compare, for example, the $`n=2`$ values in the lower and upper tables). One simply cannot use the large $`\beta _0`$ limit directly to provide an estimate121212This negative statement does not apply to physical quantities such as the first few moments of the structure functions. In this case estimates based on “naive non-Abelianization” are reasonable . of $`𝒜`$, $``$ or $`𝒦`$ because these functions do not contain renormalons, and thus, there is no reason that the large-$`\beta _0`$ terms will have large coefficients. In fact, it is known that the series for the cusp anomalous dimension $`𝒜`$ converges in this limit into ,

$$𝒜|_{\mathrm{large}\beta _0}=\frac{\mathrm{sin}\pi A}{\pi }\frac{\mathrm{\Gamma }(4+2A)}{6\mathrm{\Gamma }(2+A)^2}+\mathrm{},$$
(33)
where the $`\overline{\mathrm{MS}}`$ scheme is used, and the dots stand for terms that are subleading at large $`\beta _0`$. It follows from (12) and (24) that $``$, is renormalons-free as well131313Although the coefficients of $``$ are known to all orders, a closed-form expression as a function of the coupling is not known.. Thus, also here “naive non-Abelianization” is not expected to apply.

The status of $`𝒫`$ is somewhat different. It does contains renormalons so here the terms that are leading in the large-$`\beta _0`$ limit are expected to dominate at large orders. This does not necessarily imply that “naive non-Abelianization” estimates will be reliable at low orders, but it is expected that the general trend of increase characterizing $`p_n`$ in the large-$`\beta _0`$ limit right from the start, will persist in the full theory. The first few entries in the lower Table 1 are consistent with this expectation.

While the large-$`\beta _0`$ coefficients themselves cannot be used to estimate the actual QCD ones, more general properties may be deduced. As reflected in Table 1 neither $`𝒜`$ nor $``$ has a regular sign pattern. The sign is flipped every three or four terms owing to the $`\mathrm{sin}\pi u`$ function in (12) associated with taking a final-state cut. Another property is that in the large-$`\beta _0`$ limit the coefficients of $`𝒜`$ decrease with order while those of $``$ have a general tendency to increase. Starting at $`n\begin{array}{c}>\hfill \\ \hfill \end{array}4`$ there is a clear hierarchy $`\left|a_{n+1}\right|\left|k_n\right|`$, i.e. the physical kernel is dominated by $``$: $`k_n=a_{n+1}nb_nnb_n`$. An even stronger hierarchy exists between $`𝒜`$ and $`𝒫`$: $`\left|a_{n+1}\right|\left|p_n\right|`$, as the latter contains renormalons.

Finally, to make concrete predictions for the subleading logs in (20), the knowledge of $`\stackrel{~}{p}_n\stackrel{~}{a}_{n+1}`$ and $`\stackrel{~}{p}_n`$ at higher orders ($`n3`$) is required. However, from the discussion above it is expected that at higher orders $`\left|\stackrel{~}{a}_{n+1}\right|\left|\stackrel{~}{p}_n\right|`$, and thus the knowledge of $`𝒜`$ to higher orders is absolutely inessential. As usual, statements about large orders do not necessarily apply for any $`n3`$.

Indeed, as shown in Table 2 the expected hierarchy is not realised for $`n<3`$. Yet, the $`n=3`$ entries in Table 1 are reassuring: they suggest an order of magnitude difference between $`\stackrel{~}{p}_3`$ and $`\stackrel{~}{a}_4`$. We will therefore orient ourselves in the next section to estimating $`𝒦`$ and thus $`𝒫`$ in full QCD while for $`𝒜`$ we will simply use the first known coefficients ($`\stackrel{~}{a}_2`$ and $`\stackrel{~}{a}_3`$) supplemented by the large-$`\beta _0`$ contribution (although the latter does not dominate!) for the higher orders. An alternative is to neglect $`\stackrel{~}{a}_n`$ for $`n4`$ altogether.

### 2.5 Going beyond the large-$`\beta _0`$ limit

A major advantage of the large-$`\beta _0`$ result for $`𝒦`$ is that its analytic structure in the Borel plane is known. We will use this structure as a basis for modeling the Borel function in the full theory. We shall require consistency with the first few orders that are known, such that NNLL accuracy is guaranteed. By examining different models we would like to check the sensitivity to yet unknown higher-order coefficients, in particular in the context of open theoretical problems concerning renormalons. There is no attempt here to cover the infinite set of functions that obey the constraints.

We will assume that the singularity structure of the physical kernel is similar to that obtained in the large-$`\beta _0`$ limit from (12) and (2.2). In particular, in our models there will be two renormalons at $`u=1`$ and $`u=2`$, but no higher ones. Thus, non-perturbative contributions will be parametrized by (13). This is not to say that higher renormalon singularities do not exist in the full theory, but rather that with the information available there is no way we can control their residues. Thus, in our models, similarly to (12), $`B[𝒦](u)`$ does not vanish at $`u=1,2`$, but it does vanish at all higher integers.

We will also restrict ourselves to the case where $`B[𝒦](u)`$ is regular at $`u=1,2`$ so the renormalon singularities of the physical kernel in our models appear as simple poles. As mentioned in section 2.2, the nature of the singularities is related with the anomalous dimension of higher-twist operators. Here, we keep $`B[𝒦](u)`$ renormalon free, thus assuming that the anomalous dimension at any twist coincides at large $`N`$ with that of twist two. This issue certainly deserves further investigation in the future.

We investigate two types of models for $`B[𝒦](u)`$. The first is based on modifying the large-$`\beta _0`$ function (12) by a multiplicative factor and the second – by an additive term. A simple example with a multiplicative factor is

$`B[𝒦](u)|_{\mathrm{model}a}={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{sin}\pi u}{\pi u}}\mathrm{e}^{\frac{5}{3}u}\left[{\displaystyle \frac{1}{1u}}+{\displaystyle \frac{1}{1u/2}}\right]\times \mathrm{e}^{u\left(a_2\frac{5}{3}\right)}(1+H^{(a)}u^2),`$
$`H^{(a)}={\displaystyle \frac{3}{4}}a_2{\displaystyle \frac{5}{8}}{\displaystyle \frac{1}{2}}a_2^2+{\displaystyle \frac{1}{6}}\pi ^2+{\displaystyle \frac{1}{2}}a_3^{\overline{\text{MS}}}+{\displaystyle \frac{1}{2}}\delta _2^{\overline{\text{MS}}}{\displaystyle \frac{1}{2}}\delta a_2b_2`$ (34)
where $`H^{(a)}`$ is fixed to reproduce the exact $`\stackrel{~}{k}_1`$ and $`\stackrel{~}{k}_2`$ of (30). Another example (model $`b`$) is provided by replacing $`(1+H^{(a)}u^2)`$ by $`\mathrm{e}^{H^{(b)}u^2}`$, with $`H^{(b)}=H^{(a)}`$. Next, a couple of simple examples with additive (regular) terms are

$`B[𝒦](u)|_{\mathrm{model}c}={\displaystyle \frac{\mathrm{sin}\pi u}{\pi u}}\left\{{\displaystyle \frac{1}{2}}\mathrm{e}^{\frac{5}{3}u}\left[{\displaystyle \frac{1}{1u}}+{\displaystyle \frac{1}{1u/2}}\right]+u\left(a_2{\displaystyle \frac{5}{3}}\right)+u^2H^{(c)}\right\},`$
$`H^{(c)}={\displaystyle \frac{235}{72}}+{\displaystyle \frac{1}{6}}\pi ^2+{\displaystyle \frac{1}{2}}a_3^{\overline{\text{MS}}}+{\displaystyle \frac{1}{2}}\delta _2^{\overline{\text{MS}}}{\displaystyle \frac{1}{2}}\delta a_2b_2`$ (35)
and

$`B[𝒦](u)|_{\mathrm{model}d}={\displaystyle \frac{\mathrm{sin}\pi u}{\pi u}}\left\{{\displaystyle \frac{1}{2}}\mathrm{e}^{\frac{5}{3}u}\left[{\displaystyle \frac{1}{1u}}+{\displaystyle \frac{1}{1u/2}}\right]+\mathrm{e}^{u(a_2\frac{5}{3})+u^2H^{(d)}}1\right\},`$
$`H^{(d)}={\displaystyle \frac{335}{72}}+{\displaystyle \frac{1}{6}}\pi ^2+{\displaystyle \frac{1}{2}}a_3^{\overline{\text{MS}}}+{\displaystyle \frac{1}{2}}\delta _2^{\overline{\text{MS}}}{\displaystyle \frac{1}{2}}\delta a_2b_2{\displaystyle \frac{1}{2}}a_2^2+{\displaystyle \frac{5}{3}}a_2.`$ (36)
The main difference between the two sets of models, is that models $`c`$ and $`d`$ keep the residues of the renormalon singularities at their large-$`\beta _0`$ values, whereas in models $`a`$ and $`b`$ these residues are modified at NLL by purely non-Abelian contributions, proportional to $`C_A/\beta _0`$, and, starting from NNLL by correlated multi-gluon emission even in the Abelian case. In general the residues should be modified by such contributions, but the way the functional form varies is not yet under theoretical control.

As shown in figure 4 different models extrapolating from the large-$`\beta _0`$ limit have very different behaviour indeed. The knowledge of the value of $`B[𝒦](u)`$ at $`u=0`$ and the first two derivative there cannot constrain it away from this point.

Note that models with an additive correction ($`c`$ and $`d`$) follow closely the large-$`\beta _0`$ functional form away from $`u=0`$, as their residues at $`u=1`$ and at $`u=2`$ are fixed by this limit. On the other hand, models with a multiplicative correction ($`a`$ and $`b`$) can have a completely different form and even the sign of the first renormalon residue varies.

In general, the detailed behaviour of $`B[𝒦](u)`$ away from the origin, has a minor significance concerning the value of the physical kernel itself, or $`J_N(Q^2;\mu _F^2)`$. To see the sensitivity to the behaviour of $`B[𝒦](u)`$, we computed $`\mathrm{ln}J_N(Q^2;\mu _F^2)`$ with the different models described above. The results are shown in Fig. 5 order by order at increasing logarithmic accuracy. Differences are rather significant. Models $`a`$ and $`b`$ tend to have a moderate (although not negligible) contributions beyond the NNLL, quite a favourable scenario for perturbation theory. On the other hand, models $`c`$ and $`d`$ indicate huge corrections beyond the NNLL. The resolution of the theoretical questions allowing to pick the right model for subleading logarithmic corrections thus becomes crucial for any analysis oriented at power accuracy.

One typical feature of models $`c`$ and $`d`$ is that contributions of consequential orders are not monotonically decreasing up to the minimal term, but instead, their magnitude oscillates. This, however, should not be too surprising given that in the case under consideration oscillations in magnitude of consequential terms occur at the first few order: the NLL contribution is larger than the LL one, while the NNLL contribution is smaller than both.

As mentioned above, already when constructing a model for $`B[𝒦](u)`$ which is consistent with NLL accuracy in the non-Abelian theory, there are various possibilities. Nevertheless, in this case it seems natural to use a multiplicative modification of the large-$`\beta _0`$ result, such that $`C=5/3`$ at the exponent in (12) is effectively replaced for the full $`a_2`$ coefficient (31). In fact, the choice made when promoting the large-$`\beta _0`$ result to NLL accuracy is the most significant one. This can be verified by examining the model,

$$B[𝒦](u)|_{\mathrm{model}e}=\frac{\mathrm{sin}\pi u}{\pi u}\left\{\frac{1}{2}\mathrm{e}^{a_2u}\left[\frac{1}{1u}+\frac{1}{1u/2}\right]+\mathrm{e}^{u^2H^{(e)}}1\right\}$$
(37)
with $`H^{(e)}=H^{(a)}`$, where $`5/3`$ is replaced by the full $`a_2`$ in the first term, while further corrections (NNLL) are incorporated as an additive term. This intermediate step between the two classes of models described above, turns out to be quite close to the ones with a multiplicative correction, $`a`$ and $`b`$. For example, for $`N=10`$, $`Q^2=20`$ GeV2, upon truncation at the minimal term model $`a`$ yields $`\mathrm{ln}J=0.713`$, model $`d`$ yields $`\mathrm{ln}J=0.868`$ and model $`e`$ yields $`\mathrm{ln}J=0.755`$.

Indeed, in previous applications of DGE , the procedure chosen to extrapolate from the large-$`\beta _0`$ limit and comply with NLL accuracy was the replacement $`5/3a_2`$. This follows the spirit of  who first noted the special rôle of the “gluon bremsstrahlung” effective charge, (defined by the cusp anomalous dimension $`𝒜`$) in this context. This natural generalization is sufficient at NLL, but not beyond it. In any case, the deeper theoretical problem of how purely non-Abelian contributions as well as correlations in multiple gluon emission modify the renormalon structure of the Sudakov exponent remains open.

With this caution in mind, we will proceed here on good faith with models with a multiplicative modification of (12), and in particular, with model $`a`$, to perform data analysis. A multiplicative modification can always be recast as a modification of the coupling which is integrated over by redefining $`T(u)`$ in (15), and thus it can be regarded as the most natural generalization of the large-$`\beta _0`$ limit result. Based on Fig. 5 it is also an optimistic choice, leading to rather moderate higher-order corrections.

## 3 Data analysis

In the previous section we summarized the state-of-the-art knowledge (and ignorance) on the twist-two coefficient function of $`F_2`$ at large $`x`$. We also have a picture  concerning the structure of higher twist, with a concrete formula for their parametrization (13). We now turn to investigate the implications this has on data analysis. The following questions will be emphasized:

* Are the data and the theoretical description of $`F_2`$ consistent?
* Is there any experimental indication for higher-twist contribution to $`F_2`$?
* To what extent can one determine the relevant parameters from the data. This includes $`\alpha _s`$, the non-singlet quark distribution $`q_N(\mu _F^2)`$ and the higher twist parameters $`\omega _i`$.
* What is the significance of Sudakov resummation as compared with the NNLO result for the coefficient function?
* What is the significance of renormalon resummation in the Sudakov exponent and the current uncertainty about subleading logs (beyond NNLL)?

In order to systematically treat target-mass corrections and simplify the resummation procedure, we work here in moment space. For this purpose we compute the Nachtmann moments from the available experimental data on $`F_2^p(x,Q^2)`$. We also include the contributions at $`x=1`$ from the proton electric and magnetic form factors. The data we use are from a series of SLAC experiments on $`ep`$ scattering  and the BCDMS collaboration using $`\mu p`$ scattering .

Before coming to the fits, let us briefly describe the data and the way we use it. The SLAC data cover a region from $`Q^2`$ below 1 GeV2 to 20 GeV2 and $`x`$ up to $`x_{\mathrm{max}}`$ defined by the minimum value of $`W^2=(m_p+m_\pi )^2`$. The BCDMS data cover a region $`x0.75`$ for $`20<Q^2<230`$ GeV2.

The moments are computed in $`15`$ bins of $`Q^2`$ from 2 to 105 GeV2. There is a region of overlap between the SLAC and BCDMS data sets. While small it does allow a check on the relative normalisation of the two experiments. The computation involves estimating the contribution from the region between the highest $`x`$ bin at each $`Q^2`$ value and $`x=x_{\mathrm{max}}`$. We take a linear extrapolation with an arbitrary uncertainty of $`100`$%. For the SLAC data this is a very small correction. The final errors on the values of each moment include this uncertainty.

The physics analysis we apply assumes simply the non-singlet evolution of moments. This is equivalent to assuming that high enough moments of $`F_2^p(x,Q^2)`$ are dominated by just the $`u`$ and $`d`$ valence quarks. This can be checked by taking a standard global fit to the DIS data  and evaluating the relative contribution to each moment from the sea and valence quarks. We find at $`Q^2`$ = 10 GeV2, for example, that the LO analysis of MRST2001 gives valence/total for $`N=4`$ to be over 90% and for $`N5`$ to be more than 95%. We therefore choose to analyse moments $`N5`$ and keep in mind the valence approximation as a source of an overall $`\begin{array}{c}<\hfill \\ \hfill \end{array}5`$% uncertainty.

### 3.1 NNLO based fits

It is well known  that a NLO analysis of DIS data fails to describe the large-$`x`$ region. Usually the deviations are attributed to higher-twist contributions at low $`W^2`$. In the MRST2001 NLO leading-twist analysis, data from $`W^212.5`$ GeV2 were simply dropped. Even including the NNLO corrections which are now available , though helping a little, did not succeed in diminishing the discrepancy at large $`x`$.

This situation is confirmed here by attempting to fit the moments we extracted by a NNLO formula,

$$F_2^N(Q^2)=q_N(\mu _F^2)E_N(Q^2,\mu _F^2)C_N(Q^2).$$
(38)
We use here the standard $`\overline{\mathrm{MS}}`$ factorization and renormalization schemes. The NNLO evolution factor is

$$E_N(Q^2,\mu _F^2)=\mathrm{exp}\left\{\frac{1}{2}_{\mu _F^2}^{Q^2}\gamma (\mu ^2)𝑑\mu ^2\right\}=\left[\frac{\alpha _s(Q^2)}{\alpha _s(\mu _F^2)}\right]^{d_N}[1+r_N\alpha _s(Q^2)+s_N\alpha _s^2(Q^2)],$$
(39)
and the NNLO coefficient function with the factorization scale (and renormalization scale) set to $`Q^2`$ is

$$C_N^{\text{NNLO}}(Q^2)=1+c_N^{(1)}\alpha _s(Q^2)+c_N^{(2)}\alpha _s^2(Q^2).$$
(40)
In (39), $`d_N=\gamma _N^{(0)}/2\beta _0`$, $`r_N`$ involves the NLO anomalous dimension $`\gamma _N^{(1)}`$ and $`s_N`$ the NNLO anomalous dimension $`\gamma _N^{(2)}`$. We use the results of Retey and Vermaseren  for the $`\gamma _N^{(2)}`$ for even $`N`$ and a smooth interpolation of these for odd $`N`$. The running of the coupling $`\alpha _s(Q^2)`$ is approximated by the NNLO $`\beta `$ function.

We fit each moment in turn with two parameters, the moment of the valence quark distribution $`q_N(\mu _F^2)`$ and $`\alpha _s^{\overline{\text{MS}}}(M_Z^2)`$. The resulting values of the latter are shown in Fig. 6 (upper box) showing two features: (i) large value of the coupling and (ii) a marked increase as $`N`$ increases. The MRST2001 NNLO fit with $`W^212.5`$ GeV2 prefers a much lower value, $`\alpha _s^{\overline{\text{MS}}}(M_Z^2)=0.1155`$. In Fig. 7 we show two NNLO descriptions, one with the pairs of parameters just described and the other corresponding to the extrapolation of the MRST2001 NNLO description to low $`W^2`$. We conclude, with no surprise, that the leading-twist NNLO formula cannot satisfactorily describe the large $`N`$ moments of $`F_2^p`$.

### 3.2 Resummation based fits

Next consider the additional effect of Sudakov resummation. This is done by modifying eq. (40) to include higher-order terms which are enhanced by powers of $`\mathrm{ln}N`$. We use, the so-called log-R matching, namely

$$\mathrm{ln}C_N(Q^2)=\mathrm{ln}C_N^{\text{NNLO}}(Q^2)+\mathrm{ln}J_N(Q^2)\mathrm{ln}C_N^{\text{overlap}}$$
(41)
where $`\mathrm{ln}J_N(Q^2)`$ is given by (15) or (18) with the factorization scale set to $`Q^2`$. The subtracted term accounts for log-enhanced contributions up to order $`\alpha _s^2`$ which are included in both the NNLO coefficient function and in the Sudakov resummed one,

$`\mathrm{ln}C_N^{\text{overlap}}`$ $`=`$ $`{\displaystyle \frac{C_F}{\beta _0}}[({\displaystyle \frac{1}{2}}L^2+(\gamma _Eb_1)L)A`$
$`+({\displaystyle \frac{1}{6}}L^3+({\displaystyle \frac{1}{2}}\gamma _E+{\displaystyle \frac{1}{2}}(a_2b_1))L^2+(\gamma _E(a_2b_1)b_2+{\displaystyle \frac{\gamma _E^2}{2}}+{\displaystyle \frac{\pi ^2}{12}})L)A^2].`$
We apply the two following criteria in selecting data for the fits. Firstly, we fit moments $`N5`$ , to ensure non-singlet dominance. Secondly, we limit $`\lambda =A\mathrm{ln}N=(\beta _0\alpha _s/\pi )\mathrm{ln}N\lambda _{\mathrm{max}}`$, to ensure the validity of the perturbative treatment. We choose $`\lambda _{\mathrm{max}}=0.4`$. Since $`\lambda `$ depends on one of our fit parameters, $`\alpha _s(M_Z^2)`$, the number of data points included in a fit could vary with this parameter. To avoid this we choose the data to be fitted for $`\lambda _{\mathrm{max}}`$ evaluated for the MRST2001 NNLO default value, $`\alpha _s(M_Z^2)=0.1155`$. The moments of the experimental data going into the fits are given in Table 3.

Note that we chose this same data set for the pure NNLO fits for the sake of comparison.

As in the NNLO fit we have two free parameters, $`q_N(\mu _F^2)`$ and $`\alpha _s^{\overline{\text{MS}}}(M_Z^2)`$. We choose three variants for Sudakov resummation. The first is the standard, NNLL accuracy resummation, and the other two include renormalon resummation according to models $`a`$ and $`d`$ discussed in section 2.5. Contrary to fixed logarithmic accuracy, renormalon resummation is oriented at power accuracy. The renormalons are regularized here simply by truncating the series (18) at the minimal term.

The results are shown in Fig. 6. In contrast to the NNLO case, the values of the coupling are practically independent of $`N`$. This is true for any of the renormalon-based models as well as for pure NNLL resummation. Thus, Sudakov resummation proves to be important: it provides the ingredient which is missing in the fixed-order perturbative description discussed above.

However, as shown in Fig. 6, the extracted value of $`\alpha _s(M_Z^2)`$ varies significantly depending on the chosen model. In pure NNLL resummation, typically $`\alpha _s(M_Z^2)=0.1155\pm 0.0010`$141414The error quoted, like the bars in Fig. 6, reflects just the propagated experimental error (taking into account the uncertainty in extrapolating the data to large $`x`$). It does not include any estimate of theoretical uncertainty, for example due to factorization or renormalization scale dependence. while with renormalon resummation $`\alpha _s(M_Z^2)`$ varies between $`0.1135\pm 0.0010`$ (model $`a`$) and $`0.1070\pm 0.0010`$ (model $`d`$). Since in any case NNLL accuracy is guaranteed, these differences are due to sub-leading logs, N3LL and beyond. Note that the result quoted for model $`a`$ actually represents a large class of models with a multiplicative modification of the large-$`\beta _0`$ all-order result. The details of this factor have a minor effect of the result. Taking the optimistic approach that these multiplicative models represent well the structure of the Sudakov exponent, the central value of the strong coupling from this analysis is $`0.1135`$. Clearly, more theoretical work is required to reduce the uncertainty.

Fig. 8 shows the results of fitting $`q_N`$ based on model $`a`$ for a fixed value of $`\alpha _s(M_Z^2)=0.1135`$. The plotted curves are almost identical to those obtained in a similar fit where $`\alpha _s`$ is free, since the variations of the best-fit value of $`\alpha _s`$ as a function of $`N`$ are small (Fig. 6). Note that in summing the series (18), the asymptotic expansion is terminated at the minimal term. A sensible estimate of the corresponding power-suppressed ambiguity is simply the size of the last term included in the sum. Fig. 8 shows the magnitudes of this ambiguity together with extrapolations of the fits to regions excluded from the fitting procedure. Under the assumption taken here that the dominant higher-twist terms are those associated with mixing with the leading twist, this power-suppressed ambiguity should provide an estimate of the magnitude of the higher twist. As clearly seen in Fig. 8 this is a small correction even in the highest moment. These results imply that the low $`W^2`$ data, i.e. $`W^212.5`$GeV2 can be successfully described upon improving the standard NNLO description by including Sudakov resummation. The resulting values of $`\alpha _s^{\overline{\text{MS}}}(M_Z^2)`$ are stable with $`N`$ and are close to the value extracted from the MRST NNLO analysis where a cut on low $`W^2`$ is imposed.

The other parameter $`q_N(\mu _F^2)`$ is shown in Fig. 9 for four values of $`\mu _F^2`$ for comparison with the corresponding values expected from the MRST NNLO analysis. For the most part the ranges of $`x`$ being explored in this analysis of the high moments correspond to ranges of $`W^2`$ below the value 12.5 GeV2 excluded by the MRST analysis. Thus Fig. 9 is an indication of how much the large $`x`$ parton distributions would have to be modified to take account of this resummation effect.

### 3.3 Fitting the higher twist

Finally, according to eq. (13), we include the simplest model for higher-twist corrections through an extra factor $`J^{\text{NP}}(N\mathrm{\Lambda }^2/Q^2)`$ where

$$J^{\text{NP}}(N\mathrm{\Lambda }^2/Q^2)=\mathrm{exp}\left\{\frac{C_{\text{HT}}N}{Q^2}\right\},$$
(43)
with $`C_{\text{HT}}\omega _1\frac{C_F}{\beta _0}\mathrm{\Lambda }^2`$. Here we neglect subleading corrections proportional to $`\omega _2`$.

We perform fits on the same set of moments as before, but with the extra parameter $`C_{\text{HT}}`$ characterising the additional $`Q^2/N`$ dependence from the higher twist. Naturally, there is a compensation between the magnitude of the parameters $`\alpha _s(M_Z^2)`$ and $`C_{\text{HT}}`$ which is shown in Fig. 10 in the lower box.

For a wide range of $`C_{\text{HT}}`$ there is little difference in the quality of the fits, see the upper box in Fig. 10. The conclusion is clear: there is no phenomenological evidence for a non-vanishing twist-four (or higher twist) contributions going beyond the resummed perturbation theory. On the other hand a positive $`C_{\text{HT}}`$ would decreases the already ‘lowish’ value of $`\alpha _s`$.

## 4 Conclusions

In this paper we investigated the large-$`x`$ limit of the structure function $`F_2(x,Q^2)`$, concentrating on the interplay between Sudakov resummation, renormalons and higher twist. The main assumption we make  is that the dominant higher-twist contribution in this region is related to the mixing of higher-twist operators with the leading twist. Consequently, our investigation focuses on constraining the twist-two coefficient function at large-$`x`$.

Our first point is that the physical kernel $`d\mathrm{ln}F_2^N/d\mathrm{ln}Q^2`$ receives log-enhanced contributions from a single source, the jet function depending on $`W^2=Q^2(1x)/x`$. This property is a consequence of factorization and it distinguishes structure functions from less inclusive quantities, such are event-shapes, which have an additional observable-dependent sensitivity to large angle soft emission.

Much is known about the structure of the Sudakov exponent of $`F_2`$ from general considerations, from the large-$`\beta _0`$ limit and from fixed-order calculations in dimensional regularization. Nevertheless, the investigation of Sec. 2.5 reveals that subleading logs, N3LL and beyond, which are not yet under theoretical control, may have a large effect. Further theoretical input on the all-order structure of the exponent is required to constrain these contributions.

While a pure NNLO formula without Sudakov resummation does not describe the large-$`x`$ data, once Sudakov resummation (NNLL) is employed the purely perturbative description of $`F_2`$ is fully consistent with the data. This is reflected in the good fit of the $`Q^2`$ dependence of each of the moments and in the stability of the extracted $`\alpha _s`$ as a function of $`N`$. This conclusion holds independently of whether renormalons are resummed and of what one assumes about subleading logs.

The theoretical uncertainty concerning the renormalon structure of the Sudakov exponent beyond the large-$`\beta _0`$ limit, as reflected in the models investigated in Sec. 2.5 and 3.2, translates into an uncertainty of order $`\pm 6`$% in the extracted value of the coupling. Thus, taking a conservative approach one would assign a rather large uncertainty to $`\alpha _s`$ from scaling violation of $`F_2^N(Q^2)`$ with $`N\begin{array}{c}>\hfill \\ \hfill \end{array}5`$. A more optimistic approach would be to assume that a multiplicative modification of the large-$`\beta _0`$ result, such as model $`a`$ in Eq. (2.5), represents well the all-order structure of the exponent. In this case the total effect of subleading logs, N3LL and beyond, is similar to the NNLL contribution (see Fig. 4). As a result the central value of $`\alpha _s(M_Z^2)`$ changes from the pure NNLL result of $`0.1155`$ to $`0.1135`$.

A priori, higher twist is expected to be important for $`F_2`$ at large $`x`$. However, as it stands, there is no experimental evidence for a non-vanishing higher-twist contribution going beyond the resummed perturbation theory. This is consistent with (at least, it does not contradict) the assumption of  that the dominant higher twist at large-$`x`$ is the one related to mixing with the leading twist. On the other hand, in the experimentally accessible range and excluding the very small $`Q^2/N`$ region ($`Q^2/N\mathrm{\Lambda }^2`$), the natural size of the higher twist under the above assumption is still significantly smaller than the theoretical uncertainty.

In addition to solving the theoretical problems mentioned above, progress in this field requires additional data. Although the SLAC and BCDMS data are reasonably precise it is desirable to extend the range they cover. Ideally, to verify or falsify Eq. (1) one would need data covering a large range of $`N`$ and $`Q^2`$ with a fixed $`Q^2/N`$.

Acknowledgements

We would like to thank Volodya Braun, Stefano Catani, Andreas Freund, Gregory Korchemsky, Douglas Ross and Bryan Webber for useful discussions. EG thanks the DFG for financial support. RGR thanks the Leverhulme Trust for the award of an Emeritus Fellowship at RAL.