text:
# KEK preprint 2001-2 Belle preprint 2001-3 Measurement of Inclusive Production of Neutral Pions from 𝚼(𝟒𝑺) Decays Belle Collaboration

Belle Collaboration

(March 24, 2001)

## Abstract

Using the Belle detector operating at the KEKB $`e^+e^{}`$ storage ring, we have measured the mean multiplicity and the momentum spectrum of neutral pions from the decays of the $`\mathrm{{\rm Y}}(4S)`$ resonance. We measure a mean of $`4.70\pm 0.04\pm 0.22`$ neutral pions per $`\mathrm{{\rm Y}}(4S)`$ decay.

This analysis presents the first results on the average multiplicity of neutral pions and their momentum spectrum in $`B`$ meson decays at the $`\mathrm{{\rm Y}}(4S)`$ resonance using the Belle detector ref-belle-detector  at KEKB ref-KEKB . The measurement is based on the sample of $`3.4\times 10^6`$ $`B`$ meson pairs collected by the Belle detector at a center-of-mass (CM) energy of $`\sqrt{s}=10.58\mathrm{GeV}`$.

The particle composition of hadronic final states in $`e^+e^{}`$ annihilation and the measurement of inclusive particle production rates have been important subjects for various energy regions. Previous studies of particle composition at the $`\mathrm{{\rm Y}}(4S)`$ resonance include measurements of charged particles ($`\pi ^\pm `$ ref-CLEO-pipm ; ref-ARGUS-pipm  and $`K^\pm `$ ref-CLEO-kpm ; ref-ARGUS-kpm ), $`\eta `$ mesons ref-CLEO-eta , and vector mesons ref-CLEO-phi ; ref-ARGUS-vect . For typical $`B`$ meson decays, the bulk of the neutral energy is carried by neutral pions. Measurements of inclusive spectra contain information on $`B`$ meson decay mechanisms, especially in the high momentum region where important rare decays may become detectable. If the details of these inclusive spectra are known more precisely, they will allow better estimation of backgrounds and better modeling of $`B`$ meson decay.

The Belle detector (see Fig. 1), described in detail in Ref. ref-belle-detector , consists of a silicon vertex detector (SVD) ref-belle-svd , central drift chamber (CDC) ref-belle-cdc , aerogel Čerenkov counter (ACC) ref-belle-acc , time of flight/trigger scintillation counter (TOF/TSC) ref-belle-tof , CsI electromagnetic calorimeter (ECL) ref-belle-ecl  and KL/muon detector (KLM) ref-belle-klm . The SVD measures the precise position of decay vertices. It consists of three layers of double-sided silicon strip detectors (DSSD) in a barrel-only design and covers 86$`\%`$ of the solid angle. The layer radii are 3.0, 4.5, and 6.0 cm. Charged tracks are reconstructed primarily by the CDC that covers the 17$`{}_{}{}^{}<\theta <`$150 polar angular region. It consists of 50 cylindrical layers of drift cells organized into 11 super-layers (axial or small-angle-stereo), each containing between three and six layers. A low $`Z`$ gas mixture (50% He, 50% C2H6) is used to minimize multiple-Coulomb scattering. The inner and outer radii of the CDC are 9 cm and 86 cm, respectively. The solenoidal magnetic field of 1.5 Tesla is chosen to optimize momentum resolution without sacrificing reconstruction efficiency for low momentum tracks. Kaon identification (KID) is provided by specific ionization ($`dE`$/$`dx`$) measurements in the CDC, Čerenkov threshold measurements in the ACC, and the cylindrical TOF scintillator barrel. The ECL is made of finely segmented CsI(T$`\mathrm{}`$) crystals 30 cm in length. The cross section of one counter is approximately $`55\times 55`$ mm2 at the front surface. The ECL crystals cover the 12$`{}_{}{}^{}<\theta <157^{}`$ angular region. The inner radius of the barrel is 1.25 m, while the annular endcaps are placed at +2.0 m and $`1.0`$ m along the beam line from the interaction point. The calibration of the calorimeter is performed using cosmic rays and Bhabha events. The KLM consists of alternating layers of glass resistive plate counters and 4.7 cm thick iron plates.

The data samples used in this analysis correspond to $`3.2\mathrm{fb}^1`$ of integrated luminosity taken at the $`\mathrm{{\rm Y}}(4S)`$ resonance and $`0.6\mathrm{fb}^1`$ taken at a CM energy 60 MeV below the resonance; the latter was used to subtract underlying continuum background. The integrated luminosity was determined from the number of Bhabha events for which we require both electron and positron in the region of $`46.7^{}<\theta ^{}<145.7^{}`$ in the center-of-mass frame.

Hadronic events are selected based on charged track information from the CDC and cluster information from the ECL. We require at least three charged tracks, that the energy sum in the calorimeter be between 10% and 80% of $`\sqrt{s}`$, and that the charged track momentum be balanced in the $`z`$ direction. This removes the majority of two photon, radiative Bhabha, and $`\tau ^+\tau ^{}`$ events where both $`\tau `$’s decay to leptons. Radiative Bhabha events with one electron outside of the ECL acceptance are removed by requiring at least one large-angle cluster in the ECL and requiring that the average cluster energy be below 1 GeV. Higher multiplicity $`\tau ^+\tau ^{}`$ events are removed if the charged and neutral energy sums in the event are consistent with a $`\tau `$ pair event and if the reconstructed invariant mass of the particles found in each of the two hemispheres perpendicular to the event thrust axis falls below the $`\tau `$ mass. Beam-gas and beam-wall backgrounds are removed by reconstructing the primary vertex of the event and requiring it to be consistent with the known location of the interaction point. The selection is 99% efficient for $`B\overline{B}`$ events and approximately 87% efficient for continuum events. The efficiency for the hadronic event selection was found by a GEANT-based Brun:1987ma  Monte Carlo simulation program. To suppress continuum, we require that the ratio $`R_2`$ of second to zeroth Fox-Wolfram moments ref-fwm , determined using charged tracks and neutral clusters, be less than 0.5.

Photons are reconstructed from neutral clusters in the ECL that have a lateral shape consistent with that of an electromagnetic shower. The energy resolution was measured to be $`\sigma _E/E=0.066\%/E0.81\%/E^{0.25}1.34\%`$ ($`E`$ in GeV) from beam tests ref-belle-ecl . To keep the combinatorial background at a reasonable level, only photons in the central barrel region (35$`{}_{}{}^{}<\theta <`$120) with $`E_\gamma 30\mathrm{MeV}`$ are used in this analysis; the endcap regions have worse energy resolution due to more intervening material and higher beam-associated background.

For each 100 MeV/$`c`$ momentum bin in the CM momentum range 0 to 3 GeV/$`c`$, the $`\gamma \gamma `$ invariant mass distribution is fit to an asymmetric (symmetric above 2 GeV/$`c`$) Gaussian (i.e., a Gaussian with different widths on either side of the mean) for the signal plus a polynomial for the combinatorial background to extract the $`\pi ^0`$ yields. Fig. 2 shows typical mass spectra obtained from the on-resonance data. For the asymmetric Gaussians, the mass resolution is defined as the mean of the left- and right-hand sigmas. An average mass resolution of 5 MeV/$`c^2`$ is obtained, dominated by energy resolution at low momenta or by angular resolution at high momenta. The mass peak is shifted slightly from the established value ref-pdg  because of the asymmetric energy response of the calorimeter due to shower leakage. The observed mass peak position and resolution are consistent with Monte Carlo expectations, as shown in Fig. 3.

To extract the $`\pi ^0`$ momentum spectrum from $`\mathrm{{\rm Y}}(4S)`$ decays, the underlying continuum in the on-resonance data is subtracted bin-by-bin using off-resonance data. The inclusive spectrum is calculated using

$$\frac{1}{\sigma _h}\frac{d\sigma _{\pi ^0}}{dp_{\pi ^0}^{}}=\frac{1}{N_h}\frac{Y_{\mathrm{on}}^i\alpha Y_{\mathrm{off}}^i}{ϵ^i\mathrm{\Delta }p_{\pi ^0}^{}},$$
(1)
where $`N_h`$ is the number of produced $`B\overline{B}`$ events, $`Y_{\mathrm{on}}`$ and $`Y_{\mathrm{off}}`$ are the background-subtracted $`\pi ^0`$ yields obtained from on- and off-resonance data fits, $`\alpha =(_{\mathrm{on}}/_{\mathrm{off}})(s_{\mathrm{off}}/s_{\mathrm{on}})`$ is the on-off scaling factor, and $`ϵ`$ is the product of acceptance and detection efficiency for each momentum bin. The average $`\pi ^0`$ multiplicity is obtained by summing the data from the measured individual momentum bins as shown in Table 1.

The $`\pi ^0`$ acceptance and detection efficiency are determined from Monte Carlo simulations of $`B\overline{B}`$ decays (equal proportions of charged and neutral $`B`$ mesons) and continuum processes. The $`\gamma \gamma `$ invariant mass distributions are fit with the same functions as used in the real data analysis. The product of acceptance and detection efficiency are defined as the ratio of the fitted $`\pi ^0`$ yield to the generated count. Efficiencies for the high momentum $`\pi ^0`$’s above the kinematic limit for $`B\overline{B}`$ decays are deduced from continuum Monte Carlo normalized in the 2.2–2.3 GeV/$`c`$ bin, and scaled efficiencies are used that account for the different acceptances.

Possible sources of systematic uncertainties and their effect on the inclusive $`\pi ^0`$ mean multiplicity measurement are summarized in Table 2. These are discussed next.

The uncertainties in the number of $`B\overline{B}`$ events and the hadronic event selection efficiency are estimated to be 1% and 1.1%, respectively, or 1.5% combined.

The effect of uncertainty in the relative luminosity ratio, used to subtract the continuum background from the on-resonance data, was studied by varying the size of the continuum subtraction by 1%, and we find a 1.5% uncertainty in the mean multiplicity.

Due to the CM energy difference between on- and off-resonance data, the characteristics of continuum events may not match, leading to a bias in the continuum subtraction from the on-resonance data. In particular, the event shape at lower CM energy is slightly less jet-like so that more continuum events survive the common $`R_2`$ cut (resulting in an oversubtraction from the on-resonance data), and the particle momenta scale with CM energy. Hadronic event selection and $`R_2`$ cut bias from the on-off energy difference were determined by using a continuum Monte Carlo Sjostrand:1992mh  sample generated at the CM energy of the off-resonance data; we estimate a 0.4% uncertainty in the $`\pi ^0`$ multiplicity due to this effect. Particle momentum scaling with CM energy was studied by comparing the multiplicity with and without momentum scaling; this leads to an uncertainty of 0.3%.

The uncertainties in the $`\pi ^0`$ detection efficiency due to the application of the shower transverse shape cut, the charged-track veto, photon energy smearing, and the ECL’s modest non-linear energy response correction, are estimated to be 0.4%, 1.6%, 2%, and 0.6%, respectively (determined by comparing the yields with and without each of these criteria).

