Abstract

We present the results of our QCD analysis for nonsinglet
unpolarized quark distributions and structure function
F2(x,Q2) up to next-to-next-to-next-to leading order(N3LO).
In this regards 4-loop anomalous dimension can be obtained from
the Padé approximations. The analysis is based on the Jacobi
polynomials expansion of the structure function. New
parameterizations are derived for the nonsinglet quark
distributions for the kinematic wide range of x and Q2. Our
calculations for nonsinglet unpolarized quark distribution
functions up to N3LO are in good agreement with available
theoretical models. The higher twist contributions of
Fp,d2(x,Q2) are extracted in the large x region in
N3LO analysis. The values of ΛQCD and
αs(M2z) are determined.

pacs:

13.60.Hb, 12.39.-x, 14.65.Bt

I Introduction

Structure functions in deep-inelastic scattering (DIS) and their
scale evolution are closely related to the origins of quantum
chromodynamics (QCD). DIS processes have played and still play a
very important role for our understanding of QCD and nucleon
structure Altarelli:2009za (). In fact, DIS structure
functions have been the subject of detailed theoretical and
experimental investigations. Today, with high-precision data from
the electron proton collider, HERA, and in view of the outstanding
importance of hard scattering processes at proton(anti)proton
colliders like the TEVATRON and the forthcoming Large Hadron
Collider (LHC) at CERN, a quantitative understanding of
deep-inelastic processes is indispensable.

To predict the rates of the various processes, a set of universal
parton distribution functions (PDF’s) is required. On the other
hand all calculations of high energy processes with initial
hadrons, whether within the standard model or exploring new
physics, require PDF’s as an essential input. The reliability of
these calculations, which underpins both future theoretical and
experimental progress, depends on understanding the uncertainties
of the PDF’s. These distribution functions can be determined by
QCD global fits to all the available DIS and related
hard-scattering data. The QCD fits can be performed at leading
order (LO), next-to-leading order (NLO), next-to-next-to-leading
order (N2LO) in the strong coupling αs.

The assessment of PDF’s, their uncertainties and extrapolation to
the kinematics relevant for future colliders such as the LHC have
been an important challenge to high energy physics in recent
years. Over the last couple of years there has been a considerable
improvement in the precision, and in the kinematic range of the
experimental measurements for many of these processes, as well as
new types of data becoming available. In addition, there have been
valuable theoretical developments, which increase the reliability
of the global analysis. It is therefore timely, particularly in
view of the forthcoming experiments at the LHC at CERN, to perform
new global analysis which incorporate all of these improvements. A
lot of efforts and challenges have been done to obtain PDF’s for
the LHC CERN:2009 () which take into account the higher order
corrections Martin:2009iq (); Thorne:2009me (); Nadolsky:2008zw ().

For quantitatively reliable predictions of DIS and hard hadronic
scattering processes, perturbative QCD corrections at the N2LO
and the next-to-next-to-next-to-leading order (N3LO) need to be
taken into account. Based on our experience obtained in a series
of LO, NLO and N2LO analysis Khorramian:2008yh () of the
nonsinglet parton distribution functions, here we extend our work
to N3LO accuracy in perturbative QCD.

In the present paper we perform a QCD analysis of the flavor
nonsinglet unpolarized deep–inelastic charged e(μ)p and
e(μ)d world data BCDMS (); SLAC (); NMC (); H1 (); ZEUS () at N3LO and
derived parameterizations of valence quark distributions
xuv(x,Q2) and xdv(x,Q2) at a starting scale Q20
together with the QCD–scale ΛQCD by using the
Jacobi polynomial expansions. We have therefore used the 3-loop
splitting functions and Padé approximations
Pade1 (); Pade2 (); Pade3 (); Baker (); BenderOrszag () for the evolution of
nonsinglet quark distributions of hadrons.

Previous 3–loop QCD analysis were mainly performed as combined
singlet and nonsinglet analysis MRST03 (); A02 (), partly based
on preliminary, approximative expression of the 3–loop splitting
functions. Other analyses were carried out for fixed moments only
in the singlet and nonsinglet case analyzing neutrino data
SY (); KAT (); SYS (). First results of the nonsinglet analysis were
published in BBG04 (). Very recently a 3–loop nonsinglet
analysis was also carried out in
Refs. Khorramian:2008yh (); Gluck:2006yz (); Blumlein:2006be (). The
results of 4–loop QCD analysis are also reported in
Blumlein:2006be (); ABKM:2009aa (). The results of the present
work are based on the Jacobi polynomials expansion of the
nonsinglet structure function.

