Nonsinglet parton distribution functions from the precise next-to-next-to-next-to leading order QCD fit

# Nonsinglet parton distribution functions from the precise next-to-next-to-next-to leading order QCD fit

Ali N. Khorramian    H. Khanpour    S. Atashbar Tehrani Physics Department, Semnan University, Semnan, Iran School of Particles and Accelerators, Institute for Research in Fundamental Sciences (IPM), P.O.Box 19395-5531, Tehran, Iran
July 21, 2019
###### Abstract

We present the results of our QCD analysis for nonsinglet unpolarized quark distributions and structure function up to next-to-next-to-next-to leading order(NLO). 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 and . Our calculations for nonsinglet unpolarized quark distribution functions up to NLO are in good agreement with available theoretical models. The higher twist contributions of are extracted in the large region in NLO analysis. The values of and 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 (NLO) in the strong coupling .

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 NLO and the next-to-next-to-next-to-leading order (NLO) need to be taken into account. Based on our experience obtained in a series of LO, NLO and NLO analysis Khorramian:2008yh () of the nonsinglet parton distribution functions, here we extend our work to NLO accuracy in perturbative QCD.

In this work this important problem is studied with the help of the method of the structure function reconstruction over their Mellin moments, which is based on the expansion of the structure function in terms of Jacobi polynomials. This method was developed and applied for different QCD analyses parisi (); Barker (); Krivokhizhin:1987rz (); Krivokhizhin:1990ct (); Chyla:1986eb (); Barker:1980wu (); Kataev:1997nc (); Kataev:1998ce (); Kataev:1999bp (); Kataev:2001kk (); Khorramian:2007zz (); Khorramian:2008zz (); Khorramian:2008zza (); AtashbarTehrani:2009zz (); Khorramian:2009zz1 (); Khorramian:2009zz2 (); Khorramian:2009zz3 (); Khorramian:2009zz4 (). The same method has also been applied in the polarized case in Refs. Leader:1997kw () and Atashbar Tehrani:2007be (); Khorramian:2007gu (); Khorramian:2007zza (); Mirjalili:2007ep (); Mirjalili:2006hf ().

In the present paper we perform a QCD analysis of the flavor nonsinglet unpolarized deep–inelastic charged and world data BCDMS (); SLAC (); NMC (); H1 (); ZEUS () at NLO and derived parameterizations of valence quark distributions and at a starting scale together with the QCD–scale 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 in Mellin- 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 data are illustrated in Sec. IV. The numerical results are illustrated in Sec. V before we summarize our findings in Sec. VI.

## Ii Theoretical formalism of the QCD analysis

In the common factorization scheme the relevant structure function as extracted from the DIS process can be written as ref3 (); ref4 (); ref5 (); Gluck:2006pm ()

 x−1F2(x,Q2) = x−1(F2,NS(x,Q2)+F2,S(x,Q2)+F2,g(x,Q2)) (1) = C2,NS(x,Q2)⊗qNS(x,Q2) + C2,S(x,Q2)⊗qS(x,Q2) + C2,g(x,Q2)⊗g(x,Q2),

here and represent the quarks and gluons distributions respectively, stands for the usual flavor nonsinglet combination and stand for the flavor-singlet quark distribution, . Also, denotes the number of effectively massless flavors. represents the average squared charge, and denotes the Mellin convolution which turns into a simple multiplication in -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, , and the singlet-quark coefficient function is decomposed into the nonsinglet and pure singlet contribution, . The coefficient functions up to NLO have been given in Vermaseren:2005qc ().

The nonsinglet structure function up to NLO and for three active (light) flavors has the representation

 x−1F2,NS(x,Q2) = [C(0)2,q+asC(1)2,NS+a2sC(2)+2,NS+a3sC(3)+2,NS]⊗[118q+8+16q+3](x,Q2) . (3)

The flavor-singlet and gluon contributions in Eq. (1) reads

 x−1F2,S(x,Q2) = 29[C(0)2,q+asC(1)2,q+a2sC(2)2,q+a3sC(3)2,q]⊗Σ(x,Q2) ; (4)
 x−1F2,g(x,Q2) = 29[asC(1)2,g+a2sC(2)2,g+a3sC(3)2,g]⊗g(x,Q2) . (5)

The symbol denotes the Mellin convolution

 [A⊗B](x)=∫10dx1∫10dx2δ(x−x1x2) A(x1)B(x2). (6)

In Eq. (3) and , where . Also in Eq. (4), . Notice that in the above equations denotes the strong coupling constant and are the Wilson coefficients Vermaseren:2005qc ().

The combinations of parton densities in the nonsinglet regime and the valence region for in LO is

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

where , and , since sea quarks can be neglected in the region . So in the -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 are also given by

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

where and 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- region, the competitive mechanism of nuclear shadowing takes place. It also depends on the size of the nucleons. According to this effect we have . To obtain the we need to know the generalized vector meson dominance (VMD) and parton mechanism at low and large values of respectively. We found that the value of is important but in low values of . For example this correction value at =10 GeV and for is too small (). So in the valence region of this analysis, this effect is negligible in large and we can use the approximately.

In the region 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 since sea quarks cannot be neglected for smaller than about 0.3.