The decay angular distribution for the photons in the $`\pi ^0`$ rest frame is expected to be isotropic, but the detected distribution of the decay angle $`\theta _d`$— defined as the angle between the photon momentum in the $`\pi ^0`$ rest frame and the $`\pi ^0`$ momentum in the laboratory frame— shows an energy dependence due to the $`\gamma `$ energy cut. Inaccuracy in the Monte Carlo simulation could manifest itself as an anisotropy in the efficiency-corrected decay angle distribution. A comparison between Monte Carlo and data for $`|\mathrm{cos}\theta _d|>0.5`$ and $`|\mathrm{cos}\theta _d|0.5`$ shows no discrepancy and is used to estimate a systematic uncertainty of 0.3%.

The $`\gamma \gamma `$ invariant mass distribution for each CM momentum bin is fit to the following functions in the mass window 80–180 MeV/$`c^2`$:

* $`p_{\gamma \gamma }^{}<1.0`$ GeV/$`c`$: 4th order polynomial plus asymmetric Gaussian;
* $`1.0p_{\gamma \gamma }^{}<2.0`$ GeV/$`c`$: 2nd order polynomial plus asymmetric Gaussian;
* $`2.0p_{\gamma \gamma }^{}<3.0`$ GeV/$`c`$: 2nd order polynomial plus Gaussian.

The fitting range and the order of the polynomial are chosen to minimize statistical fluctuation. By changing the order of the polynomial, the relative yield variation between data and Monte Carlo is taken as the systematic uncertainty in the background modeling; a 2.8% uncertainty was deduced. The weighted average of the fit errors in the efficiency estimation, $`(\delta ϵ)^2=_i\sigma _{ϵ,i}^2Y_i^{\pi ^0}/_iY_i^{\pi ^0}`$, leads to 1% uncertainty in the $`\pi ^0`$ multiplicity. Combining these two numbers gives a 3% systematic uncertainty in the $`\pi ^0`$ mean multiplicity measurement due to the fit procedure.

The measured $`\pi ^0`$ momentum spectrum is compared to the $`B\overline{B}`$ Monte Carlo ref-cleo-qq  prediction in Fig. 4. Systematic differences in the intermediate momentum region may be caused by overestimation of $`bc`$ processes in the Monte Carlo event generation: 100% $`bc`$ is assumed for generic $`B\overline{B}`$ decays in our Monte Carlo. In the high momentum region above the $`bc`$ kinematic end point, we searched for a $`\pi ^0`$ excess from charmless $`B`$ decays. We find $`N_{\mathrm{excess}}=410\pm 724`$ events in the momentum interval 2.4–2.7$`\mathrm{GeV}/c`$, which corresponds to a partial branching fraction $`(B\pi ^0X;p>2.4\mathrm{GeV}/c)=(2\pm 3.5\pm 6.7)\times 10^4`$, or less than $`1.4\times 10^3`$ at 90% confidence level. However, the statistical precision here is limited and the fluctuations are dominated by the off-resonance data used in the subtraction.

In conclusion, using 3.2 $`\mathrm{fb}^1`$ of data accumulated at the $`\mathrm{{\rm Y}}(4S)`$ resonance by the Belle detector, we have measured the inclusive spectrum of neutral pions from $`\mathrm{{\rm Y}}(4S)`$ decays. By summing the measured momentum bins, the mean multiplicity of neutral pions from $`\mathrm{{\rm Y}}(4S)`$ decays has been determined to be $`n_{\pi ^0}=4.70\pm 0.04\pm 0.22`$, corresponding to an inclusive branching fraction111 Inclusive branching fractions in $`B`$ decays have a definition in terms of multiplicity and can be greater than 100% ref-pdg . of $`(B\pi ^0X)=(235\pm 2\pm 11)`$%, where the first error is statistical and the second is systematic.

###### Acknowledgements.

We gratefully acknowledge the efforts of the KEKB group in providing us with excellent luminosity and running conditions and the members of the KEK computing research center for their assistance with our computing and network systems. We thank the staffs of KEK and collaborating institutions for their contributions to this work, and acknowledge support from the Ministry of Education, Science, Sports and Culture of Japan and the Japan Society for the Promotion of Science; the Australian Research Council and the Australian Department of Industry, Science and Resources; the Department of Science and Technology of India; the BK21 program of the Ministry of Education of Korea and the Basic Science program of the Korea Science and Engineering Foundation; the Polish State Committee for Scientific Research under contract No.2P03B 17017; the Ministry of Science and Technology of Russian Federation; the National Science Council and the Ministry of Education of Taiwan; the Japan-Taiwan Cooperative Program of the Interchange Association; and the U.S. Department of Energy.


format: text

text:
# 1 Introduction

## 1 Introduction

Structure functions in inclusive deep-inelastic lepton-nucleon scattering (DIS) are among the observables best suited for precise determinations of the strong coupling constant $`\alpha _s`$. At present their experimental uncertainties result in an error $`\mathrm{\Delta }_{\mathrm{exp}}\alpha _s(M_Z^2)0.002`$ at the mass of the $`Z`$-boson . A further reduction of this error can be expected, especially from measurements at the electron-proton collider HERA after the forthcoming luminosity upgrade. The standard next-to-leading order (NLO) approximation of perturbative QCD summarized in ref. , on the other hand, leads to a theoretical error $`\mathrm{\Delta }_{\mathrm{th}}\alpha _s(M_Z^2)0.005`$. This error is dominated by the uncertainty due to the truncation of the perturbation series as estimated from the renormalization-scale dependence. Hence calculations beyond NLO are required to make full use of the present and forthcoming data on structure functions.

The ingredients necessary for next-to-next-to-leading order (NNLO) analyses of the structure functions in Bjorken-$`x`$ space111See refs.  for NNLO analyses based on the integer Mellin-$`N`$ results of refs.  only. have not been completed up to now: Unlike the two-loop coefficient functions which were calculated some time ago  (and completely checked recently ), only partial results  have been obtained for the three-loop splitting functions so far. However, we have recently demonstrated  that the uncertainties resulting from the incompleteness of this information are entirely negligible at $`x>0.05`$. Moreover, these uncertainties are small even at much lower $`x`$, down to $`x10^4`$ at not too small scales, $`Q^2\stackrel{>}{}10\text{ GeV}^2`$ . Thus analyses of structure functions in DIS (and of total cross sections for Drell-Yan lepton-pair production ) can be promoted to NNLO over a wide kinematic region. Besides more accurate determinations of the parton densities, such analyses facilitate a considerably improved theoretical accuracy $`\mathrm{\Delta }_{\mathrm{th}}\alpha _s(M_Z^2)0.002`$ of the determinations of the strong coupling.

In the present article we extend, for $`x>10^2`$, our treatment  of the flavour non-singlet (NS) sector dominating $`\alpha _s`$-extractions from DIS to the next-to-next-to-next-to-leading order (N3LO). This extension is facilitated by two circumstances: The first is the existence of constraints on the three-loop coefficient functions which prove to be sufficiently restrictive in this region of $`x`$. The seven lowest even-integer and odd-integer moments have been computed  for the structure function $`F_{2,\mathrm{NS}}`$ in electromagnetic DIS and $`F_3^{\nu +\overline{\nu }}`$ in charged-current DIS, respectively. Furthermore the four leading large-$`x`$ terms of these functions are known from the soft-gluon resummation . The second circumstance is the rapid convergence of the splitting-function expansion in the usual $`\overline{\text{MS}}`$ factorization scheme also employed in refs. . Already the impact of the three-loop splitting functions is small at $`x>10^2`$, in absolute size (less than 1% on $`\alpha _s(M_Z^2)`$) as well as compared to the two-loop coefficient functions . Hence one can safely expect that the effect of the unknown four-loop splitting functions on determinations of $`\alpha _s(M_Z^2)`$ will be well below the 1% accuracy we are aiming at.

As demonstrated below, the N3LO approximation suffices for achieving this accuracy in the region $`x\stackrel{<}{}\mathrm{\hspace{0.17em}0.75}`$ usually covered by analyses of DIS data . Terms beyond this order are relevant at $`x\stackrel{>}{}\mathrm{\hspace{0.17em}0.8}`$, on the other hand, mainly due to the presence of large soft-gluon logarithms up to $`[\mathrm{ln}^{2l1}(1x)]/(1x)`$ in the $`l`$-loop coefficient functions. The resummation of these logarithms  has been extended to the next-to-next-to-leading logarithmic (NNLL) accuracy recently . Here we will study the predictions of this resummation for the factorization-scheme independent (‘physical’) kernels governing the scaling violations (‘evolution’) of the non-singlet structure functions. Other approaches to estimate higher-order corrections to these kernels, not restricted to very large $`x`$, include Padé summations of the perturbation series  as well as renormalization scheme optimizations such as the principle of minimal sensitivity (PMS)  and the effective charge (ECH) method . We will also compare these estimates to the full NNLO and N3LO evolution kernels, and investigate the resulting predictions at order $`\alpha _s^5`$ (N4LO) and beyond.

The outline of this article is as follows: In section 2 we express the physical evolution kernels, up to N4LO and NNLL accuracy, in terms of the corresponding splitting functions and coefficient functions. The information on the three-loop coefficient functions for $`F_{2,\mathrm{NS}}`$ and $`F_3`$ discussed above is employed in section 3 to derive approximate expressions for their $`x`$-dependence. Besides these functions the N3LO evolution kernels also involve convolutions of lower-order coefficient functions for which we provide compact expressions in section 4. These results are put together in section 5 to study the effects of the N3LO terms on the evolution of the structure functions and on the resulting determinations of the strong coupling constant. In section 6 we discuss the predictions of the soft-gluon resummation and of the Padé, PMS and ECH approximations. Finally our results are summarized in section 7. Some relations for the convolutions in section 4 and for the Padé approximations in section 6 can be found in the appendix.

## 2 Fixed-order and resummed evolution kernels

For the choice $`\mu _r^2=\mu _f^2=Q^2`$ of the renormalization and mass-factorization scales, the structure functions

$$_1=2F_{1,\mathrm{NS}},_2=\frac{1}{x}F_{2,\mathrm{NS}},_3=F_3^{\nu \pm \overline{\nu }}$$
(2.1)
are in perturbative QCD given by

$`_a(x,Q^2)`$ $`=`$ $`\left[C_a(Q^2)q_{a,\mathrm{NS}}^{}(Q^2)\right](x)`$ (2.2)
$`=`$ $`{\displaystyle \underset{l=0}{}}a_s^l(Q^2)\left[c_{a,l}q_{a,\mathrm{NS}}^{}(Q^2)\right](x)`$
$`=`$ $`{\displaystyle \underset{l=0}{}}a_s^l(Q^2){\displaystyle _x^1}{\displaystyle \frac{dy}{y}}c_{a,l}(y)q_{a,\mathrm{NS}}({\displaystyle \frac{x}{y}},Q^2).`$
Here $`c_{a,l}`$ represents the $`l`$-loop non-singlet coefficient functions with $`c_{a,0}(x)=\delta (1x)`$, and $`q_{a,\mathrm{NS}}`$ stands for the respective combinations of the quark densities. The scale dependence of the running coupling of QCD, in this article normalized as