The plan of the paper is to recall the theoretical formalism of
the QCD analysis for calculating nonsinglet sector of proton
structure function F2 in Mellin-N space in Sec. II.
Section III explains the Padé approximations and 4-loop
anomalous dimensions. A description of the Jacobi polynomials
and procedure of the QCD fit of F2 data are illustrated in Sec. IV.
The numerical results are illustrated in Sec. V before we
summarize our findings in Sec. VI.

here qi and g represent the quarks and gluons distributions
respectively, qNS stands for the usual flavor nonsinglet
combination and qS stand for the flavor-singlet quark
distribution, qS=∑nfi=1(qi+¯qi). Also, nf
denotes the number of effectively massless flavors. <e2>
represents the average squared charge, and ⊗ denotes the
Mellin convolution which turns into a simple multiplication in
N-space.

The perturbative expansion of the coefficient functions can be
written as

C2,i(x,αs(Q2))=∑n=0(αs(Q2)4π)nC(n)2,i(x).

(2)

In LO, C(0)2,NS(x)=δ(x), C(0)2,PS(x)=C(0)2,g(x)=C(1)2,PS(x)=0
and the singlet-quark coefficient function is decomposed into the
nonsinglet and pure singlet contribution, C(n)2,q≡C(n)2,S=C(n)2,NS+C(n)2,PS. The
coefficient functions C(n)2,i up to N3LO have been
given in Vermaseren:2005qc ().

The nonsinglet structure function F2,NS(x,Q2) up to N3LO and for three
active (light) flavors has the representation

In Eq. (3) q+3=u+¯u−(d+¯d)=uv−dv and
q+8=u+¯u+d+¯d−2(s+¯s)=uv+dv+2¯u+2¯d−4¯s, where s=¯s. Also in
Eq. (4),
Σ(x,Q2)≡Σq=u,d,s(q+¯q)=uv+dv+2¯u+2¯d+2¯s.
Notice that in the above equations
as=as(Q2)≡αs(Q2)/4π denotes the strong coupling
constant and Ci,j are the Wilson coefficients
Vermaseren:2005qc ().

The combinations of parton densities in the nonsinglet regime and
the valence region x≥0.3 for Fp2 in LO is

1xFp2(x,Q2)=[118q+NS,8+16q+NS,3](x,Q2)+29Σ(x,Q2),

(7)

where q+NS,3=uv−dv, q+NS,8=uv+dv and
Σ=uv+dv, since sea quarks can be neglected in the region
x≥0.3. So in the x-space we have

Fp2(x,Q2)=(518xq+NS,8+16xq+NS,3)(x,Q2)

=

49xuv(x,Q2)+19xdv(x,Q2).

(8)

In the above region the combinations of parton densities for
Fd2 are also given by

Fd2(x,Q2)=(518xq+NS,8)(x,Q2)

=

518x(uv+dv)(x,Q2),

(9)

where q+NS,3=uv−dv and Fd2=(Fp2+Fn2)/2 if we
ignore the nuclear effects here. It is important to stress that
the shadowing effect as a nuclear effect may affect our analysis.
The shadowing effect Barone:1992ej (); Kwiecinski:1987tb ()
arising from the gluon recombination and in the small-x region,
the competitive mechanism of nuclear shadowing takes place. It
also depends on the size of the nucleons. According to this effect
we have Fd2=(Fp2+Fn2)/2+δFd2. To obtain the
δFd2 we need to know the generalized vector meson
dominance (VMD) and parton mechanism at low and large values of
Q2 respectively. We found that the value of δFd2 is
important but in low values of x. For example this correction
value at Q2=10 GeV2 and for x>0.1 is too small (∼10−4).
So in the valence region of this analysis, this effect
is negligible in large x and we can use the
Fd2=(Fp2+Fn2)/2 approximately.

In the region x≤0.3 for the difference of the proton and deuteron data we
use

FNS2(x,Q2)

≡

2(Fp2−Fd2)(x,Q2)

(10)

=

13xq+NS,3(x,Q2)=13x(uv−dv)(x,Q2)+23x(¯u−¯d)(x,Q2),