The first clear evidence for the flavor asymmetry combination of light parton distributions in nature came from the analysis of NMC at CERN to study of the Gottfried sum rule Amaudruz:1991at (). In our calculation we supposed the distribution Gluck:2006yz (); Blumlein:2006be (); ref19 (); Blumlein:2004pr ()

 x(¯d−¯u)(x,Q20)=1.195x1.24(1−x)9.10(1+14.05x−45.52x2) , (11)

at GeV 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 can transform to Mellin- space by using the Mellin transform to obtain the moments of the structure function as ,

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

here denotes the three above cases, i.e. . 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 and denotes the strong coupling constant in the scale of and respectively. and also denotes the three above cases, i.e. proton, deuteron and nonsinglet structure function. are the nonsinglet Wilson coefficients in which can be found in FP (); NS2 (); Vermaseren:2005qc () and denote also the Mellin transforms of the 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 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 . To continue, let’s go to Mellin- space.

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

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

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

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

Here and , where stands for the maximal order in at which the expansion coefficients have been determined from an exact calculation. The functions and are determined from these known coefficients by expanding Eq. (18) in powers of . This expansion then also provides the Padé approximate for the ()-th order quantities .

In this way it is easy to obtain the following results for and for

 ˜γ[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 plays a more central role in the present paper to the evolution of parton densities. At the scale dependence of is given by

 dasdlnQ2=βNmLO(as)=−m∑k=0ak+2sβk. (20)

The expansion coefficients of the -function of QCD are known up to , i.e., NLO Tarasov:1980au (); Larin:1993tp ()

 β0 = 11−2/3 nf, β1 = 102−38/3 nf, β2 = 2857/2−5033/18 nf+325/54 n2f, β3 = 29243.0−6946.30 nf+405.089 n2f+1093/729 n3f , (21)

here stands for the number of effectively massless quark flavors and denote the coefficients of the usual four-dimensional 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Λ)2 b1lnLΛ (22) + 1(β0LΛ)3[b21(ln2LΛ−lnLΛ−1)+b2] + 1(β0LΛ)4[b31(−ln3LΛ+52ln2LΛ+2lnLΛ−12)−3b1b2lnLΛ+b32] ,

where , , 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 . The functional form of , 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 , while at large 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 and with GeV and 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 -dependence of structure function into the weight function parisi ().

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

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

where is the number of polynomials, denotes the three cases, i.e. . Jacobi polynomials of order Parisi:1978jv (), , satisfy the orthogonality condition with the weight function

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

In the above, are the coefficients expressed through -functions and satisfying the orthogonality relation in Eq. (24) and are the moments determined in the previous section. , 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 with the required accuracy. In our analysis we use , and . The same method has been applied to calculate the nonsinglet structure function from their moments Kataev:1997nc (); Kataev:1998ce (); Kataev:1999bp (); Kataev:2001kk () and for polarized structure function Atashbar Tehrani:2007be (); Leader:1997kw (); Khorramian:2007gu (). Obviously the -dependence of the polarized structure function is defined by the -dependence of the moments.

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

In the present analysis we choose the following parametrization for the valence quark densities in the input scale of GeV

 xqv(x,Q20)=Nq xaq(1−x)bq(1+cq√x+dq x) , (25)

where and the normalization factors and are fixed by and , respectively. By QCD fits of the world data for , 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 : , in the nonsinglet regime and the valence quark region and in the region .

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 were included in the analysis and a cut in the hadronic mass of 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 , 232 for , and 208 for . By considering the additional cuts on the BCDMS () and on the NMC data( GeV) the total number of data points available for the analysis reduce from 762 to 551, because we have 227 data points for , 159 for , and 165 for .

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 Stump: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 experiment, , , and denote the data value, measurement uncertainty (statistical and systematic combined) and theoretical value for the data point. is the experimental normalization uncertainty and is an overall normalization factor for the data of experiment . The factor is a possible weighting factor (with default value 1). However, we allowed for a relative normalization shift between the different data sets within the normalization uncertainties quoted by the experiments. For example the normalization uncertainty of the NMC(combined) data is estimated to be 2.5%. The normalization shifts were fitted once and then kept fixed.

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

The one error for the parton density 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 are the derivatives of with respect to the fit parameter , and are the elements of the covariance matrix. The derivatives can be calculated analytically at the input scale . Their values at are given by evolution which is performed in Mellin- space.

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

 qv(N,aq,bq,cq,dq) = Nq M(n,aq,bq,cq,dq) , (28)

with the normalization constant

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

Here is the respective number of valence quarks, i.e. =2 and =1. In the above 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 denotes the Euler beta function for complex arguments. The general form of the derivative of the Mellin moment with respect to the parameter is given by

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

In this analysis only the parameters and have been fitted for both the and 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 and parton densities with respect to parameter and . 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 and with respect to parameters and

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

also is Euler’s -function.

To obtain the error calculation of the structure functions , , and 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 Khorramian:2008yh (); Blumlein:2006be ().

## V Results

In the QCD analysis of the present paper we used three data sets: the structure functions and in the region of and the combination of these structure functions in the region of  . Notice that we take into account the cuts GeV, GeV for our QCD fits to determine some unknown parameters. In Fig.(2) the proton, deuteron and nonsinglet data for , and are shown in the nonsinglet regime and the valence quark region indicating the above cuts by a vertical dashed line. The solid lines correspond to the NLO 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 ). At intermediate and low this assumption will begin to break down and the moments will be subject to potentially significant power corrections, of order , where 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 are negligible for the relevant region. By inserting Eq. (35) in Eq. (23) we have

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