$`a_s{\displaystyle \frac{\alpha _s}{4\pi }},`$ (2.3)
is governed by

$$\frac{da_s}{d\mathrm{ln}Q^2}=\beta (a_s)=\underset{l=0}{}a_s^{l+2}\beta _l.$$
(2.4)
Besides $`\beta _0`$ and $`\beta _1`$  also the coefficients $`\beta _2`$ and $`\beta _3`$ have been computed  in the $`\overline{\text{MS}}`$ renormalization scheme adopted throughout this study. All these four coefficients,

$`\beta _0`$ $`=`$ $`\mathrm{\hspace{0.25em}11}\mathrm{\hspace{0.25em}2}/3N_f`$
$`\beta _1`$ $`=`$ $`10238/3N_f`$
$`\beta _2`$ $`=`$ $`\mathrm{\hspace{0.17em}2857}/25033/18N_f+\mathrm{\hspace{0.17em}325}/54N_f^2`$
$`\beta _3`$ $`=`$ $`29243.0\mathrm{\hspace{0.25em}6946.30}N_f+405.089N_f^2+1093/729N_f^3,`$ (2.5)
are required for N3LO calculations. The irrational coefficients in Eq. (2) and in Eqs. (2.15) and (2.18) below have been truncated to six digits for brevity. $`N_f`$ denotes the number of effectively massless flavours (mass effects are not considered in this article). Finally the evolution equations for the parton densities in Eq. (2.2) read

$`{\displaystyle \frac{d}{d\mathrm{ln}Q^2}}q_{a,\mathrm{NS}}^{}(x,Q^2)`$ $`=`$ $`\left[P_a(Q^2)q_{a,\mathrm{NS}}^{}(Q^2)\right](x)`$ (2.6)
$`=`$ $`{\displaystyle \underset{l=0}{}}a_s^{l+1}\left[P_{a,l}q_{a,\mathrm{NS}}^{}(Q^2)\right](x),`$
where $``$ abbreviates the Mellin convolution written out in the third line of Eq. (2.2). Like the coefficient functions $`c_{a,l}(x)`$, the $`(l+1)`$-loop splitting functions $`P_{a,l}(x)`$ are scale independent for the above choice of $`\mu _r`$ and $`\mu _f`$.

Explicit expressions up to order $`\alpha _s^2`$ can be found in refs.  and  for the non-singlet splitting functions and coefficient functions, respectively. For the third-order splitting functions $`P_{a,2}(x)`$ we will employ our approximate expressions of ref. . The three-loop coefficient functions $`c_{a,3}(x)`$ are the subject of section 3 below.

It is convenient to express the scaling violations of the non-singlet structure functions in terms of these structure functions themselves, thus explicitly eliminating any dependence on the factorization scheme and the scale $`\mu _f`$. The corresponding ‘physical’ evolution kernels $`K_a`$ for $`\mu _r^2=Q^2`$ can be derived by differentiating Eq. (2.2) with respect to $`Q^2`$ by means of Eqs. (2.4) and (2.6), and finally eliminating $`q_{a,\mathrm{NS}}^{}`$ using the inverse of Eq. (2.2). Suppressing the dependences on $`x`$ and $`Q^2`$ one arrives at the evolution equations

$`{\displaystyle \frac{d}{d\mathrm{ln}Q^2}}_a`$ $`=`$ $`\left(P_a+\beta {\displaystyle \frac{d\mathrm{ln}C_a}{da_s}}\right)_a`$ (2.7)
$`=`$ $`K_a_a={\displaystyle \underset{l=0}{}}a_s^{l+1}K_{a,l}_a`$
$`=`$ $`\left\{a_sP_0+{\displaystyle \underset{l=1}{}}a_s^{l+1}\left(P_{a,l}{\displaystyle \underset{k=0}{\overset{l1}{}}}\beta _k\stackrel{~}{c}_{a,lk}\right)\right\}_a`$
with

$`\stackrel{~}{c}_{a,1}`$ $`=`$ $`c_{a,1}`$
$`\stackrel{~}{c}_{a,2}`$ $`=`$ $`2c_{a,2}c_{a,1}^2`$
$`\stackrel{~}{c}_{a,3}`$ $`=`$ $`3c_{a,3}3c_{a,2}c_{a,1}+c_{a,1}^3`$ (2.8)
$`\stackrel{~}{c}_{a,4}`$ $`=`$ $`4c_{a,4}4c_{a,3}c_{a,1}2c_{a,2}^2+4c_{a,2}c_{a,1}^2c_{a,1}^4.`$
In Eq. (LABEL:eq26) we have used the abbreviation $`f^l`$ for the $`(l1)`$-fold convolution of a function $`f(x)`$ with itself, i.e., $`f^2=ff`$ etc. The generalizations $`𝒦_{a,l}`$ of the kernels $`K_{a,l}`$ in Eq. (2.7) to $`\mu _r^2Q^2`$ can be obtained by expanding 222Up to the fifth order this expansion can be read off from the $`K_{a,0}`$ terms of Eq. (2). $`a_s(Q^2)`$ in terms of $`a_s(\mu _r^2)`$ and $`L=\mathrm{ln}(Q^2/\mu _r^2)`$, yielding

$`𝒦_{a,0}`$ $`=`$ $`K_{a,0}`$
$`𝒦_{a,1}`$ $`=`$ $`K_{a,1}\beta _0LK_{a,0}`$
$`𝒦_{a,2}`$ $`=`$ $`K_{a,2}2\beta _0LK_{a,1}(\beta _1L\beta _0^2L^2)K_{a,0}`$
$`𝒦_{a,3}`$ $`=`$ $`K_{a,3}3\beta _0LK_{a,2}(2\beta _1L3\beta _0^2L^2)K_{a,1}`$
$`\left(\beta _2L{\displaystyle \frac{5}{2}}\beta _1\beta _0L^2+\beta _0^3L^3\right)K_{a,0}`$
$`𝒦_{a,4}`$ $`=`$ $`K_{a,4}4\beta _0LK_{a,3}(3\beta _1L6\beta _0^2L^2)K_{a,2}`$
$`(2\beta _2L7\beta _1\beta _0L^2+4\beta _0^3L^3)K_{a,1}`$
$`\left(\beta _3L3\beta _2\beta _0L^2{\displaystyle \frac{3}{2}}\beta _1^2L^2+{\displaystyle \frac{13}{3}}\beta _1\beta _0^2L^3\beta _0^4L^4\right)K_{a,0}.`$
At N$`{}_{}{}^{m}{}_{}{}^{}`$LO the terms up to $`l=m`$ are included in Eq (2.7). Hence Eqs. (LABEL:eq26) and (2) formally specify the evolution kernels $`𝒦_a`$ up to N$`{}_{}{}^{4}{}_{}{}^{}`$LO. Their extension to higher orders is straightforward but irrelevant for the time being, as at least the coefficient functions beyond four loops will not be calculated in the foreseeable future.

The leading terms of the coefficient functions for $`x1`$, however, are known to all orders from the soft-gluon resummation . Switching to Mellin moments,

$$f^N=_0^1𝑑xx^{N1}f(x),$$
(2.10)
for the remainder of this section, the large-$`N`$ (large-$`x`$) behaviour of the coefficient functions in Eq. (2.2) takes the form

$$C_{\mathrm{res}}^N=g_0(a_s)\mathrm{exp}\left\{\mathrm{ln}Ng_1(\lambda )+g_2(\lambda )+a_sg_3(\lambda )+𝒪(a_s^2f(\lambda ))\right\}$$
(2.11)
up to terms which vanish for $`N\mathrm{}`$. Here we have used the abbreviation

$$\lambda =a_s\beta _0\mathrm{ln}N=\frac{\alpha _s}{4\pi }\beta _0\mathrm{ln}N,$$
(2.12)
and we have again put $`\mu _r^2=\mu _f^2=Q^2`$. By virtue of the first line of (2.7), Eq. (2.11) leads to the following expression for the resummed kernel up to next-to-next-to-leading logarithmic (NNLL) accuracy :

$`K_{\mathrm{res}}^N`$ $`=`$ $`(A_1a_s+A_2a_s^2+A_3a_s^3)\mathrm{ln}N\left(1+{\displaystyle \frac{\beta _1}{\beta _0}}a_s+{\displaystyle \frac{\beta _2}{\beta _0}}a_s^2\right)\lambda ^2{\displaystyle \frac{dg_1}{d\lambda }}`$ (2.13)
$`(a_s\beta _0+a_s^2\beta _1)\lambda {\displaystyle \frac{dg_2}{d\lambda }}a_s^2\beta _0{\displaystyle \frac{d}{d\lambda }}\left(\lambda g_3(\lambda )\right)+𝒪(a_s^3(f(\lambda )).`$
Thus the leading logarithmic (LL), next-to-leading logarithmic (NLL) and NNLL large-$`N`$ contributions to the physical evolution kernels are of the form $`(a_s\mathrm{ln}N)^n`$, $`a_s(a_s\mathrm{ln}N)^n`$ and $`a_s^2(a_s\mathrm{ln}N)^n`$, respectively. This is in contrast to the coefficient functions which receive contributions up to $`(a_s\mathrm{ln}^2N)^n`$. The constants $`A_l`$ in Eq. (2.13) are the coefficients of the leading  large-$`x`$ terms $`1/[1x]_+`$ of the $`l`$-loop $`\overline{\text{MS}}`$ splitting functions — recall that

$$f^N=\mathrm{ln}N+𝒪(1)\text{ for }f(x)=\mathrm{\hspace{0.25em}1}/[1x]_+.$$
(2.14)
As in Eq. (2) inserting the numerical values for the QCD colour factors $`C_A`$ and $`C_F`$, these constants are given by

$$A_1=\mathrm{\hspace{0.25em}16}/3,A_2=\mathrm{\hspace{0.25em}66.4732}160/27N_f,$$
(2.15)
and the yet approximate, but sufficiently accurate three-loop result

$$A_3=(1178.8\pm 11.5)(183.95\pm 0.85)N_f64/81N_f^2.$$
(2.16)
Inserting the explicit form of the functions $`g_1`$, $`g_2`$  and $`g_3`$  in Eq. (2.13) and restoring the dependence on $`L=\mathrm{ln}(Q^2/\mu _r^2)`$ leads to