where now q+NS,3=uv−dv+2(¯u−¯d) since
sea quarks cannot be neglected for x smaller than about 0.3.

at Q20=4 GeV2 which gives a good description of the
Drell-Yan dimuon production data E866 (). In this analysis,
like other analyses
Khorramian:2008yh (); Khorramian:2007zz (); Gluck:2006yz (); Blumlein:2006be (); ref19 (); Blumlein:2004pr (),
we used the above distribution for considering the symmetry
breaking of sea quarks. Although, in fact, this parametrization
plays a marginal role in our analysis, in order to find the impact
effect of this distribution, which is essentially used in the
paper, it is desirable to study the QCD fits by varying this
distribution with another asymmetry sea quark distribution which
is derived in other analyses. In Sec. VI we will discuss our
outputs when we change the above sea distribution.

Now these results in the physical region 0<x≤1 can transform
to Mellin-N space by using the Mellin transform to obtain the
moments of the structure function as 1xFk2,

Fk2(N,Q2)≡M[Fk2,N]=∫10dxxN−11xFk2(x,Q2),

(12)

here k denotes the three above cases, i.e. k=p,d,NS. One of
the advantages of Mellin-space calculations is the fact that the
Mellin transform of a convolution of functions in
Eqs. (3,4,5) reduces to a
simple product

M[A⊗B,N]=M[A,N]M[B,N]=A(N)B(N)

(13)

By using the solution of the nonsinglet evolution equation for the
parton densities to 4− loop order, the nonsinglet structure
functions are given by Blumlein:2006be ()

Fk2(N,Q2)

=

(1+asC(1)2,NS(N)+a2sC(2)2,NS(N)+a3sC(3)2,NS(N))Fk2(N,Q20)

(14)

×(asa0)−^P0(N)/β0{1−1β0(as−a0)[^P+1(N)−β1β0^P0(N)]

−12β0(a2s−a20)[^P+2(N)−β1β0^P+1(N)+(β21β20−β2β0)^P0(N)]

+12β20(as−a0)2(^P+1(N)−β1β0^P0(N))2

−13β0(a3s−a30)[^P+3(N)−β1β0^P+2(N)+(β21β20−β2β0)^P+1(N)

+(β31β30−2β1β2β20+β3β0)^P0(N)]

+12β20(as−a0)(a20−a2s)(^P+1(N)−β1β0^P0(N))

×[^P2(N)−β1β0^P1(N)−(β21β20−β2β0)^P0(N)]

−16β30(as−a0)3(^P+1(N)−β1β0^P0(N))3}.

Here as(=αs/4π) and a0 denotes the strong coupling
constant in the scale of Q2 and Q20 respectively. k=p,d
and NS also denotes the three above cases, i.e. proton, deuteron
and nonsinglet structure function. C(m)2,NS(N) are the
nonsinglet Wilson coefficients in O(ams) which can be
found in FP (); NS2 (); Vermaseren:2005qc () and ^Pm denote
also the Mellin transforms of the (m+1)− loop splitting
functions.

Iii Padé approximations and 4-loop anomalous dimensions

In spite of the unknown 4-loop anomalous dimensions, one can
obtain the nonsinglet parton distributions and ΛQCD by
estimating uncalculated fourth-order corrections to the nonsinglet
anomalous dimension. On the other hand the 3–loop Wilson
coefficients are known Vermaseren:2005qc () and now it is
possible to know, which effect has the 4-loop anomalous dimension
if compared to the Wilson coefficient. In this case the 4-loop
anomalous dimension may be obtained from Padé approximations.

Padé approximations have proved to be useful in many physical
applications. Padé approximations may be used either to predict
the next term in some perturbative series, called a Padé
approximation prediction, or to estimate the sum of the entire
series, called Padé summation.

For this purpose we use the Padé approximations of the
perturbative series, discussed in detail for QCD, e.g., in
Refs. Pade1 (); Pade2 (); Pade3 (). Padé approximations
Baker (); BenderOrszag () are rational functions chosen to equal
the perturbative series to the order calculated:

[N/M]=a0+a1x+...+aNxN1+b1x+...+bMxM,

(15)

to the series

S=S0+S1x+...+SN+MxN+M,

(16)

where we set

[N/M]=S+O(xN+M+1),

and write an equation for the coefficients of each power of x.
To continue, let’s go to Mellin-N space.

A generic QCD anomalous dimension expansion in term of as then
may be written in the form

γ(N)=∞∑l=0al+1sγ(l)(N).

(17)

In Mellin-N space and by using this approach we can replace
γ(N) by a rational function in asVermaseren:2005qc (),

˜γ[N/M](N)≡[N/M](N)=p0+asp1(N)+…+aNspN(N)1+asq1(N)+…+aMsqM(N).

(18)

Here M≥1 and N+M=n,
where n stands for the maximal order in as at which the
expansion coefficients γ(n)(N) have been determined from
an exact calculation. The functions pi(N) and qj(N) are
determined from these known coefficients by expanding
Eq. (18) in powers of as. This expansion then also
provides the [N/M] Padé approximate for the
(n+1)-th order quantities γ(n+1).

In this way it is easy to obtain the following results for M=N=1 and for M=0,N=2

˜γ[1/1](N)

≡

[1/1](N)=γ(2)2(N)γ(1)(N),

˜γ[0/2](N)

≡

[0/2](N)=2γ(1)(N)γ(2)(N)γ(0)(N)−γ(1)3(N)γ(0)2(N).

(19)

The strong coupling constant as plays a more central role in
the present paper to the evolution of parton densities. At
NmLO the scale dependence of as is given by

here nf stands for the number of effectively massless quark
flavors and βk denote the coefficients of the usual
four-dimensional ¯¯¯¯¯¯¯¯¯MS beta function of QCD. In complete
4-loop approximation and using the Λ-parametrization, the
running coupling is given by Vogt:2004ns (); Chetyrkin:1997sg ():

as(Q2)

=

1β0LΛ−1(β0LΛ)2b1lnLΛ

(22)

+

1(β0LΛ)3[b21(ln2LΛ−lnLΛ−1)+b2]

+

1(β0LΛ)4[b31(−ln3LΛ+52ln2LΛ+2lnLΛ−12)−3b1b2lnLΛ+b32],

where LΛ≡ln(Q2/Λ2), bk≡βk/β0, and Λ is the QCD scale parameter. The first
line of Eq. (22) includes the 1- and the 2-loop
coefficients, the second line is the 3-loop and the third line
denotes the 4-loop correction. Equation (22) solves the
evolution equation (20) only up to higher orders in
1/LΛ. The functional form of αs(Q2), in 4-loop
approximation and for 6 different values of Λ, is
displayed in Fig. 1. The slope and dependence on the
actual value of Λ is especially pronounced at small Q2,
while at large Q2 both the energy dependence and the dependence
on Λ becomes increasingly feeble. To be able to compare
with other measurements of Λ we adopt the matching of
flavor thresholds at Q2=m2c and Q2=m2b with mc=1.5
GeV and mb=4.5 GeV as described in BAR (); R1998 ().

Iv Jacobi polynomials and the procedure of QCD fits

One of the simplest and fastest possibilities in the structure
function reconstruction from the QCD predictions for its Mellin
moments is Jacobi polynomials expansion. The Jacobi polynomials
are especially suitable for this purpose since they allow one to
factor out an essential part of the x-dependence of structure
function into the weight function parisi ().

According to this method, one can relate the F2 structure
function with its Mellin moments

Fk,Nmax2(x,Q2)

=

xβ(1−x)αNmax∑n=0Θα,βn(x)n∑j=0c(n)j(α,β)Fk2(j+2,Q2),

(23)

where Nmax
is the number of polynomials, k denotes the three cases, i.e.
k=p,d,NS. Jacobi polynomials of order nParisi:1978jv (), Θα,βn(x), satisfy the
orthogonality condition with the weight function wαβ=xβ(1−x)α

∫10dxwαβΘα,βk(x)Θα,βl(x)=δk,l.

(24)

In the above, c(n)j(α,β) are the coefficients
expressed through Γ-functions and satisfying the
orthogonality relation in Eq. (24) and F2(j+2,Q2) are
the moments determined in the previous section. Nmax,
α and β have to be chosen so as to achieve the
fastest convergence of the series on the right-hand side of
Eq. (23) and to reconstruct F2 with the required
accuracy. In our analysis we use Nmax=9, α=3.0 and
β=0.5. The same method has been applied to calculate the
nonsinglet structure function xF3 from their moments
Kataev:1997nc (); Kataev:1998ce (); Kataev:1999bp (); Kataev:2001kk ()
and for polarized structure function xg1Atashbar Tehrani:2007be (); Leader:1997kw (); Khorramian:2007gu (). Obviously the
Q2-dependence of the polarized structure function is defined by
the Q2-dependence of the moments.