$`𝒦_{\mathrm{res}}^N`$ $`=`$ $`{\displaystyle \frac{A_1}{\beta _0}}\mathrm{ln}(1\lambda )+a_sLA_1{\displaystyle \frac{\lambda }{1\lambda }}+a_s^2L^2A_1\beta _0{\displaystyle \frac{\lambda (\lambda 2)}{2(1\lambda )^2}}`$
$`+a_s\left\{A_1\beta _1\left(\lambda +\mathrm{ln}(1\lambda )\right)+(B_1A_1\gamma _e)\beta _0^2\lambda A_2\beta _0\lambda \right\}{\displaystyle \frac{1}{\beta _0^2(1\lambda )}}`$
$`+a_s^2L\left\{\left[(B_1\gamma _eA_1)\beta _0A_2\right]{\displaystyle \frac{\lambda (\lambda 2)}{(1\lambda )^2}}A_1\beta _1{\displaystyle \frac{\mathrm{ln}(1\lambda )}{\beta _0(1\lambda )^2}}\right\}`$
$`+a_s^2\{[A_1(\gamma _e^2+\zeta _2)\beta _0^2+2A_2\gamma _e\beta _0+A_32B_1\gamma _e\beta _0^22B_2\beta _0]{\displaystyle \frac{\lambda (\lambda 2)}{2\beta _0(1\lambda )^2}}`$
$`+[\mathrm{\hspace{0.17em}2}(A_1\gamma _eB_1)\beta _1\beta _0^2\mathrm{ln}(1\lambda )+A_1\beta _1^2(\lambda ^2\mathrm{ln}^2(1\lambda ))`$
$`A_1\beta _2\beta _0\lambda ^2A_2\beta _1\beta _0(\lambda ^22\lambda 2\mathrm{ln}(1\lambda ))]{\displaystyle \frac{1}{2\beta _0^3(1\lambda )^2}}`$
$`+2D_2{\displaystyle \frac{\lambda (1\lambda )}{(12\lambda )^2}}\}+𝒪(a_s^3f(\lambda ,L))`$
$``$ $`{\displaystyle \underset{l=0}{\overset{\mathrm{}}{}}}a_s^{l+1}𝒦_{\mathrm{res},l}^N.`$
Here $`\zeta _2=\pi ^2/6`$, and $`\gamma _e`$ represents the Euler-Mascheroni constant, $`\gamma _e0.577216`$. Furthermore $`B_1=4`$ , and the constants $`B_2`$ and $`D_2`$ are related by

$$B_2+D_2=\mathrm{\hspace{0.25em}36.2657}+6.34888N_f.$$
(2.18)
We will return to the latter coefficients at the end of the next section.

After subtracting the terms up to order $`a_s^{m+1}`$ in $`𝒦_{\mathrm{res}}^N`$ already taken into account in the NmLO terms (2), Eq. (2) can be added to these fixed-order results to obtain the (NmLO + resummed) approximation for the non-singlet evolution kernels,

$$𝒦_{\mathrm{N}^m\mathrm{LO}+\mathrm{res}}^N=\underset{l=0}{\overset{m}{}}a_s^{l+1}(𝒦_{a,l}^N𝒦_{\mathrm{res},l}^N)+𝒦_{\mathrm{res}}^N.$$
(2.19)
Due to the renormalon singularities at $`\lambda =1`$ and $`\lambda =1/2`$ in Eq. (2) the resummed evolution equations cannot be uniquely inverted to $`x`$-space, unlike the fixed-order approximations $`𝒦_{a,l}^N_a^N`$. Note that strength of these singularities — located at $`N2000`$ and $`N45`$ for $`\lambda =1`$ and $`\lambda =1/2`$, respectively, at $`\alpha _s=0.2`$ and $`N_f=4`$ — increases with the order of the soft-gluon expansion: the behaviour is logarithmic at the leading-log level, but involves poles of order $`k`$ for the NkLL approximations. For our numerical study of the all-order case at the end of section 6 we will use the standard ‘minimal prescription’ contour  for the Mellin inversion. This contour runs to the left of the renormalon singularities, but to the right of all other poles in the $`N`$-plane.

## 3 The 3-loop non-singlet coefficient functions

The $`l`$-loop coefficient functions $`c_{a,l}`$ for the non-singlet structure functions $`_{a=1,2,3}`$ defined in Eq. (2.1) can be represented as