The evolution equations allow one to calculate the
Q2-dependence of the parton distributions provided at a
certain reference point Q20. These distributions are usually
parameterized on the basis of plausible theoretical assumptions
concerning their behavior near the end points x=0,1.

In the present analysis we choose the following parametrization
for the valence quark densities in the input scale of Q20=4
GeV2

xqv(x,Q20)=Nqxaq(1−x)bq(1+cq√x+dqx),

(25)

where q=u,d and the normalization factors Nu and
Nd are fixed by ∫10uvdx=2 and ∫10dvdx=1, respectively. By QCD fits of the world data for
Fp,d2, we can extract valence quark densities using the
Jacobi polynomials method. For the nonsinglet QCD analysis
presented in this paper we use the structure function data
measured in charged lepton-proton and deuteron deep-inelastic
scattering. The experiments contributing to the statistics are
BCDMS BCDMS (), SLAC SLAC (), NMC NMC (),
H1 H1 (), and ZEUS ZEUS (). In our QCD analysis we use
three data samples : Fp2(x,Q2), Fd2(x,Q2) in the
nonsinglet regime and the valence quark region x≥0.3 and
FNS2=2(Fp2−Fd2) in the region x<0.3.

The valence quark region may be parameterized by the nonsinglet
combinations of parton distributions, which are expressed through
the parton distributions of valence quarks. Only data with Q2>4GeV2 were included in the analysis and a cut in the hadronic
mass of W2≡(1x−1)Q2+m2N>12.5GeV2 was applied in order to widely eliminate higher
twist (HT) effects from the data samples. After these cuts we are
left with 762 data points, 322 for Fp2, 232 for Fd2, and
208 for FNS2. By considering the additional cuts on the
BCDMS (y>0.35) and on the NMC data(Q2>8 GeV2) the total
number of data points available for the analysis reduce from 762
to 551, because we have 227 data points for Fp2, 159 for
Fd2, and 165 for FNS2.

For data used in the global analysis, most experiments combine
various systematic errors into one effective error for each data
point, along with the statistical error. In addition, the fully
correlated normalization error of the experiment is usually
specified separately. For this reason, it is natural to adopt the
following definition for the effective χ2Stump:2001gu (); Khorramian:2008yh ()

χ2global

=

∑nwnχ2n,(nlabels the different %
experiments)

χ2n

=

(1−NnΔNn)2+∑i(NnFdata2,i−Ftheor2,iNnΔFdata2,i)2.

(26)

For the nth experiment, Fdata2,i, ΔFdata2,i, and Ftheor2,i denote the data value, measurement uncertainty
(statistical and systematic combined) and theoretical value for
the ith data point. ΔNn is the
experimental normalization uncertainty and Nn is an
overall normalization factor for the data of experiment n. The
factor wn is a possible weighting factor (with default value
1). However, we allowed for a relative normalization shift Nn between the different data sets within the normalization
uncertainties ΔNn quoted by the experiments.
For example the normalization uncertainty of the NMC(combined)
data is estimated to be 2.5%.
The normalization shifts Nn were fitted
once and then kept fixed.

Now the sums in χ2global run over all data
sets and in each data set over all data points. The minimization
of the above χ2 value to determine the best parametrization
of the unpolarized parton distributions is done using the program
MINUITMINUIT ().

The one σ error for the parton density xqv as given by
Gaussian error propagation is Blumlein:2006be ()

σ(xqv(x))2=np∑i=1np∑j=1(∂xqv∂pi)(∂xqv∂pj)cov(pi,pj),

(27)

where the sum runs over all fitted parameters. The functions
∂xqv/∂pi are the derivatives of xqv with respect to the fit parameter pi,
and cov(pi,pj) are the elements of the covariance
matrix. The derivatives ∂xqv/∂pi can be
calculated analytically at the input scale Q20. Their values
at Q2 are given by evolution which is performed in Mellin-N
space.

Now we need to discuss the derivatives in Mellin-N space a bit
further. The Mellin-N moment for complex values of N
calculated at the input scale Q20 for the parton density
parameterized as in Eq. (25) is given by

qv(N,aq,bq,cq,dq)

=

NqM(n,aq,bq,cq,dq),

(28)

with the normalization constant

Nq=CqvM(1,aq,bq,cq,dq).

(29)

Here Cqv is the respective number of valence quarks, i.e.
Cuv=2 and Cdv=1. In the above
M(n,aq,bq,cq,dq) is given by

M(n,aq,bq,cq,dq)=B[aq+n−1,bq+1]+cqB[a+n+1/2,b+1]+duB[aq+n,bq+1],

where B[a,b] denotes the Euler beta function for complex
arguments. The general form of the derivative of the Mellin moment
qv with respect to the parameter p is given by

∂qv(N,p)∂p=M(n,p)∂Nq∂p+Nq∂M(n,p)∂p.

(31)

In this analysis only the parameters aq and bq have been
fitted for both the xuv and xdv parametrization while the
other parameters involved are kept fixed after a first
minimization in the MINUIT program, since their errors turned out
to be rather large compared to the central values. Here we want to
show the derivatives uv and dv parton densities with respect
to parameter aq and bq. For example:

f(n,aq)

≡

∂M(n,aq)∂aq=B[aq+n−1,bq+1](ψ[aq+n−1]−ψ[aq+bq+n])+

(32)

cqB[aq+n−1/2,bq+1](ψ[aq+n−1/2]−ψ[a+b+n+1/2])+

dqB[aq+n,bq+1](ψ[aq+n]−ψ[aq+bq+n+1]),

f(n,bq)

≡

∂M(n,bq)∂bq=B[aq+n−1,bq+1](ψ[bq+1]−ψ[aq+bq+n])+

(33)

cqB[aq+n−1/2,bq+1](ψ[1+bq]−ψ[aq+bq+n+1/2])+

dqB[aq+n,bq+1](ψ[bq+1]−ψ[aq+bq+n+1]),

and now we can reach the below derivatives for uv(N) and
dv(N) with respect to parameters aq and bq

∂qv(N,p)∂p=Nq(f(n,p)−f(1,p)M(n,p)/M(1,p)),

(34)

also ψ[n]=dlnΓ(n)/dn is Euler’s ψ-function.

To obtain the error calculation of the structure functions
Fp2, Fd2 , and FNS2 the relevant gradients of the
PDF’s in Mellin space have to be multiplied with the corresponding
Wilson coefficients. This yields the errors as far as the QCD
parameter Λ is fixed and regarded uncorrelated. The error
calculation for a variable Λ is done numerically due to
the nonlinear relation and required iterative treatment in the
calculation of αs(Q2,Λ)Khorramian:2008yh (); Blumlein:2006be ().

V Results

In the QCD analysis of the present paper we used three data sets:
the structure functions Fp2(x,Q2) and Fd2(x,Q2) in
the region of x≥0.3 and the combination of these structure
functions FNS2(x,Q2) in the region of x<0.3 .
Notice that we take into account the cuts Q2>4 GeV2,
W2>12.5 GeV2 for our QCD fits to determine some unknown
parameters. In Fig.(2) the proton, deuteron and
nonsinglet data for Fp2(x,Q2), Fd2(x,Q2) and
FNS2(x,Q2) are shown in the nonsinglet regime and the
valence quark region x≥0.3 indicating the above cuts by a
vertical dashed line. The solid lines correspond to the N3LO
QCD fit. Now, it is possible to take into account the target mass
effects in our calculations. The perturbative form of the moments
is derived under the assumption that the mass of the target hadron
is zero (in the limit Q2→∞). At intermediate
and low Q2 this assumption will begin to break down and the
moments will be subject to potentially significant power
corrections, of order O(m2N/Q2), where mN
is the mass of the nucleon. These are known as target mass
corrections (TMCs) and when included, the moments of flavor
nonsinglet structure function have the form
Georgi:1976ve (); Gluck:2006yz ()

Fk2,TMC(n,Q2)

≡

∫10xn−11xFk2,TMC(x,Q2)dx

(35)

=

Fk2(n,Q2)+n(n−1)n+2(m2NQ2)Fk2(n+2,Q2)

+(n+2)(n+1)n(n−1)2(n+4)(n+3)(m2NQ2)2Fk2(n+4,Q2)+O(m2NQ2)3,

where higher powers than (m2N/Q2)2 are negligible for
the relevant x<0.8 region. By inserting Eq. (35) in
Eq. (23) we have