$`c_{a,l}(x,N_f)`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{2l1}{}}}A_l^{(m)}𝒟_m+\stackrel{~}{B}_l\delta (1x)+c_{a,l}^{\mathrm{smooth}}(x,N_f)`$ (3.1)
$`+`$ $`{\displaystyle \underset{n=1}{\overset{2l1}{}}}\left(C_{a,l}^{(n)}L_1^n+D_{a,l}^{(n)}L_0^n\right).`$
Here $`A_l^{(m)}`$, $`\stackrel{~}{B}_l`$, $`C_{a,l}^{(n)}`$ and $`D_{a,l}^{(n)}`$ are numerical coefficients which in general depend on the number of flavours $`N_f`$, and we have employed the abbreviations

$$𝒟_k=\left[\frac{\mathrm{ln}^k(1x)}{1x}\right]_+,L_1=\mathrm{ln}(1x),L_0=\mathrm{ln}x$$
(3.2)
for the +-distributions (see Eqs. (3.3) and (3.4) below) and the end-point logarithms. The functions $`c_{a,l}^{\mathrm{smooth}}(x,N_f)`$ collect all contributions which are finite for $`0x1`$. This regular term constitutes the mathematically complicated part of Eq. (3.1), it involves higher transcendental functions like the harmonic polylogarithms introduced in ref. . As usual, the +-distributions are defined via

$$_0^1𝑑xa(x)_+f(x)=_0^1𝑑xa(x)\{f(x)f(1)\}$$
(3.3)
where $`f(x)`$ is a regular function. The convolutions with the distributions occurring in Eq. (3.1) can be written as 333The second line of Eq. (3.4) is given by $`xf(x)_0^x𝑑ya(y)`$ for a general +-distribution $`[a(x)]_+`$.

$`x[𝒟_kf](x)`$ $`=`$ $`{\displaystyle _x^1}𝑑y{\displaystyle \frac{\mathrm{ln}^k(1x)}{1x}}\left\{{\displaystyle \frac{x}{y}}f\left({\displaystyle \frac{x}{y}}\right)xf(x)\right\}`$ (3.4)
$`+xf(x){\displaystyle \frac{1}{k+1}}\mathrm{ln}^{k+1}(1x).`$
As already indicated in Eq. (3.1), the coefficients of $`𝒟_m`$ and of $`\delta (1x)`$ are independent of the choice of the structure function.

The three-loop contributions $`c_{S,3}`$ known from the soft-gluon resummation read

$`c_{S,3}(x,3)`$ $`=`$ $`512/27(𝒟_5L_1^5)14400/81𝒟_4+264.062𝒟_3+1781.704𝒟_2`$
$`c_{S,3}(x,4)`$ $`=`$ $`512/27(𝒟_5L_1^5)13760/81𝒟_4+188.210𝒟_3+1962.178𝒟_2`$ (3.5)
$`c_{S,3}(x,5)`$ $`=`$ $`512/27(𝒟_5L_1^5)13120/81𝒟_4+113.938𝒟_3+2131.195𝒟_2,`$
where we have again truncated the irrational coefficients and restricted ourselves to the practically relevant cases $`N_f=3`$, 4 and 5. Besides the $`𝒟_m`$-terms determined in ref. , Eq. (3) also includes the leading integrable large-$`x`$ logarithm. The general relation between the coefficients of this term and the leading +-distribution has been conjectured in ref. . Eqs. (3) complement the main present constraints on $`c_{a,3}(x,N_f)`$ provided by the computation  of the seven lowest even-integer and odd-integer moments (2.10), respectively, for electromagnetic (e.m.) DIS and the charged-current (CC) combination $`F_3^{\nu +\overline{\nu }}`$. Note that the coefficients of the leading small-$`x`$ logarithms are presently unknown here, unlike for the splitting functions and the singlet coefficient functions .

We use this information for approximate reconstructions of $`c_{2,3}(x,N_f)`$ and $`c_{3,3}(x,N_f)`$ at $`N_f=3`$, 4 and 5. Our method is analogous to the treatment of the three-loop splitting functions in refs. : A simple ansatz is chosen for $`c_{a,3}^{\mathrm{smooth}}`$ in Eq. (3.1), and its free parameters are determined from the available moments together with a reasonably balanced subset of the coefficients $`A_3^{(0,1)}`$, $`\stackrel{~}{B}_3`$, $`C_{a,3}^{(n)}`$ and $`D_{a,3}^{(n)}`$. The ansatz for $`c_{a,3}^{\mathrm{smooth}}`$ and the choice of the non-vanishing end-point parameters are then varied in order to estimate the residual uncertainty of $`c_{a,3}`$. Specifically we keep $`A_3^{(1)}`$; one of each pair $`A_3^{(0)}`$ and $`\stackrel{~}{B}_3`$, $`C_{a,3}^{(4)}`$ and $`C_{a,3}^{(3)}`$, $`C_{a,3}^{(2)}`$ and $`C_{a,3}^{(1)}`$; one or two of the $`D_{a,3}^{(n<4)}`$; and one or two parameters of a polynomial up to second order in $`x`$ representing $`c_{a,3}^{\mathrm{smooth}}`$. For a few combinations the resulting system of linear equations which fixes these parameters by the seven moments becomes almost singular, resulting in exceptionally large numerical coefficients. After rejecting those about 5% of the combinations for which the modulus of at least one parameter exceeds $`10^5`$, we are left with about 90 approximations for each case.

Before we present the approximate results for $`c_{a,3}(x,N_f)`$, it is appropriate to illustrate our procedure by applying it to a known quantity, for which we choose the two-loop e.m. coefficient function $`c_{2,2}(x,N_f=4)`$. Adopting the coefficients of $`𝒟_3`$ and $`𝒟_2`$ defined in Eq. (3.2) from the soft-gluon exponentiation, the procedure described in the preceding paragraph is applied to this function with the small adjustment that two of the $`C_{2,2}^{(1,2,3)}`$ are kept as $`C^{(4)}`$ does not occur at two loops according to Eq. (3.1). Also here we reject a couple of combinations, those with parameter(s) of modulus above $`310^3`$. The remaining about 70 approximations are compared to the exact result of ref.  in Fig. 1.

The seven lowest even-integer moments supplemented by the soft-gluon coefficients $`A_2^{(m>1)}`$ prove to constrain $`c_{2,2}(x)`$ rather tightly at $`x\stackrel{>}{}0.3`$. The region $`0.1\stackrel{<}{}x\stackrel{<}{}0.3`$ is less accurately covered, and at $`x<0.1`$ the lack of small-$`x`$ information mentioned above becomes very prominent. Also shown in Fig. 1 is the exact result  for $`c_{2,2}`$ in charged-current DIS. The difference to the electromagnetic case — originating in a sign difference of the contributions from $`\gamma /W+qq+q+\overline{q}`$ with identical quarks in the final-state — is clearly visible only at $`x\stackrel{<}{}0.2`$. The effect of this difference on the evolution of the structure functions at NNLO is unnoticeable at $`x\stackrel{>}{}0.1`$, and amounts to less that 1% for $`x>0.01`$, see Fig. 11 of ref. . We expect that the corresponding three-loop effect will at least not be larger. Hence the approximations for $`c_{2,2}(x,N_f)`$, constructed for the e.m. case, should be applicable also for neutrino DIS without introducing any relevant error.

The corresponding approximations for $`c_{2,3}`$ and $`c_{3,3}`$ are shown in Fig. 2 for $`N_f=4`$ (concerning the scale of the ordinate recall the rather small expansion parameter (2.3) ). As expected, the accuracy pattern is qualitatively similar to the two-loop case of Fig. 2. The uncertainty of $`c_{3,3}`$ is smaller than that of $`c_{2,3}`$ at small $`x`$, since for $`c_{3,3}`$ the lowest calculated moment, $`N=1`$, is closer to the location of the rightmost pole at $`N=0`$. For both functions two representatives, denoted by $`A`$ and $`B`$, are selected which rather completely cover the uncertainty bands. With $`c_{S,3}`$ of Eq. (3) these representatives read

$`c_{2,3}^A(x,4)`$ $`=`$ $`c_{S,3}(x,4)6456.231𝒟_11085.97\delta (1x)+258.876L_1^4`$
$`22430.79L_174705.15x^24062.14(2x)313.0L_0^3`$
$`c_{2,3}^B(x,4)`$ $`=`$ $`c_{S,3}(x,4)5081.227𝒟_1+5028.23𝒟_0+1059.423L_1^4`$
$`7292.84L_1^217741.28(2x)18154.5L_0+1168.02L_0^3`$
and

$`c_{3,3}^A(x,4)`$ $`=`$ $`c_{S,3}(x,4)6940.648𝒟_11526.23𝒟_042.598L_1^4`$
$`33562.64L_191639.87x^25898.84+424.49L_0^2`$
$`c_{3,3}^B(x,4)`$ $`=`$ $`c_{S,3}(x,4)4907.988𝒟_1+5587.906𝒟_0+1092.436L_1^4`$
$`8267.99L_1^218120.787083.63L_0+283.59L_0^3.`$
The uncertainty bands of Fig. 2 do not directly indicate the range of applicability of our approximations, as the coefficient function enter the structure functions and their evolution only via the smoothening convolution (2.2) with non-perturbative initial distributions. In Fig. 3 we therefore present the convolutions of the results (3) and (3) with a typical non-singlet shape. This illustration shows that the residual uncertainties of $`c_{a,3}`$ do not lead to any relevant effects for $`x0.1`$. The situation at smaller $`x`$ depends on the numerical size of the $`c_{a,3}`$ contributions to the evolution kernels given by Eqs. (2.7) and (LABEL:eq26). Anticipating our findings in section 5, we note that these contributions are actually unproblematically small for $`x>10^2`$.

The results for $`N_f=3`$ and $`N_f=5`$ are similar to those presented in Fig. 2 and Fig. 3. For brevity they are not displayed graphically in this article. The selected approximations for $`N_f=3`$ are given by

$`c_{2,3}^A(x,3)`$ $`=`$ $`c_{S,3}(x,3)6973.782𝒟_13929.78\delta (1x)+232.676L_1^4`$
$`35949.75L_193141.53x^210283.51418.43L_0^3`$
$`c_{2,3}^B(x,3)`$ $`=`$ $`c_{S,3}(x,3)5575.903𝒟_1+5474.48𝒟_0+927.478L_1^4`$
$`4646.38L_1^223345.8511094.92L_0+759.69L_0^3`$
and

$`c_{3,3}^A(x,3)`$ $`=`$ $`c_{S,3}(x,3)7018.496𝒟_12957.33\delta (1x)+162.667L_1^4`$
$`28363.91L_174640.87x^25720.48(1+x)+330.21L_0^2`$
$`c_{3,3}^B(x,3)`$ $`=`$ $`c_{S,3}(x,3)5160.335𝒟_1+6896.360𝒟_0+1109.512L_1^4`$
$`7715.43L_1^220541.227595.83L_0+290.34L_0^3.`$
The corresponding functions for $`N_f=5`$ read

$`c_{2,3}^A(x,5)`$ $`=`$ $`c_{S,3}(x,5)5951.174𝒟_1391.37𝒟_02341.422L_1^3`$
$`+19986.58L_1+5517.39(2x)+5969.63L_0284.23L_0^3`$
$`c_{2,3}^B(x,5)`$ $`=`$ $`c_{S,3}(x,5)4802.695𝒟_1+3784.97𝒟_0+1041.041L_1^4`$
$`8021.15L_1^215556.5(2x)16445.21L_0+1084.36L_0^3`$
and

$`c_{3,3}^A(x,5)`$ $`=`$ $`c_{S,3}(x,5)6560.902𝒟_12412.54𝒟_0+98.499L_1^3`$
$`27899.69L_182015.71x^22983.8261.43L_0^3`$
$`c_{3,3}^B(x,5)`$ $`=`$ $`c_{S,3}(x,5)4637.854𝒟_1+4317.29𝒟_0+1070.036L_1^4`$
$`8767.685L_1^215676.576543.88L_0+276.32L_0^3.`$
In all cases the average $`1/2(c_{a,3}^A+c_{a,3}^B)`$ represents our central result.

We conclude this section by returning to the coefficients $`B_2`$ and $`D_2`$ entering the NNLL soft-gluon resummation of the quark coefficient functions (2.11) and structure-function evolution kernels (2). If only one of these constants were present, say $`D_2`$, then this constant would be fixed by the consistency of Eq. (2.11) with the soft-gluon part $`c_{S,2}`$ of the NNLO coefficient functions of ref. , more precisely by the coefficient $`A_2^{(0)}`$ in Eq. (3.1). Digressing for a moment, we note that this situation is actually realized for the (very closely related) NNLL soft-gluon resummations of the quark-antiquark annihilation contribution to the Drell-Yan cross section  and of the cross section for Higgs production via gluon-gluon fusion in the heavy top-quark limit. For these two processes the NNLO soft- and virtual-gluon contribution have been computed in refs.  and , respectively. In the present DIS case, however, this consistency conditions only implies the constraint (2.18). An exact result for the coefficient $`A_2^{(3)}`$ of the three-loop coefficient functions (3.1) would suffice to determine $`B_2`$ and $`D_2`$, as this coefficient is related to the combination $`B_2+2D_2`$ independent from Eq. (2.18).

As discussed in the paragraph below that of Eq. (3), all our approximations (shown for $`N_f=4`$ in Fig. 2) include $`A_2^{(3)}`$. For about 95% of these approximations this coefficient falls into the range

$$A_2^{(3)}\left\{\begin{array}{c}6800\mathrm{}5800\\ 6350\mathrm{}5500\\ 5950\mathrm{}5200\end{array}\right\}\text{ for }N_f=\left\{\begin{array}{c}3\\ 4\\ 5\end{array}\right\}.$$
(3.12)
The comparison of these results to the expansion of Eq. (2.11) (using $`g_3(\lambda )`$ of ref. ) leads to the rather weak constraints

$$B_2\left\{\begin{array}{c}32\mathrm{}\mathrm{\hspace{0.25em}87}\\ 42\mathrm{}\mathrm{\hspace{0.25em}93}\\ 49\mathrm{}\mathrm{\hspace{0.25em}98}\end{array}\right\}\text{ for }N_f=\left\{\begin{array}{c}3\\ 4\\ 5\end{array}\right\}$$
(3.13)
which can to sufficient accuracy be combined to the estimate

$$B_2P_1^\delta +\frac{1}{3}\xi \beta _0P_0^\delta \text{ with }8\stackrel{<}{}\xi \stackrel{<}{}\mathrm{\hspace{0.25em}12}.$$
(3.14)
Here $`P_{l1}^\delta `$ represent the coefficients of $`\delta (1x)`$ in the $`l`$-loop quark splitting functions. Retaining the colour factors $`C_A=N_c=3`$, $`C_F=(N_c^21)/(2N_c)`$ these coefficients read

$`P_0^\delta `$ $`=`$ $`3C_F`$
$`P_1^\delta `$ $`=`$ $`C_F^2\left({\displaystyle \frac{3}{2}}+24\zeta _312\zeta _2\right)+C_FC_A\left({\displaystyle \frac{17}{6}}12\zeta _3+{\displaystyle \frac{44}{3}}\zeta _2\right)`$
$`C_FN_f\left({\displaystyle \frac{1}{3}}+{\displaystyle \frac{8}{3}}\zeta _2\right)`$
for our normalization (2.3) of the expansion parameter.

## 4 Convolutions for the N3LO evolution kernels

Besides the $`(l+1)`$-loop splitting functions and the $`l`$-loop coefficient functions, the NlLO evolution kernels (2.7) involve simple and multiple convolutions of the coefficient functions of lower order. Required at N3LO are the simple and double convolutions of the one-loop coefficient functions $`c_{a,1}`$ with themselves, $`c_{a,1}^2`$ and $`c_{a,1}^3`$, and the convolutions of the one- with the two-loop coefficient functions, $`c_{a,1}c_{a,2}`$, see Eq. (LABEL:eq26). Especially the latter lead to rather complex exact expressions. These terms, however, do not require any attention if the evolution is carried out using the moment-space technique , as in $`N`$-space the convolutions reduce to products. On the other hand, many analyses of data on structure functions are performed using ‘brute-force’ $`x`$-space programs for solving the evolution equations. For application in such programs we provide compact and accurate parametrizations of the convolution contributions to the evolution kernels up to N3LO.

These approximations are derived analogously to those of the two-loop coefficient functions in ref. : The +-distribution parts are treated exactly (truncating irrational coefficients), see the appendix. The integrable $`x<1`$ terms are fitted to the exact results for $`x10^6`$. Finally the coefficients of $`\delta (1x)`$ are slightly adjusted from their exact values using the lowest integer moments. The resulting parametrizations deviate from the exact results by no more than a few permille. This accuracy applies directly to Eqs. (4.1)–(4) as well as to their convolutions with typical hadronic input distributions.

Using the abbreviations (3.2) the simple convolutions of the one-loop coefficient functions for $`F_{2,\mathrm{NS}}`$ and $`F_3`$ can be written as

$`c_{2,1}^2(x)`$ $`=`$ $`256/9𝒟_364𝒟_2283.157𝒟_1+304.751𝒟_0+346.213\delta (1x)`$ (4.1)
$`26.51L_1^3+192.9L_1^2+198.2L_1+113.0L_1^2L_0`$
$`1.230L_0^3+9.466L_0^2+32.45L_0483.3x410.5`$
and

$`c_{3,1}^2(x)`$ $`=`$ $`256/9𝒟_364𝒟_2283.157𝒟_1+304.751𝒟_0+345.993\delta (1x)`$ (4.2)
$`27.09L_1^3+162.1L_1^2+248.0L_1+91.79L_1^2L_0`$
$`1.198L_0^3+3.054L_0^2+65.54L_0305.3x335.7.`$
The corresponding parametrizations for the double convolutions read

$`c_{2,1}^3(x)`$ $`=`$ $`1024/9𝒟_51280/3𝒟_42757.883𝒟_3+9900.585𝒟_2`$
$`+3917.516𝒟_112573.13𝒟_02851.0\delta (1x)`$
$`151.4L_1^5+118.9L_1^46155L_1^347990L_1^230080L_1+6423L_1^2L_0`$
$`0.35L_0^54.30L_0^4106.7L_0^31257L_0^24345L_0+3618x+8547`$
and

$`c_{3,1}^3(x)`$ $`=`$ $`1024/9𝒟_51280/3𝒟_42757.883𝒟_3+9900.585𝒟_2`$
$`+3917.516𝒟_112573.13𝒟_02888.1\delta (1x)`$
$`138.4L_1^5+409.0L_1^41479L_1^324700L_1^2+9646L_110080L_1^2L_0`$
$`0.119L_0^5+3.126L_0^4+84.84L_0^3+288.7L_0^2+264.87L_0+15410x+15890.`$
The convolutions $`c_{a,2}c_{a,1}`$ for $`F_{2,\mathrm{NS}}`$ in electromagnetic DIS and for the charged-current combination $`F_3^{\nu +\overline{\nu }}`$ are parametrized as

$`[c_{2,2}c_{2,1}](x)=`$
$`1536/27𝒟_5343.702𝒟_4633.29𝒟_3+5958.86𝒟_26805.10𝒟_12464.47𝒟_0`$
$`101.7L_1^5155.1L_1^46553L_1^323590L_1^210620L_1+9290L_1^2L_00.35L_0^5`$
$`+0.64L_0^4+92.93L_0^3+761.9L_0^2+2450L_01251x+6286+8609.2\delta (1x)`$
$`+N_f\{7.0912𝒟_455.3087𝒟_3+18.629𝒟_2+619.865𝒟_1584.260𝒟_0`$
$`11.71L_1^4+60.82L_1^3618.0L_1^21979L_1919.6L_1^2L_0+0.48L_0^4`$
$`1.08L_0^343.83L_0^2125.5L_0295.1x+522.4809.14\delta (1x)\}`$
and

$`[c_{3,2}c_{3,1}](x)=`$
$`1536/27𝒟_5343.702𝒟_4633.29𝒟_3+5958.86𝒟_26805.10𝒟_12464.47𝒟_0`$
$`77.39L_1^5+289.1L_1^42823L_1^312500L_1^2+25420L_1+9515L_1^2L_00.524L_0^5`$
$`6.104L_0^4+39.23L_0^3+553.5L_0^2+1393L_0+20080x+2548+8478.2\delta (1x)`$
$`+N_f\{7.0912𝒟_455.3087𝒟_3+18.629𝒟_2+619.865𝒟_1584.260𝒟_0`$
$`14.30L_1^4+10.47L_1^3775.2L_1^22458L_1392.9L_1^2L_0+0.482L_0^4`$
$`+2.541L_0^341.04L_0^2223.9L_0891.9x+468.0803.43\delta (1x)\}.`$
The corresponding results for the charged-current quantities $`F_{2,\mathrm{NS}}`$ and $`F_3^{\nu \overline{\nu }}`$ are obtained by replacing the second and third line of Eqs. (4) and (4), respectively, by

$`109.2L_1^5243.4L_1^46890L_1^324000L_1^2+10840L_1+9144L_1^2L_00.45L_0^5`$
$`+1.80L_0^4+114.0L_0^3+856.6L_0^2+2602L_0711.6x+6298+8569.2\delta (1x)`$
and

$`77.39L_1^5+295.5L_1^42587L_1^310580L_1^2+30580L_1+6461L_1^2L_00.404L_0^5`$
$`5.525L_0^4+23.80L_0^3+484.7L_0^2+1577L_0+22220x+3349+8485.0\delta (1x).`$
## 5 Numerical results for the scaling violations

In this section we illustrate the effect of the next-to-next-to-next-to-leading order (N3LO) contributions to the physical evolution kernels (2.7)–(2) for the electromagnetic structure function $`F_{2,\mathrm{NS}}`$ and the charged-current combination $`F_3^{\nu +\overline{\nu }}`$ henceforth simply denoted $`F_3`$. Specifically, we will discuss the logarithmic derivatives $`\dot{F}_ad\mathrm{ln}F_a/d\mathrm{ln}Q^2`$ calculated at a fixed reference scale $`Q^2=Q_0^2`$ for the initial conditions

$$F_{2,\mathrm{NS}}(x,Q_0^2)=xF_3(x,Q_0^2)=x^{0.5}(1x)^3.$$
(5.1)
The simple model shape (5.1) incorporates the most important features of non-singlet $`x`$-distributions of nucleons. Its overall normalization is irrelevant for the logarithmic derivatives considered here. The reference scale $`Q_0^2`$ is specified via

$$\alpha _s(\mu _r^2=Q_0^2)=\mathrm{\hspace{0.25em}0.2}$$
(5.2)
irrespective of the order of the expansion. Eq. (5.2) corresponds to $`Q_0^230`$ GeV2, a scale typical for fixed-target DIS, for $`\alpha _s(M_Z^2)0.116`$ beyond leading order (LO). The same input (5.1) and (5.2) is chosen at all orders and for both structure functions in order to facilitate a direct comparison of the various contributions to the evolution kernels (2.7). All graphical illustrations below refer to $`N_f=4`$ effectively massless quark flavours.

Before turning to the numerical results we have to specify our treatment of the four-loop splitting functions $`P_{a,3}`$. These functions enter Eq. (2.7) at order $`\alpha _s^4`$ (N3LO) together with the three-loop coefficient functions (3) – (3) and the convolutions (4) – (4). As already mentioned in the introduction, the size of the two- and three-loop terms in the expansion of the non-singlet splitting functions strongly indicates that the effects of $`P_{a,3}`$ are very small in the $`x`$-region addressed by the present study, $`x>10^2`$. Hence a rather rough estimate of these quantities is sufficient here. We have checked that the \[0/1\] Pade approximation444 A brief discussion of the Padé summations and the resulting higher-order approximations can be found in the next section. gives a reasonable, though not particularly accurate estimate of the three-loop non-singlet splitting functions $`P_{a,2}^N`$ in $`N`$-space for $`N_f=3\mathrm{}5`$. Thus we choose the Mellin inverse of

$$P_{a,3}^N\eta [P_{a,3}^N]_{[1/1]\mathrm{Pad}\stackrel{´}{\mathrm{e}}},\eta =\mathrm{\hspace{0.25em}0}\mathrm{}\mathrm{\hspace{0.17em}2}$$
(5.3)
as our estimate of $`P_{a,3}(x)`$, i.e., we assign a 100% error to the predictions of the \[1/1\] Padé summation (the \[0/2\] Padé results are similar). The results obtained by combining Eq. (5.3) for $`\eta =2`$ with $`c_{a,3}^A`$ in Eqs. (3) – (3) are denoted by N3LOA in the figures below, those using $`\eta =0`$ and $`c_{a,3}^B`$ by N3LOB. As in section 3 the central predictions $`1/2(\mathrm{N}^3\mathrm{LO}_A+\mathrm{N}^3\mathrm{LO}_B)`$ are not shown separately.

The logarithmic scale derivatives $`\dot{F}_{2,\mathrm{NS}}`$ and $`\dot{F}_3`$ resulting from Eqs. (5.1), (5.2) are shown in the left parts of Fig. 4 and Fig. 5, respectively, for the standard choice $`\mu _r^2=Q^2`$ of the renormalization scale. At $`x\stackrel{<}{}\mathrm{\hspace{0.25em}0.5}`$ the size of the NNLO and N3LO corrections can barely be read on the scale of these graphs, therefore the differences $`\dot{F}_a(\dot{F}_a)_{\mathrm{NNLO}}`$ are displayed on a larger scale in the left parts of both figures. The difference of the N3LOA and N3LOB results is very small down to $`x10^2`$ even on this enlarged scale, demonstrating that our approximations (3) – (3) and (5.3) are completely sufficient in this region of $`x`$. For both structure functions the N3LO corrections are large only at very large $`x`$, where the kernels are dominated by the universal soft-gluon contributions. Towards smaller $`x`$ the N3LO effects rapidly decrease, e.g., from 6% at $`x=0.85`$ to 2% at $`x=0.65`$. The corresponding NNLO contributions amount to about 12% and 6%, respectively. At $`10^2\stackrel{<}{}x\stackrel{<}{}\mathrm{\hspace{0.17em}0.6}`$ the N3LO corrections are particularly small for $`F_{2,\mathrm{NS}}`$.

The dependence of these scale derivatives on the renormalization scale555As already mentioned below Eq. (2.4), we use the $`\overline{\text{MS}}`$ renormalization scheme throughout this study. is illustrated in Fig. 6 and Fig. 7. In the former figure the consequences of varying $`\mu _r`$ are shown for $`\dot{F}_{2,\mathrm{NS}}`$ at six representative values of $`x`$ (note that the scales of the ordinates are different in all six parts). Here we vary $`\mu _r`$ over a rather wide range, $`\frac{1}{8}Q^2\mu _r^28Q^2`$, corresponding to $`\mathrm{\hspace{0.17em}0.29}\stackrel{>}{}\alpha _s(\mu _r^2)\stackrel{>}{}\mathrm{\hspace{0.17em}0.15}`$ for the initial condition (5.2). In the latter figure we display the absolute scale uncertainties of $`\dot{F}_{2,\mathrm{NS}}`$ and $`\dot{F}_3`$ at $`Q^2=Q_0^2`$, estimated by

$$\mathrm{\Delta }_{\mathrm{abs}}\dot{F}_a\frac{1}{2}\left\{\mathrm{max}\left[\dot{F}_a(\mu _r^2=\frac{1}{4}Q^2\mathrm{}4Q^2)\right]\mathrm{min}\left[\dot{F}_a(\mu _r^2=\frac{1}{4}Q^2\mathrm{}4Q^2)\right]\right\},$$
(5.4)
i.e., using the smaller conventional interval $`\frac{1}{2}Q\mathrm{}2Q`$ for $`\mu _r`$. Also shown here are the further improvements resulting from including the approximate $`\alpha _s^5`$ (N4LO) contributions to Eqs. (2.7) – (2) discussed in the next section.

Our new N3LO results represent a clear improvement over the NNLO stability  for all $`x`$-values of Fig. 6 except for $`x=0.05`$ (here, however, the absolute spread is very small, see Fig. 7), where the difference between the N3LOA and N3LOB results at small $`\mu _r`$ becomes comparable to the $`\mu _r`$-variation at NNLO. This enhanced sensitivity at small scales is due to the larger values of $`\alpha _s`$ up to almost 0.3, which enter the approximate contributions to Eq. (2.7) as $`\alpha _s^4`$. The present approximation uncertainties of the N3LO results are actually dominated by the estimate (5.3) for the four-loop splitting functions, not by the residual uncertainties of the three-loop coefficient functions quantified in Eqs. (3)–(3).

The results shown in Fig. 7 correspond to relative uncertainties $`(\mathrm{\Delta }_{\mathrm{abs}}\dot{F}_a)/\dot{F}_a`$ of 8% at NNLO and 5% and 3% at N3LO and N4LO, respectively, for both $`F_{2,\mathrm{NS}}`$ and $`F_3`$ at $`x=0.85`$. The corresponding figures at $`x=0.65`$ read 5% (NNLO), 2% (N3LO) and 1% (N4LO). These scale uncertainties are rather similar to the relative size of the highest-order contributions at $`\mu _r^2=Q^2`$ in Fig. 4 and Fig. 5 (for the N4LO contribution see Fig. 11 in section 6). Hence the $`\mu _r`$-variation (5.4) and the size of the last contribution included in Eq. (2.7) yield consistent uncertainty estimates in this region of $`x`$. At smaller $`x`$ the absolute scale uncertainties are very small at N3LO and N4LO. For $`\mathrm{\Delta }_{\mathrm{abs}}\dot{F}_{2,\mathrm{NS}}`$ values even below 0.001 are reached for $`x0.5`$ at N3LO and $`x0.6`$ at N4LO.

We conclude this section by illustrating the effect of the higher-order terms in Eq. (2.7) on the determination of $`\alpha _s`$ from the scaling violations of non-singlet structure functions. For this illustration we assume that $`F_{2,\mathrm{NS}}`$ and $`F_3`$ at $`Q_0^230`$ GeV2 are given by Eq. (5.1) with negligible uncertainty. The resulting average N3LO predictions of $`\dot{F}_{2,\mathrm{NS}}`$ and $`\dot{F}_3`$ for $`\mu _r^2=Q_0^2`$ and $`\alpha _s=0.2`$ are employed as model data at $`x_k=0.1k0.05`$ with $`k=1,\mathrm{},8`$. Roughly following the experimental pattern, we assign errors of 0.005 for $`k=2,\mathrm{}\mathrm{\hspace{0.17em}6}`$, of 0.01 for $`k=1,7`$ and of 0.02 for $`k=8`$ to these data points (for Eq. (5.9) only the relative size of these errors is relevant). Again already including the N4LO estimate obtained in the next section, the fits of $`\alpha _s(Q_0^2)`$ to these model data yield

$`\alpha _s(Q_0^2)_{\mathrm{NLO}}^{}`$ $`=`$ $`\mathrm{\hspace{0.25em}0.2080}\begin{array}{c}+\mathrm{\hspace{0.25em}\hspace{0.25em}0.021}\hfill \\ \mathrm{\hspace{0.25em}\hspace{0.25em}0.013}\hfill \end{array},0.2035\begin{array}{c}+\mathrm{\hspace{0.25em}\hspace{0.25em}0.019}\hfill \\ \mathrm{\hspace{0.25em}\hspace{0.25em}0.011}\hfill \end{array}`$ (5.9)
$`\alpha _s(Q_0^2)_{\mathrm{NNLO}}^{}`$ $`=`$ $`\mathrm{\hspace{0.25em}0.2010}\begin{array}{c}+\mathrm{\hspace{0.25em}\hspace{0.25em}0.008}\hfill \\ \mathrm{\hspace{0.25em}\hspace{0.25em}0.0025}\hfill \end{array},0.1995\begin{array}{c}+\mathrm{\hspace{0.25em}\hspace{0.25em}0.0065}\hfill \\ \mathrm{\hspace{0.25em}\hspace{0.25em}0.0015}\hfill \end{array}`$ (5.14)
$`\alpha _s(Q_0^2)_{\mathrm{N}^3\mathrm{LO}}^{}`$ $`=`$ $`\mathrm{\hspace{0.25em}0.2000}\begin{array}{c}+\mathrm{\hspace{0.25em}\hspace{0.25em}0.003}\hfill \\ \mathrm{\hspace{0.25em}\hspace{0.25em}0.001}\hfill \end{array},0.2000\begin{array}{c}+\mathrm{\hspace{0.25em}\hspace{0.25em}0.0025}\hfill \\ \mathrm{\hspace{0.25em}\hspace{0.25em}0.0005}\hfill \end{array}`$ (5.19)
$`\alpha _s(Q_0^2)_{\mathrm{N}^4\mathrm{LO}}^{}`$ $`=`$ $`\mathrm{\hspace{0.25em}0.2000}\begin{array}{c}+\mathrm{\hspace{0.25em}\hspace{0.25em}0.0015}\hfill \\ \mathrm{\hspace{0.25em}\hspace{0.25em}0.0005}\hfill \end{array},0.2005\begin{array}{c}+\mathrm{\hspace{0.25em}\hspace{0.25em}0.0015}\hfill \\ \mathrm{\hspace{0.25em}\hspace{0.25em}0.0005}\hfill \end{array}`$ (5.24)
where the first column refers to $`F_{2,\mathrm{NS}}`$ and the second to $`F_3`$. The central values represent the respective results for $`\mu _r^2=Q_0^2`$, and the errors are due to the renormalization scale variation $`\frac{1}{4}Q^2\mu _r^24Q^2`$, for the N3LO and N4LO cases combined with the approximation uncertainties. Unlike the NNLO terms, the N3LO and N4LO corrections do not cause significant shifts of the central values, but just lead to a reduction of the $`\mu _r`$ uncertainties which reach about $`\pm 1\%`$ at N3LO. The difference of the NLO and NNLO central results for $`F_3`$ is half as large as that for $`F_{2,\mathrm{NS}}`$. This effect is due to larger positive corrections to the logarithmic derivative at $`x<0.4`$ in the former case (see Figs. 4 and 5), which counteract the effect of the negative corrections at large $`x`$ in the fit. As far as Eq. (5.9) can be compared to the fits of real data in refs.  (where higher-twist contributions affecting the central values are included), our finding are consistent with those results.

## 6 Resummations and optimizations

Finally we address the predictions of the soft-gluon resummation, the ECH and PMS scheme optimizations, and the Padé approximations for the physical evolution kernel $`K_a`$ in Eq. (2.7). Being especially interested in the region of large-$`x/`$large-$`N`$, where the higher-order corrections are large but similar for $`F_{2,\mathrm{NS}}`$ and $`F_3`$, we will for brevity focus on the former, more accurately measured quantity. Also in this section the numerical results are given for $`N_f=4`$ and the initial conditions (5.1) and (5.2). We will mainly consider the predictions of the above-mentioned approaches at fixed order in $`\alpha _s`$, and only at the end briefly turn to the all-order results for the soft-gluon exponentiation and the Padé summations.

The NlLO predictions of the soft-gluon resummation for the kernels (2.7) are given by the terms $`a_s^{l+1}𝒦_{\mathrm{res},l}`$ in Eq. (2). Recall that the leading logarithmic (LL), next-to-leading logarithmic (NLL) and next-to-next-to-leading logarithmic (NNLL) contributions behave as $`\mathrm{ln}^{l+1}N`$, $`\mathrm{ln}^lN`$ and (at $`l2`$) $`\mathrm{ln}^{l1}N`$, respectively. The terms in the $`l`$-loop coefficient functions $`c_{a,l}`$ proportional to $`\mathrm{ln}^kN`$ with $`k=l+2,\mathrm{},2l`$ cancel in the combinations (LABEL:eq26) for $`l2`$. This implies that, from $`l=5`$ onwards, actually none of the four leading $`\mathrm{ln}^kN`$ terms of $`c_{a,l}`$ presently fixed by the soft-gluon exponentiation (2.11) contributes to the NlLO kernels (2.7). Consequently, we expect a pattern for the numerical soft-gluon approximations to the physical kernels which is rather different from that discussed in ref.  for the $`\overline{\text{MS}}`$ coefficient functions.

The cumulative effect of soft-gluon terms at NNLO ($`\alpha _s^3`$) and N3LO ($`\alpha _s^4`$) is compared to the (approximate) full results in Fig. 8. The two solid curves in the right plot refer to the N3LOA and N3LOB approximation discussed below Eq. (5.3), the two NNLL results to $`\xi =8`$ and $`\xi =12`$ in Eq. (3.14). For the (undisplayed) NLO contribution the LL and NNLL predictions are considerably smaller and larger, respectively, than the full result, whereas the inclusion of also the $`a_s^2N^0`$ term arising from soft and virtual emissions leads to a reasonable approximation. Combined with this situation, the results of Fig. 8 indicate that the number of soft-gluon logarithms required for a realistic approximation at NlLO systematically increases with the order $`l`$: The full NLO, NNLO, and N3LO curves run between the LL and NLL, close to the NLL, and between the NLL and NNLL results, respectively. The NNLL soft-gluon contribution may thus be expected to represent a reasonable estimate for the N4LO ($`\alpha _s^5)`$ term of Eq. (2.7) at large $`x/`$large $`N`$. As shown in Fig. 9 666Already Fig. 3 demonstrates that the universal soft-gluon terms do not provide a good approximation for $`x<0.7`$, hence the corresponding $`x`$-space results in Figs. 9 and 12 are shown only at larger $`x`$., however, the spread due to the present uncertainty (3.14) of the parameters $`B_2`$ and $`D_2`$ entering Eq. (2) is unfortunately rather large. Moreover, even if with this uncertainty removed, e.g., by a future exact calculation of the three-loop coefficient functions, Figs. 8 and 9 indicate that the soft-gluon resummation can hardly be expected to provide accurate information on the N4LO term, even for moments as large as $`N30`$.

The renormalization-scheme optimizations assume that that the higher-order corrections to the NlLO physical kernels $`𝒦_a^{(l)N}`$ in $`N`$-space given by

$$\frac{d\mathrm{ln}_a^N}{d\mathrm{ln}Q^2}=𝒦_a^{(l)N}=a_sK_0^N(1+a_sr_{a,1}^N+\mathrm{}+a_s^lr_{a,l}^N)$$
(6.1)
are small in a certain ‘optimal’ scheme. The principle of minimal sensitivity (PMS) proposed in ref.  selects this scheme by the requirement

$$\frac{d^{}𝒦_a^N(\mathrm{RS})}{d(\mathrm{RS})}=\mathrm{\hspace{0.25em}0},$$
(6.2)
where $`d/d(\mathrm{RS})`$ abbreviates the derivatives with respect to the $`l`$ independent parameters specifying the renormalization scheme at NlLO. In the effective charge (ECH) method of ref. , on the other hand, these parameters are chosen such that

$$r_{a,1}^N=\mathrm{}=r_{a,l}^N=\mathrm{\hspace{0.25em}0}.$$
(6.3)
Assuming that in these schemes the next terms $`r_{a,l+1}^N`$ are not just small but vanishing, the transformation back to $`\overline{\text{MS}}`$ (or any other scheme) leads to the respective PMS and ECH predictions for this quantity in terms of $`r_{a,1}^N\mathrm{}r_{a,l}^N`$ and the coefficients (2.4) of the $`\beta `$-function. Up to $`r_4`$ these predictions are explictly given in Eqs. (6) – (11) and (13) – (17) of ref. , thus we refrain from repeating them here.

Another approach for estimating the higher-order corrections is provided by the Padé summation of the perturbation series, for QCD in detail discussed, e.g., in refs. . In this method $`𝒦_a^{(l)N}`$ in Eq. (6.1) is replaced by

$$𝒦_{a,[𝒩/𝒟]}^N=a_sK_0\frac{1+a_sp_{a,1}^N+\mathrm{}+a_s^𝒩p_{a,𝒩}^N}{1+a_sq_{a,1}^N+\mathrm{}+a_s^𝒟q_{a,𝒟}^N}$$
(6.4)
with

$$𝒟\mathrm{\hspace{0.25em}1}\text{and}𝒩+𝒟=l.$$
(6.5)
The determination of the parameters $`p_i`$ and $`q_j`$ from the $`r_1\mathrm{}r_l`$ of Eq. (6.1) are automatized in programs for symbolic manipulation such as Maple . Expanding $`𝒦_{a,[𝒩/𝒟]}^N`$ to order $`l+1`$ then yields the $`[𝒩/𝒟]`$ Padé predictions for the N$`{}_{}{}^{l+1}LO`$ coefficients $`r_{a,l+1}^N`$. Also these predictions need not to be written down here. Beyond the second-order results there is no obvious relation between the predictions of the scheme optimizations and those of the Padé approximations. Consistent result of these methods for $`r_{l>2}`$ are thus usually considered as evidence of the approximate correctness of these predictions .

The PMS, ECH and Padé results for the NNLO and N3LO $`N`$-space kernels (6.1) are compared in Fig. 10 to the (approximate) full results already shown, on the same scale, in Fig. 8. Disregarding large relative, but small absolute deviations at NNLO for small $`N`$, the PMS and ECH results (which are very similar at NNLO and identical at N3LO) represent good approximations at both orders. The Padé approximations are somewhat smaller, however, this offset seems to decrease with the order in $`\alpha _s`$. In the left part of Fig. 11 we present the corresponding N4LO predictions. The inner three curves have been derived from the central N3LO results of section 5. The PMS and ECH are again very similar, they are not shown separately. The impact of the present uncertainty of the N3LO kernels, dominated by the estimate (5.3) of the four-loop splitting functions, is included in the two dotted curves which represent our final estimate for the N4LO term and its uncertainty. In the right part of Fig. 10 the N4LO corrections to the results of Fig. 4 are compared to the N3LO contribution. Within the large uncertainties of the latter, these results are consistent with the NNLL soft-gluon prediction shown in Fig. 9. The consequences of including the N4LO estimates have been presented in Fig. 7 and Eq. (5.9), respectively, for the renormalization-scale stability and the determination of $`\alpha _s`$.

Finally the $`\alpha _s^{l>5}`$ infinite-order predictions of the soft-gluon resummation (2) (using the minimal prescription contour ) and the Padé approximations (6.4) are compared in Fig. 12. For the present uncertainties (3.14) — the curves in the figure refer to $`\xi =8`$ (upper), $`D_2=0`$ in Eq. (2.18) (middle) and $`\xi =12`$ (lower) — it is not possible to draw any conclusions from the soft-gluon result. The Padé summation, on the other hand, provides rather definite predictions: The terms beyond $`\alpha _s^5`$ can be expected to have a very small impact at $`x\stackrel{<}{}\mathrm{\hspace{0.17em}0.75}`$. For our standard reference value $`\alpha _s=0.2`$ their effect reaches about the size of the N4LO and N3LO contributions at $`x0.9`$ and $`x0.95`$, respectively.

## 7 Summary

We have investigated the predictions of massless perturbative QCD for the scaling violations of the most important non-singlet structure functions in unpolarized DIS, extending our previous NNLO results  to N3LO and N4LO for the region $`x>10^2`$. The main objective of this extension is to reduce the theoretical uncertainty of determinations of $`\alpha _s`$ from inclusive DIS to about $`1\%`$, an accuracy which is sufficient to make full use of present and future structure function measurements. Our results also facilitate improved determinations of power-suppressed contributions to the structure-function evolution by fits to data, especially at large $`x`$ where the uncertainties are still sizeable at NNLO.

The new ingredients entering the N3LO physical evolution kernels are the four-loop splitting functions and the three-loop coefficient functions. The impact of the former quantities is expected to be very small in the $`\overline{\text{MS}}`$ scheme at $`x>10^2`$; it has been estimated by a Padé approximation assigned a 100% uncertainty. For the latter quantities we have derived approximate expressions based on the available integer moment  and soft gluon  results. The effect of these functions is very well under control at $`x>10^2`$, and almost perfectly at $`x\stackrel{>}{}\mathrm{\hspace{0.17em}0.1}`$. In fact, the uncertainty of the splitting functions dominates the small residual uncertainties of the evolution kernels. Hence the accuracy of our present N3LO results will be superseded only by a future four-loop calculation.

We have also studied the predictions of the NNLL soft-gluon resummation  and of the Padé, PMS and ECH approximations . Presently the predictions of the resummation for the physical evolution kernels beyond N3LO (in any case applicable only at $`x\stackrel{>}{}0.8`$) suffer from the incomplete determination of the soft-gluon parameters $`B_2`$ and $`D_2`$, a problem which will be removed by forthcoming exact calculation of the three-loop coefficient functions . The Padé, PMS and ECH approximations are found to agree rather well with the NNLO and N3LO results for the evolution kernels; these approaches seem to provide reliable predictions of the effects at N4LO and beyond.

For $`\alpha _s\stackrel{<}{}\mathrm{\hspace{0.17em}0.2}`$ the N3LO and N4LO corrections at $`\mu _rQ`$ are very small at $`x<0.6`$ and $`x<0.8`$, respectively, especially for the most accurately measured structure function $`F_2`$. Consequently the central values of $`\alpha _s`$ determined from the non-singlet scaling violations hardly change any more once the larger NNLO terms have been included, see also ref. . The scale uncertainty of the resulting $`\alpha _s(M_Z^2)`$ is reduced to the unproblematic level of less than 1% at N3LO and 0.5% at N4LO. In order to ensure an overall theoretical accuracy of about 1% also the heavy quark (especially charm) mass effects need to be controlled with this precision. We will address this point in a forthcoming publication.

Fortran subroutines of our approximations of the three-loop coefficient functions in section 3 and of the parametrizations of the convolutions entering the evolution kernels in section 4 can be found at http://www.lorentz.leidenuniv.nl/$``$avogt.

## Acknowledgment

We are grateful to J. Vermaseren for communicating the fourteenth moment of the three-loop coefficient function for $`F_2`$ to us prior to publication. This work has been supported by the European Community TMR research network ‘Quantum Chromodynamics and the Deep Structure of Elementary Particles’ under contract No. FMRX–CT98–0194.

## Appendix: Convolutions of +-distribution

The convolutions $`𝒟_i𝒟_j`$ of the +-distributions (3.2) for $`i+j4`$ are given by

$`𝒟_0𝒟_0`$ $`=`$ $`2𝒟_1\zeta _2\delta (1x)`$
$`𝒟_1𝒟_0`$ $`=`$ $`{\displaystyle \frac{3}{2}}𝒟_2\zeta _2𝒟_0+\zeta _3\delta (1x)`$
$`𝒟_1𝒟_1`$ $`=`$ $`𝒟_32\zeta _2𝒟_1+2\zeta _3𝒟_0{\displaystyle \frac{1}{10}}\zeta _2^2\delta (1x)`$
$`𝒟_2𝒟_0`$ $`=`$ $`{\displaystyle \frac{4}{3}}𝒟_32\zeta _2𝒟_1+2\zeta _3𝒟_0{\displaystyle \frac{4}{5}}\zeta _2^2\delta (1x)`$
$`𝒟_2𝒟_1`$ $`=`$ $`{\displaystyle \frac{5}{6}}𝒟_43\zeta _2𝒟_2+6\zeta _3𝒟_1\zeta _2^2𝒟_0+\left(4\zeta _52\zeta _2\zeta _3\right)\delta (1x)`$
$`𝒟_2𝒟_2`$ $`=`$ $`{\displaystyle \frac{2}{3}}𝒟_54\zeta _2𝒟_3+12\zeta _3𝒟_24\zeta _2^2𝒟_1+\left(16\zeta _58\zeta _2\zeta _3\right)𝒟_0`$
$`+\left(4\zeta _3^2{\displaystyle \frac{46}{35}}\zeta _2^3\right)\delta (1x)`$
$`𝒟_3𝒟_0`$ $`=`$ $`{\displaystyle \frac{5}{4}}𝒟_43\zeta _2𝒟_2+6\zeta _3𝒟_1{\displaystyle \frac{12}{5}}\zeta _2^2𝒟_0+6\zeta _5\delta (1x)`$
$`𝒟_3𝒟_1`$ $`=`$ $`{\displaystyle \frac{3}{4}}𝒟_54\zeta _2𝒟_3+12\zeta _3𝒟_2{\displaystyle \frac{27}{5}}\zeta _2^2𝒟_1+\left(18\zeta _56\zeta _2\zeta _3\right)𝒟_0`$
$`+\left(3\zeta _3^2{\displaystyle \frac{36}{35}}\zeta _2^3\right)\delta (1x)`$
$`𝒟_4𝒟_0`$ $`=`$ $`{\displaystyle \frac{6}{5}}𝒟_54\zeta _2𝒟_3+12\zeta _3𝒟_2{\displaystyle \frac{48}{5}}\zeta _2^2𝒟_1+24\zeta _5𝒟_0{\displaystyle \frac{192}{35}}\zeta _2^3\delta (1x)`$
up to integrable contributions dealt with numerically in section 4. Here $`\zeta _l`$ stand for the Riemann $`\zeta `$-function, and $`\zeta _4`$ and $`\zeta _6`$ have been expressed in terms of $`\zeta _2`$ and $`\zeta _3`$, respectively. A convenient method to derive (and to extend, if required) eqs. (Appendix: Convolutions of +-distribution) is by using the relation between $`𝒟_i`$ and the harmonic sums discussed, for example, in ref. .