Nonsinglet parton distribution functions from the precise next-to-next-to-next-to leading order QCD fit
Ali N. Khorramian $^{a,b}$
We present the results of our QCD analysis for nonsinglet unpolarized quark distributions and structure function $F_2(x,Q^2)$ up to next-to-next-to-next-to leading order(N$^3$LO). In this regards 4-loop anomalous dimension can be obtained from the Pad\'e 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 $Q^2$. Our calculations for nonsinglet unpolarized quark distribution functions up to N$^3$LO are in good agreement with available theoretical models. The higher twist contributions of $F_2^{p,d}(x,Q^2)$ are extracted in the large $x$ region in N$^3$LO analysis. The values of $\Lambda_{QCD}$ and $\alpha_s(M_z^2)$ are determined.
PRD2010-vII.tex
%\documentclass[twocolumn,showpacs,preprintnumbers,amsmath,amssymb]{revtex4}
\documentclass[preprint,showpacs,preprintnumbers,amsmath,amssymb]{revtex4}

% Some other (several out of many) possibilities
%\documentclass[preprint,aps]{revtex4}
%\documentclass[preprint,aps,draft]{revtex4}
%\documentclass[prb]{revtex4}% Physical Review B

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math


\newcommand\GeV{\,\mbox{GeV}}



%\nofiles

\begin{document}
\title{Nonsinglet parton distribution functions from the precise next-to-next-to-next-to leading order
QCD fit}
\author{Ali N. Khorramian $^{a,b}$}
\email{Khorramiana@theory.ipm.ac.ir}
\homepage{http://particles.ipm.ir/}
\author{H. Khanpour $^{a}$}
\email{hamzeh_khanpour@nit.ac.ir}
\author{S. Atashbar Tehrani $^{b}$}
\email{Atashbar@ipm.ir}
\affiliation{%
$^{(a)}$ Physics Department, Semnan University, Semnan, Iran}
\affiliation{%
$^{(b)}$ School of Particles and Accelerators, Institute for
Research in Fundamental Sciences (IPM), P.O.Box 19395-5531,
Tehran, Iran}
%This line break forced with \textbackslash\textbackslash
%}%
%\preprint{IPM/07-05-2009-SPA}
\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\begin{abstract}
We present the results of our QCD analysis for nonsinglet
unpolarized quark distributions and structure function
$F_2(x,Q^2)$ up to next-to-next-to-next-to leading order(N$^3$LO).
In this regards 4-loop anomalous dimension can be obtained from
the Pad\'e 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 $Q^2$. Our
calculations for nonsinglet unpolarized quark distribution
functions up to N$^3$LO are in good agreement with available
theoretical models.  The higher twist contributions of
$F_2^{p,d}(x,Q^2)$ are extracted in the large $x$ region in
N$^3$LO analysis. The values of $\Lambda_{QCD}$ and
$\alpha_s(M_z^2)$ are determined.
\end{abstract}

\pacs{13.60.Hb, 12.39.-x, 14.65.Bt}% PACS, the Physics and Astronomy
                             % Classification Scheme.
%\keywords{Suggested keywords}%Use showkeys class option if keyword
                              %display desired
\maketitle

\section{\label{sec:Introduction}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 \cite{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 (N$^2$LO) in the strong coupling $\alpha_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 \cite{CERN:2009} which take into account the higher order
corrections \cite{Martin:2009iq,
%Thorne:2009ky,
Thorne:2009me,Nadolsky:2008zw}.

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

%The earlier studies on this topic [2] have already demonstrated
%that ...

 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
\cite{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. \cite{Leader:1997kw} and \cite{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 $e(\mu) p$ and
$e(\mu) d$ world data \cite{BCDMS,SLAC,NMC,H1,ZEUS} at N$^3$LO and
derived parameterizations of valence quark distributions
$xu_v(x,Q^2)$ and $xd_v(x,Q^2)$ at a starting scale $Q_0^2$
together with the QCD--scale $\Lambda_{\rm QCD}$ by using the
Jacobi polynomial expansions.  We have therefore used the 3-loop
splitting functions and Pad\'e approximations
\cite{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 \cite{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
\cite{SY,KAT,SYS}. First results of the nonsinglet analysis were
published in \cite{BBG04}. Very recently a 3--loop nonsinglet
analysis was also carried out in
Refs.~\cite{Khorramian:2008yh,Gluck:2006yz,Blumlein:2006be}. The
results of 4--loop QCD analysis are also reported in
\cite{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 $F_2$ in Mellin-$N$ space in Sec.~II.
Section~III explains the Pad\'e approximations and 4-loop
anomalous dimensions. A description of the Jacobi polynomials
 and procedure of the QCD fit of $F_2$ data are illustrated in Sec.~IV.
The numerical results are illustrated in Sec.~V before we
summarize our findings in Sec.~VI.
\section{\label{sec:QCD Fit}Theoretical formalism of the QCD analysis}
In the common $\overline{\rm MS}$ factorization scheme the
relevant $F_2$ structure function as extracted from the DIS $ep$
process can be written as \cite{ref3,ref4,ref5,Gluck:2006pm}
\begin{eqnarray}
\label{Eq:F2-tot} x^{-1}F_2(x,Q^2) &=& x^{-1}\left(F_{2,{\rm
NS}}(x,Q^2)+ F_{2,S}(x,Q^2)+ F_{2,g}(x,Q^2)\right) \nonumber \\&=&
C_{2,{\rm NS}}(x,Q^2)\otimes q_{\rm NS}(x,Q^2)\nonumber \\&+&<e^2>
C_{2,{\rm S}}(x,Q^2)\otimes q_{\rm S}(x,Q^2)\nonumber
\\&+&<e^2>C_{2,{\rm g}}(x,Q^2)\otimes g(x,Q^2)\;,
\end{eqnarray}
here $q_i$ and $g$ represent the quarks and gluons distributions
respectively, $q_{\rm NS}$ stands for the usual flavor nonsinglet
combination and $q_S$ stand for the flavor-singlet quark
distribution, $q_S=\sum_{i=1}^{n_f}(q_i+\bar q_i)$. Also, $n_f$
denotes the number of effectively massless flavors. $<e^2>$
represents the average squared charge, and $\otimes$ denotes the
Mellin convolution which turns into a simple multiplication in
$N$-space.

The perturbative expansion of the coefficient functions can be
written as
\begin{eqnarray}
\label{Eq:WillsonExp}  C_{2,{\rm
i}}(x,\alpha_s(Q^2))=\sum_{n=0}\left(\frac{\alpha_s(Q^2)}{4\pi}\right)^n\;C_{\rm
2,i}^{(n)}(x)~.
\end{eqnarray}
In LO, $C_{\rm 2,NS}^{(0)}(x)=\delta(x)$, $C_{\rm
2,PS}^{(0)}(x)=C_{{\rm {2,}}g}^{(0)}(x)=C_{\rm 2,PS}^{(1)}(x)=0$
and the singlet-quark coefficient function is decomposed into the
nonsinglet and pure singlet contribution, $C_{\rm 2,q}^{(n)}\equiv
C_{\rm 2,S}^{(n)}=C_{\rm 2,NS}^{(n)}+C_{\rm 2,PS}^{(n)}$. The
coefficient functions $C_{\rm 2,i}^{(n)}$ up to N$^3$LO have been
given in \cite{Vermaseren:2005qc}.

 The nonsinglet structure function $F_{2,NS}(x,Q^2)$ up to N$^3$LO and for three
active (light) flavors has the representation
% with the nonsinglet contribution for three active
%(light) flavors being given by
\begin{eqnarray}
\label{Eq:F2-NS}x^{-1}F_{2,{\rm NS}}(x,Q^2)
&=&\Big[C_{2,q}^{(0)}+a_sC_{2,{\rm NS}}^{(1)}+a_s^2C_{2,{\rm
NS}}^{(2)+}+a_s^3C_{2,{\rm NS}}^{(3)+} \Big]\otimes \left[
\frac{1}{18}\, q_8^+ +\frac{1}{6}\, q_3^+\right](x,Q^2)~.
\end{eqnarray}
The flavor-singlet and gluon contributions in
Eq.~(\ref{Eq:F2-tot}) reads
\begin{eqnarray}
\label{Eq:F2-S}x^{-1}F_{2,S}(x,Q^2)&=&\frac{2}{9}  \left[
C_{2,q}^{(0)} +
a_sC_{2,q}^{(1)}+a_s^2C_{2,q}^{(2)}+a_s^3C_{2,q}^{(3)}\right]
\otimes \Sigma (x,Q^2)~;
\end{eqnarray}
\begin{eqnarray}
\label{Eq:F2-g} x^{-1} F_{2,g}(x,Q^2)&=&\frac{2}{9}
\left[a_sC_{2,g}^{(1)}+a_s^2C_{2,g}^{(2)}+a_s^3C_{2,g}^{(3)}\right]
\otimes g (x,Q^2)~.
\end{eqnarray}
The symbol $\otimes$ denotes the Mellin convolution
\begin{equation}
[A\otimes B](x)=\int_0^1 dx_1\int_0^1
dx_2\;\delta(x-x_1x_2)~A(x_1)B(x_2)\;.
\end{equation}
In Eq.~(\ref{Eq:F2-NS}) $q_3^+=u+\bar{u}-(d+\bar{d})=u_v-d_v$ and
$q_8^+= u+\bar{u}+d+\bar{d}-2(s+\bar{s}) =
u_v+d_v+2\bar{u}+2\bar{d}-4\bar{s}$, where $s=\bar{s}$. Also in
Eq.~(\ref{Eq:F2-S}),
$\Sigma(x,Q^2)\equiv\Sigma_{q=u,d,s}(q+\bar{q})=u_v+d_v+2\bar{u}+2\bar{d}+2\bar{s}$.
Notice that in the above equations
$a_s=a_s(Q^2)\equiv\alpha_s(Q^2)/4\pi$ denotes the strong coupling
constant and $C_{i,j}$ are the Wilson coefficients
\cite{Vermaseren:2005qc}.

The combinations of parton densities in the nonsinglet regime and
the valence region $x\geq0.3$ for $F_2^p$ in LO is
\begin{equation}
\frac{1}{x}\,F_2^{p}(x,Q^2) = \left[ \frac{1}{18}\, q_{{\rm
NS,}8}^+ +\frac{1}{6}\, q_{{\rm NS,}3}^+
\right](x,Q^2)+\frac{2}{9} \Sigma(x,Q^2)~,
\end{equation}
where $q_{{\rm NS,}3}^+=u_v-d_v$, $q_{{\rm NS,}8}^+=u_v+d_v$ and
$\Sigma=u_v+d_v$, since sea quarks can be neglected in the region
$x\geq0.3$. So in the $x$-space we have
\begin{eqnarray}
  F_2^{p}(x,Q^2) = \left (\frac{5}{18}\, x\, q_{{\rm NS,}8}^+
  + \frac{1}{6}\, x\, q_{{\rm NS,}3}^+\right) (x,Q^2)& = &
       \frac{4}{9}\,  x\, u_v(x,Q^2)+\frac{1}{9}\, x\, d_v(x,Q^2)~.
\end{eqnarray}
In the above region the combinations of parton densities for
$F_2^d$ are also given by
\begin{eqnarray}
  F_2^{d}(x,Q^2) = \left (\frac{5}{18}\, x\, q_{{\rm NS,}8}^+\right) (x,Q^2)& = &
      \frac{5}{18}\, x(u_v+d_v)(x,Q^2)~,
\end{eqnarray}
where $q_{{\rm NS,}3}^+=u_v-d_v$ and $F_2^d=(F_2^p+F_2^n)/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 \cite{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 $F_2^d=(F_2^p+F_2^n)/2+\delta F_2^d$. To obtain the
$\delta F_2^d$ we need to know the generalized vector meson
dominance (VMD) and parton mechanism at low and large values of
$Q^2$ respectively. We found that the value of $\delta F_2^d$ is
important but in low values of $x$. For example this correction
value at $Q^2$=10 GeV$^2$ and for $x>0.1$ is too small ($\sim
10^{-4}$).
 So in the valence region of this analysis, this effect
is negligible in large $x$ and we can use the
$F_2^d=(F_2^p+F_2^n)/2$ approximately.


 In the region $x
\leq 0.3$ for the difference of the proton and deuteron data we
use
\begin{eqnarray}
F_2^{NS}(x,Q^2)&\equiv& 2 (F_2^{p}-F_2^{d})(x,Q^2)\nonumber \\
&=&\frac{1}{3}\, x\, q_{{\rm NS,}3}^+(x,Q^2)= \frac{1}{3}\,
x(u_v-d_v)(x,Q^2)
  +\frac{2}{3}\, x(\bar{u}-\bar{d})(x,Q^2)~,
\end{eqnarray}
where now $q_{{\rm NS,}3}^+ = u_v-d_v+2(\bar{u}-\bar{d})$ since
sea quarks cannot be neglected for $x$ smaller than about 0.3.

The first clear evidence for the flavor asymmetry combination of
light parton distributions $x(\overline{d}-\overline{u})$ in
nature came from the analysis of NMC at CERN to study of the
Gottfried sum rule \cite{Amaudruz:1991at}. In our calculation we
supposed the $\bar{d}-\bar{u}$ distribution
\cite{Gluck:2006yz,Blumlein:2006be,ref19,Blumlein:2004pr}
%eq.(4)
\begin{equation}
x(\bar{d}-\bar{u})(x,Q_0^2) = 1.195 x^{1.24}(1-x)^{9.10}
  (1 + 14.05x - 45.52 x^2)~,
  \label{Asy}
\end{equation}
at $Q_0^2=4$ GeV$^2$ which gives a good description of the
Drell-Yan dimuon production data \cite{E866}. In this analysis,
like other analyses
\cite{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\leq 1$ can transform
to Mellin-$N$ space  by using the Mellin transform to obtain the
moments of the structure function as $\frac{1}{x}F_2^k$,
\begin{equation}
F_2^k(N,Q^2)
\equiv{\bf{M}}[F_2^k,N]=\int_0^1dx~x^{N-1}\frac{1}{x}F_2^k(x,Q^2)~,
  \label{Asy}
\end{equation}
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.~(\ref{Eq:F2-NS},\ref{Eq:F2-S},\ref{Eq:F2-g}) reduces to a
simple product

\begin{equation}
{\bf{M}}[A\otimes B,N]={\bf{M}}[A,N]{\bf{M}}[B,N]=A(N)B(N)
  \label{Momconv}
\end{equation}

By using the solution of the nonsinglet evolution equation for the
parton densities to 4$-$ loop order, the nonsinglet structure
functions are given by  \cite{Blumlein:2006be}
\begin{eqnarray}
F_2^k(N,Q^2)&=&\left(1+a_s\;C_{2,{\rm
NS}}^{(1)}(N)+a_s^2\;C_{2,{\rm NS}}^{(2)}(N)+a_s^3\;C_{2,{\rm
NS}}^{(3)}(N)\right) F_2^k(N,Q_0^2)
 \nonumber
\\
&&\times\left(\frac{a_s}{a_0}\right)^{-\hat{P}_0(N)/{\beta_0}}\Biggl\{1
- \frac{1}{\beta_0} (a_s - a_0) \left[\hat{P}_1^+(N)
- \frac{\beta_1}{\beta_0} \hat{P}_0(N) \right] \nonumber\\
& & - \frac{1}{2 \beta_0}\left(a_s^2 - a_0^2\right)
\left[\hat{P}_2^+(N) - \frac{\beta_1}{\beta_0} \hat{P}_1^+(N) +
\left( \frac{\beta_1^2}{\beta_0^2} - \frac{\beta_2}{\beta_0}
\right) \hat{P}_0(N)   \right] \nonumber\\ & & + \frac{1}{2
\beta_0^2} \left(a_s - a_0\right)^2 \left(\hat{P}_1^+(N) -
\frac{\beta_1}{\beta_0} \hat{P}_0(N) \right)^2 \nonumber\\ & &
%---
- \frac{1}{3 \beta_0} \left(a_s^3 - a_0^3\right)
\Biggl[\hat{P}_3^+(N) - \frac{\beta_1}{\beta_0} \hat{P}_2^+(N) +
\left(\frac{\beta_1^2}{\beta_0^2} - \frac{\beta_2}{\beta_0}\right)
\hat{P}_1^+(N) \nonumber\\ & & +\left(\frac{\beta_1^3}{\beta_0^3}
-2 \frac{\beta_1 \beta_2}{\beta_0^2} + \frac{\beta_3}{\beta_0}
\right) \hat{P}_0(N)  \Biggr] \nonumber\\ & & + \frac{1}{2
\beta_0^2} \left(a_s-a_0\right)\left(a_0^2 - a_s^2\right)
\left(\hat{P}_1^+(N)-\frac{\beta_1}{\beta_0} \hat{P}_0(N) \right)
\nonumber\\ & & \times \left[\hat{P}_2(N) -
\frac{\beta_1}{\beta_0} \hat{P}_1(N) -
\left(\frac{\beta_1^2}{\beta_0^2} - \frac{\beta_2}{\beta_0}
\right) \hat{P}_0(N)
 \right]
\nonumber\\ && - \frac{1}{6 \beta_0^3} \left(a_s-a_0\right)^3
\left(\hat{P}_1^+(N)-\frac{\beta_1}{\beta_0} \hat{P}_0(N)
\right)^3 \Bigg\}~.
\end{eqnarray}
Here $a_s(=\alpha_s/4\pi)$ and  $a_0$ denotes the strong coupling
constant in the scale of $Q^2$ and $Q_0^2$ respectively. $k=p,d$
and $NS$ also denotes the three above cases, i.e. proton, deuteron
and nonsinglet structure function. $C_{2,NS}^{(m)}(N)$ are the
nonsinglet Wilson coefficients in ${\it{O}}(a_s^m)$ which can be
found in \cite{FP,NS2,Vermaseren:2005qc} and $\hat{P}_m$ denote
also the Mellin transforms of the $(m+1)-$ loop splitting
functions.


\section{\label{sec:QCD Fit}Pad\'e approximations and 4-loop anomalous dimensions}

In spite of the unknown 4-loop anomalous dimensions, one can
obtain the nonsinglet parton distributions and $\Lambda_{QCD}$ by
estimating uncalculated fourth-order corrections to the nonsinglet
anomalous dimension. On the other hand the 3--loop Wilson
coefficients are known \cite{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\'e approximations.

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


For this purpose we use the Pad\'e approximations of the
perturbative series, discussed in detail for QCD, e.g., in
Refs.~\cite {Pade1,Pade2,Pade3}. Pad\'e approximations
\cite{Baker,BenderOrszag} are rational functions chosen to equal
the perturbative series to the order calculated:
\begin{eqnarray}
 [{\cal N}/{\cal M}]=\frac{a_0 + a_1x + ...
+a_{\cal N}x^{\cal N}}{1 + b_1x + ... + b_{\cal M}x^{\cal M}}~,
\label{PadeDef}
\end{eqnarray}
to the series
\begin{eqnarray}
 \label{Spade}
 S=S_0+S_1x+...+S_{{\cal N}+{\cal M}}x^{{{\cal N}+{\cal M}}}~,
\end{eqnarray}
where we set
\begin{eqnarray}
 \label{formII pade}
[{\cal N}/{\cal M}] = S + {\cal O}(x^{{\cal N}+{\cal M}+1})~,
\phantom{aaa} \nonumber
\end{eqnarray}
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 $a_s$ then
may be written in the form
\begin{eqnarray}
 \label{Anom pade}
\gamma(N)=\sum_{l=0}^{\infty}a_s^{l+1}\gamma^{(l)}(N)~.
\phantom{aaa}
\end{eqnarray}
In Mellin-$N$ space and by using this approach we can replace
$\gamma(N)$ by a rational function in $a_s$
\cite{Vermaseren:2005qc},
\begin{eqnarray} \label{eq:pade}
 \widetilde{\gamma}^{\:[\cal{N}/\cal{M}]}(N) \equiv [{\cal{N}/\cal{M}}](N) =
 \frac{p_{0} + a_s p_{1}(N) + \ldots + a_s^{\cal N} p_{\cal N}(N)}
 {1 + a_s q_{1}(N) + \ldots + a_s^{\cal M} q_{\cal M}(N)} \:\: .
\end{eqnarray}
Here ${\cal M} \,\geq\, 1$ and ${\cal N} + {\cal M} \, = \, n$,
where $n$ stands for the maximal order in $a_s$ at which the
expansion coefficients $\gamma^{(n)}(N)$ have been determined from
an exact calculation. The functions $p_i(N)$ and $q_{\!j}(N)$ are
determined from these known coefficients by expanding
Eq.~(\ref{eq:pade}) in powers of $a_s$. This expansion then also
provides the $[\cal{N}/\cal{M}]$ Pad\'e approximate for the
($n\!+\!1$)-th order quantities $\gamma^{(n\!+\!1)}$.

In this way it is easy to obtain the following results for ${\cal
M}={\cal N}=1$ and for  ${\cal M}=0,{\cal N}=2$
\begin{eqnarray}
 \label{pade11}
 \widetilde{\gamma}^{\:[1/1]}(N)&\equiv&[1/1](N)=\frac{\gamma^{(2)^2}(N)}{\gamma^{(1)}(N)}~,\nonumber
 \\
 \widetilde{\gamma}^{\:[0/2]}(N)&\equiv&[0/2](N)=\frac{2\gamma^{(1)}(N)\gamma^{(2)}(N)}
 {\gamma^{(0)}(N)}-\frac{\gamma^{(1)^3}(N)}{\gamma^{(0)^2}(N)}~.
\end{eqnarray}

The strong coupling constant $a_{s}$ plays a more central role in
the present paper to the evolution of parton densities. At
$N^{m}LO$ the scale dependence of $a_{s}$ is given by
\begin{eqnarray}
 \label{as-eqn}
  \frac{d\, a_{s}}{d \ln Q^2} \; = \; \beta_{N^mLO}(a_{s})
  \; = \; - \sum_{k=0}^m \, a_{s}^{k+2} \,\beta_k \;.
\end{eqnarray}
The expansion coefficients $\beta_k$ of the $\beta$-function of
QCD are known up to $k=3$, i.e., N$^3$LO
\cite{Tarasov:1980au,Larin:1993tp}
\begin{eqnarray}
 \label{beta-exp}
  \beta_0 &=& 11-2/3~n_f\;,
  \nonumber \\
  \beta_1 &=& 102-38/3~n_f\;,
  \nonumber \\
  \beta_2 &=& 2857/2-5033/18~n_f+325/54~n_f^2\;,
  \nonumber \\
  \beta_3 &=& 29243.0 - \: 6946.30~n_f + 405.089~n_f^2
                  + 1093/729~n_f^3~,
\end{eqnarray}
here $n_f$ stands for the number of effectively massless quark
flavors and $\beta_k$ denote the coefficients of the usual
four-dimensional $\overline{MS}$ beta function of QCD. In complete
4-loop approximation and using the $\Lambda$-parametrization, the
running coupling is given by \cite{Vogt:2004ns,Chetyrkin:1997sg}:
\begin{eqnarray}
\label{as-exp1}
 a_s(Q^2)
 &=& \frac{1}{\beta_0{L_{\Lambda}}}  -
    \frac{1}{(\beta_0{L_{\Lambda}})^2}~b_1 \ln {L_{\Lambda}} \nonumber\\  &+&
    \frac{1}{(\beta_0{L_{\Lambda}})^3} \left[b_1^2 \left(\ln^2 {L_{\Lambda}}-\ln {L_{\Lambda}}-1
    \right) + b_2\right] \nonumber\\
 &+& \frac{1}{(\beta_0{L_{\Lambda}})^4}\left[b_1^3 \left(-\ln^3 {L_{\Lambda}} +
    \frac{5}{2} \ln^2{L_{\Lambda}} +2 \ln{L_{\Lambda}}-\frac{1}{2}\right) - 3b_1 b_2\ln{L_{\Lambda}}
    +\frac{b_3}{2} \,\right]~,
\end{eqnarray}
where $L_{\Lambda}\equiv ln (Q^2/\Lambda^2)$, $b_k\equiv \beta_k
/\beta_0$, and $\Lambda$ is the QCD scale parameter. The first
line of Eq.~(\ref{as-exp1}) 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~(\ref{as-exp1}) solves the
evolution equation (\ref{as-eqn}) only up to higher orders in
$1/L_{\Lambda}$. The functional form of $\alpha_s(Q^2)$, in 4-loop
approximation and for 6 different values of $\Lambda$, is
displayed in Fig.~\ref{fig:0}. The slope and dependence on the
actual value of $\Lambda$ is especially pronounced at small $Q^2$,
while at large $Q^2$ both the energy dependence and the dependence
on $\Lambda$ becomes increasingly feeble. To be able to compare
with other measurements of $\Lambda$ we adopt the matching of
flavor thresholds at $Q^2=m_c^2$ and $Q^2=m_b^2$ with $m_c=1.5$
GeV and $m_b=4.5$ GeV as described in \cite{BAR,R1998}.




\section{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 \cite{parisi}.

According to this method, one can relate the $F_2$ structure
function with its Mellin moments
\begin{eqnarray} F_{2}^{~k,N_{max}}(x,Q^2)&=&x^{\beta}(1-x)^{\alpha}
\sum_{n=0}^{N_{max}}\Theta_n ^{\alpha,
\beta}(x)\sum_{j=0}^{n}c_{j}^{(n)}{(\alpha ,\beta )}
F_{2}^k(j+2,Q^2), \label{eg1Jacob} \end{eqnarray} where $N_{max}$
is the number of polynomials, $k$ denotes the three  cases, i.e.
$k=p,d, NS$.  Jacobi polynomials of order $n$
\cite{Parisi:1978jv}, $\Theta_n ^{\alpha, \beta}(x)$, satisfy the
orthogonality condition with the weight function $w^{\alpha
\beta}=x^{\beta}(1-x)^{\alpha}$
\begin{equation}
\int_{0}^{1}dx\;w^{\alpha \beta} \Theta_{k} ^{\alpha , \beta}(x)
\Theta_{l} ^{\alpha , \beta}(x)=\delta_{k,l}\ .\label{e8}
\end{equation}
In the above, $c_{j}^{(n)}{(\alpha ,\beta )}$ are the coefficients
expressed through $\Gamma$-functions and satisfying the
orthogonality relation in Eq.~(\ref{e8}) and $F_{2}(j+2,Q^2)$ are
the moments determined in the previous section. $N_{max}$,
$\alpha$ and $\beta$  have to be chosen so as to achieve the
fastest convergence of the series on the right-hand side of
Eq.~(\ref{eg1Jacob}) and to reconstruct $F_2$ with the required
accuracy. In our analysis we use $N_{max}=9$, $\alpha=3.0$ and
$\beta=0.5$. The same method has been applied to calculate the
nonsinglet structure function $xF_3$ from their moments
\cite{Kataev:1997nc,Kataev:1998ce, Kataev:1999bp,Kataev:2001kk}
and for polarized structure function $xg_1$ \cite{Atashbar
Tehrani:2007be,Leader:1997kw,Khorramian:2007gu}. Obviously the
$Q^2$-dependence of the polarized structure function is defined by
the $Q^2$-dependence of the moments.



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

%\section{\label{sec:QCD Fit}The  Procedure of the QCD Fits of  $F_2$ Data}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In the present analysis we choose the following parametrization
for the valence quark densities in the input scale of $Q^2_0=4$
GeV$^2$
\begin{eqnarray}
\label{equ:param} x q_v(x,Q^2_0) = {{\cal N}}_q~x^{a_q}(1-x)^{b_q}
(1 + c_q \sqrt{x} + d_q~x)~,
\end{eqnarray}
where $q=u,d$ and the normalization factors ${{\cal N}}_u$ and
${{\cal N}}_d$ are fixed by $\int_0^1 u_v dx=2$ and $\int_0^1 d_v
dx=1$, respectively.  By QCD fits of the world data for
$F_2^{p,d}$, 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~\cite{BCDMS}, SLAC~\cite{SLAC}, NMC~\cite{NMC},
H1~\cite{H1}, and ZEUS~\cite{ZEUS}. In our QCD analysis we use
three data samples~: $F_2^p(x,Q^2)$, $F_2^d(x,Q^2)$ in the
nonsinglet regime and the valence quark region $x \geq 0.3$ and
$F_2^{NS} = 2 (F_2^p - F_2^d)$ 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 $Q^2 >
4~\GeV^2$ were included in the analysis and a cut in the hadronic
mass of $W^2\equiv (\frac{1}{x}-1)\, Q^2+m_{\rm N}^2
> 12.5~\GeV^2$ 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 $F_2^p$, 232 for $F_2^d$, and
208 for $F_2^{NS}$. By considering the additional cuts on the
BCDMS ($y> 0.35$) and on the NMC data($Q^2 > 8$ GeV$^2$) the total
number of data points available for the analysis reduce from 762
to 551, because we have 227 data points for $F_2^p$, 159 for
$F_2^d$, and 165 for $F_2^{NS}$.

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 $\chi^2$
\cite{Stump:2001gu,Khorramian:2008yh}
\begin{eqnarray}
\chi _{\mathrm{global}}^{2} &=&
\sum_{n} w_{n} \chi _{n}^{2}\;,\qquad (n\;%
\mbox{labels the different experiments})\nonumber
\label{eq:Chi2global}
\\
\chi _{n}^{2} &=&\left(\frac{1-{\cal N}_{n}}{\Delta{\cal
N}_{n}}\right)^{2} +\sum_{i}\left( \frac{{\cal
N}_{n}F_{2,i}^{data}-F_{2,i}^{theor}}{{\cal N}_{n}\Delta
F_{2,i}^{data}} \right)^{2}\;. \label{eq:Chi2n}
\end{eqnarray}

For the $n^{\mathrm{th}}$ experiment, $F_{2,i}^{data}$, $\Delta
F_{2,i}^{data}$, and $%
F_{2,i}^{theor}$ denote the data value, measurement uncertainty
(statistical and systematic combined) and theoretical value for
the $i^{\mathrm{th}}$ data point. ${\Delta{\cal N}_{n}}$ is the
experimental normalization uncertainty and ${\cal N}_{n}$ is an
overall normalization factor for the data of experiment $n$. The
factor $w_{n}$ is a possible weighting factor (with default value
1).  However, we allowed for a relative normalization shift ${\cal
N}_{n}$ between the different data sets within the normalization
uncertainties ${\Delta{\cal N}_{n}}$ quoted by the experiments.
For example the normalization uncertainty of the NMC(combined)
data is estimated to be 2.5\%.
 The normalization shifts ${\cal N}_{n}$ were fitted
once and then kept fixed.

Now the sums in $\chi _{\mathrm{global}}^{2}$ run over all data
sets and in each data set over all data points. The minimization
of the  above $\chi^2$ value to determine the best parametrization
of the unpolarized parton distributions is done using the program
{\tt MINUIT} \cite{MINUIT}.

The  one $\sigma$ error for the parton density $xq_v$ as given by
Gaussian error propagation is \cite{Blumlein:2006be}
\begin{eqnarray}
\sigma( xq_v(x))^2 =  \sum_{i=1}^{n_p}\sum_{j=1}^{n_p}
                \left( \frac{\partial xq_v}{\partial p_i}\right)
                \left(\frac{\partial xq_v}{\partial p_j} \right)
                \textrm{cov}(p_i,p_j)~,
\end{eqnarray}
 where the sum runs over all fitted parameters. The functions
$\partial xq_v /
\partial p_i$ are the derivatives of $xq_v$ with respect to the fit parameter $p_i$,
and $\textrm{cov}(p_i,p_j)$ are the elements of the covariance
matrix. The derivatives $\partial xq_v / \partial p_i$ can be
calculated analytically at the input scale $Q_0^2$. Their values
at $Q^2$ 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 $Q_0^2$ for the parton density
parameterized as in Eq.~(\ref{equ:param}) is given by
\begin{eqnarray}
q_{v}(N,a_{q},b_{q},c_{q},d_{q})&=&{\cal{N}}_{q}~{\bf{M}}(n,a_{q},b_{q},c_{q},d_{q})~,
\end{eqnarray}
with the normalization constant
\begin{equation}
{\cal{N}}_{q}=
\frac{C_{q_v}}{{\bf{M}}(1,a_{q},b_{q},c_{q},d_{q})}~.
\end{equation}
Here $C_{q_v}$ is the respective number of valence quarks, i.e.
$C_{u_v}$=2 and $C_{d_v}$=1. In the above
${\bf{M}}(n,a_{q},b_{q},c_{q},d_{q})$  is given by
\begin{eqnarray}
{\bf{M}}(n,a_{q},b_{q},c_{q},d_{q})=B[a_{q}+n-1,b_{q}+1]+c_{q}B[a+n+1/2,b+1]+d_{u}B[a_{q}+n,b_{q}+1]~,\nonumber \\
\end{eqnarray}
where $B[a,b]$ denotes the Euler beta function for complex
arguments. The general form of the derivative of the Mellin moment
$q_{v}$ with respect to the parameter $p$ is given by
\begin{equation}
\frac{\partial q_{v}(N,p)}{\partial p}={\bf{M}}(n,p)\frac{\partial {\cal{N}}_{q}}{%
\partial p}+{\cal{N}}_{q}\frac{\partial {\bf{M}}(n,p)}{\partial
p}~.
\end{equation}
In this analysis only the parameters $a_q$ and  $b_q$ have been
fitted for both the $xu_v$ and $xd_v$ 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 $u_v$ and $d_v$ parton densities with respect
to parameter $a_q$ and $b_q$. For example:

\begin{eqnarray}
f(n,a_q) &\equiv&\frac{\partial {\bf{M}}(n,a_q)}{\partial a_q}=B[a_{q}+n-1,b_{q}+1](%
\psi [a_{q}+n-1]-\psi [a_{q}+b_{q}+n])+ \nonumber\\
&&c_{q}B[a_{q}+n-1/2,b_{q}+1](\psi [a_{q}+n-1/2]-\psi
[a+b+n+1/2])+
\nonumber \\
&&d_{q}B[a_{q}+n,b_{q}+1](\psi [a_{q}+n]-\psi [a_{q}+b_{q}+n+1])~,
\end{eqnarray}


\begin{eqnarray}
f(n,b_q) &\equiv&\frac{\partial {\bf{M}}(n,b_q)}{\partial b_{q}}=B[a_{q}+n-1,b_{q}+1](%
\psi [b_{q}+1]-\psi [a_{q}+b_{q}+n])+  \nonumber\\
&&c_{q}B[a_{q}+n-1/2,b_{q}+1](\psi [1+b_{q}]-\psi
[a_{q}+b_{q}+n+1/2])+
\nonumber \\
&&d_{q}B[a_{q}+n,b_{q}+1](\psi [b_{q}+1]-\psi [a_{q}+b_{q}+n+1])~,
\end{eqnarray}
and now we can reach the below derivatives for $u_v(N)$ and
$d_v(N)$ with respect to parameters $a_q$ and $b_q$
\begin{equation}
\frac{\partial q_{v}(N,p)}{\partial p}%
={\cal{N}}_{q}\left(f(n,p)-f(1,p){\bf{M}}(n,p)/{\bf{M}}(1,p)\right)~,
\end{equation}
also $\psi [n]={d~ln~\Gamma(n)}/{dn}$ is Euler's  $\psi$-function.

To obtain the error calculation of the structure functions
$F_2^p$,  $F_2^d$ , and  $F_2^{NS}$ 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 $\Lambda$ is fixed and regarded uncorrelated. The error
calculation for a variable $\Lambda$  is done numerically due to
the non�linear relation and required iterative treatment in the
calculation of $\alpha_s(Q^2,\Lambda)$
\cite{Khorramian:2008yh,Blumlein:2006be}.

\section{Results}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In the QCD analysis of the present paper we used three data sets:
the structure functions $F_2^{p}(x,Q^2)$  and $F_2^{d}(x,Q^2)$ in
the region  of $x \geq 0.3$ and the combination of these structure
functions $F_2^{\rm NS}(x,Q^2)$ in the region  of $x < 0.3$~.
Notice that we take into account the cuts $Q^2>4$ GeV$^2$,
$W^2>12.5$ GeV$^2$ for our QCD fits to determine some unknown
parameters. In Fig.(\ref{fig:1}) the proton, deuteron and
nonsinglet data for $F_2^p(x,Q^2)$, $F_2^d(x,Q^2)$ and
$F_2^{NS}(x,Q^2)$ are shown in the nonsinglet regime and the
valence quark region $x \geq 0.3$ indicating the above cuts by a
vertical dashed line. The solid lines correspond to the N$^3$LO
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 $Q^{2}\rightarrow\infty$). At intermediate
and low $Q^{2}$ this assumption will begin to break down and the
moments will be subject to potentially significant power
corrections, of order ${{\cal O}}~(m_{N}^{2}/Q^{2})$, where $m_N$
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
\cite{Georgi:1976ve,Gluck:2006yz}
\begin{eqnarray}
\label{equ:TMC} F_{2,{\rm TMC}}^{k}(n,Q^2) & \equiv &
  \int_0^1 x^{n-1} \frac{1}{x}F_{2,{\rm TMC}}^{k}(x,Q^2)\, dx  \nonumber\\
& = &F_{2}^{
k}(n,Q^2)+\frac{n(n-1)}{n+2}\left(\frac{m_N^2}{Q^2}\right)\,F_{2}^{
k}(n+2,Q^2)\nonumber \\
&&+\frac{(n+2)(n+1)n(n-1)}{2(n+4)(n+3)}\left(\frac{m_N^2}{Q^2}\right)^2\,F_{2}^{
k}(n+4,Q^2)+ {\cal{O}}\left( \frac{m_N^2}{Q^2}\right)^3~,
\end{eqnarray}
where higher powers than $(m_{\rm N}^2/Q^2)^2$ are negligible for
the relevant $x <0.8$ region. By inserting Eq.~(\ref{equ:TMC}) in
Eq.~(\ref{eg1Jacob}) we have
\begin{eqnarray}
F_{2}^{~k, N_{max}}(x,Q^2)&=&x^{\beta}(1-x)^{\alpha}
\sum_{n=0}^{N_{max}}\Theta_n ^{\alpha,
\beta}(x)\times\sum_{j=0}^{n}c_{j}^{(n)}{(\alpha ,\beta )}
F_{2,{\rm TMC}}^{k}(j+2,Q^2)\;, \label{eg1JacobTMC} \end{eqnarray}
where $F_{2,{\rm TMC}}^{k}(j+2,Q^2)$ are the moments determined by
Eq.~(\ref{equ:TMC}). In Fig.(\ref{fig:1})  the dashed lines
correspond to the N$^3$LO QCD fit adding target mass corrections.


Despite the kinematic cuts ($Q^2\geq 4$ GeV$^2$, $\, W^2\equiv
(\frac{1}{x}-1)\, Q^2+m_{\rm N}^2\geq 12.5$ GeV$^2$) used for our
analysis, we also take into account higher twist corrections to
$F_2^{p}(x,Q^2)$ and $F_2^{d}(x,Q^2)$ in the kinematic region $Q^2
\geq 4 \GeV^2, 4<W^2 < 12.5 \GeV^2$ in order to learn whether
nonperturbative effects may still contaminate our perturbative
analysis. For this purpose we extrapolate the QCD fit results
obtained for $W^2 \geq 12.5 \GeV^2$ to the region $Q^2 \geq 4
\GeV^2, 4<W^2 < 12.5 \GeV^2$ and from the difference between data
and theory, applying target mass corrections in addition. Now by
considering higher twist correction
%-----------------------------------------------------------------------
\begin{eqnarray}
\label{equ:HTC} F_2^{\rm exp}(x,Q^2) = O_{\rm TMC}[F_2^{\rm
HT}(x,Q^2)] \cdot \left( 1 +
\frac{h(x,Q^2)}{Q^2[\GeV^2]}\right)\;,
\end{eqnarray}
%-----------------------------------------------------------------------
the higher twist coefficient can be extract. Here the operation
$O_{\rm TMC}[...]$ denotes taking the target mass corrections of
the twist--2 contributions to the respective structure function.
The coefficients $h(x,Q^2)$ are determined in bins of $x$ and
$Q^2$ and are then averaged over $Q^2$.  We extrapolate our QCD
fits to the region $12.5 \GeV^2 \geq W^2 \geq 4 \GeV^2$ in
Fig.(\ref{fig:1}). The dash-dotted lines in this figure
  correspond to the N$^3$LO QCD fit adding target mass and higher twist corrections.
 There, at higher values of $x$ a clear gap between the
data and the QCD fit is seen.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\


\begin{table*}
\begin{tabular}{|c|c|c|c|c|c|}
\hline \hline
           &              &     NLO     &     N$^{2}$LO   &    N$^{3}$LO Pad\'e~[1/1]   &    N$^{3}$LO Pad\'e~[0/2]  \\
\hline $u_v$      & $a_u$          &  0.7434
$\pm$ 0.009 &  0.7772 $\pm$ 0.009 & 0.79167 $\pm$ 0.0106  & 0.79176 $\pm$ 0.0099\\
           & $b_u$         &  3.8907 $\pm$ 0.040 &
 4.0034 $\pm$ 0.033  &  4.02637 $\pm$ 0.0402  & 4.02685 $\pm$ 0.0327    \\
           & $c_u  $                 &  0.1620             &
 0.1000    & 0.0940    &  0.0940         \\
           & $d_u$             & 1.2100             &
1.1400         & 1.1100    &  1.1100     \\
\hline $d_v$      & $a_d$         &  0.7369
$\pm$ 0.040 &  0.7858 $\pm$ 0.043  &  0.80927 $\pm$ 0.0621  & 0.80927 $\pm$ 0.0407  \\
           & $b_d$          &  3.5051 $\pm$ 0.225 &
 3.6336 $\pm$ 0.244   & 3.76847 $\pm$ 0.3499 & 3.76858 $\pm$ 0.2278 \\
           & $c_d  $              & 0.3899             &
0.1838       & 0.1399  & 0.1399      \\
           & $d_d$             & -1.3700             &
-1.2152     & -1.1200   &  -1.1200     \\
\hline \multicolumn{2}{|c|}{$\Lambda_{\rm QCD}^{\rm N_f=4}$, MeV}
&
263.8 $\pm$ 30 & 239.9 $\pm$ 27  &  241.44 $\pm$ 29 &  241.45 $\pm$ 27 \\
\hline \hline \multicolumn{2}{|c|}{$\chi^2 / ndf$} & 523/546 =
0.9578 &
506/546 = 0.9267 & 491.07/546 = 0.8994   &  491.12/546 = 0.8995 \\
\hline \hline
\end{tabular}
\caption{\label{tab:param}{\sf Parameter values of the NLO,
N$^2$LO from Ref.~\cite{Khorramian:2008yh} and N$^3$LO nonsinglet
QCD fit at $Q_0^2 = 4~ \mbox{GeV}^2$ for Pad\'e [1/1] and Pad\'e
[0/2].}}
%\begin{ruledtabular}
%\end{ruledtabular}
\end{table*}
\vspace{1 cm}

In Table (\ref{tab:param}) we summarize the NLO, N$^2$LO, and
N$^3$LO with using Pad\'e [1/1] and [0/2] fit results for the
parameters of the parton densities $xu_v(x,Q^2_0)$,
$xd_v(x,Q^2_0)$ and $\Lambda_{\rm QCD}^{\rm N_f =4}$. The values
without error have been fixed after a first minimization since the
data do not constrain these parameters well enough. In this table
we also  compare the N$^3$LO  results with the NLO and N$^2$LO
results from Ref.\cite{Khorramian:2008yh}. The results show a good
compatibility between Pad\'e [1/1] and [0/2] approximations in
4--loop order. The resulted value of $\chi^2/ndf$ is 0.9578 at
NLO,  0.9267 at N$^2$LO, and 0.8994 and 0.8995 for Pad\'e [1/1]
and [0/2] respectively at N$^3$LO. Our results for the covariance
matrix of the  N$^3$LO nonsinglet QCD fit for Pad\'e [1/1] and
[0/2] are presented in Table(\ref{tab:cov}).

\renewcommand{\arraystretch}{0.9}
\begin{table*}
\begin{tabular}{|c||c|c|c|c|c|}
\hline \hline
{\bf N$^3$LO} Pad\'e[1/1] & $a_{u}$ & $b_{u}$ & $a_{d}$ & $b_{d}$& $\Lambda_{\rm QCD}^{\rm N_f=4}$ \\
\hline
 $a_{u}$ &~{\bf 1.13$\times$10$^{-4}$}&  &  &  &  \\
\hline
 $b_{u}$        &~2.35$\times$10$^{-4}$ &{\bf ~1.62$\times$10$^{-3}$} &  &  &  \\
\hline
 $a_{d}$        &~1.09$\times$10$^{-4}$ &~-1.59$\times$10$^{-3}$ & {\bf ~3.86$\times$10$^{-3}$} &  & \\
\hline
 $b_{d}$        &~1.67$\times$10$^{-4}$ &~-8.84$\times$10$^{-3}$ &~2.11$\times$10$^{-2}$ &~{\bf 1.23$\times$10$^{-1}$} & \\
\hline
 $\Lambda_{QCD}^{(4)}$        &~1.71$\times$10$^{-4}$ &~-3.49$\times$10$^{-4}$ &~5.04$\times$10$^{-4}$ ~&2.61$\times$10$^{-3}$ &~{\bf 8.65$\times$10$^{-4}$} \\
\hline \hline
{\bf N$^3$LO} Pad\'e[0/2]  & $a_{u}$ & $b_{u}$ & $a_{d}$ & $b_{d}$& $\Lambda_{\rm QCD}^{\rm N_f=4}$ \\
\hline
 $a_{u}$ & {\bf ~0.98$\times$10$^{-4}$} &  &  &  &  \\
\hline
 $b_{u}$        &~1.83$\times$10$^{-4}$ &~{\bf 1.07$\times$10$^{-3}$} &  &  &  \\
\hline
 $a_{d}$        &~-5.07$\times$10$^{-5}$ &~-6.01$\times$10$^{-4}$ &~{\bf 1.66$\times$10$^{-3}$} &  & \\
\hline
 $b_{d}$        &~-1.11$\times$10$^{-4}$ &~-3.30$\times$10$^{-3}$ &~8.58$\times$10$^{-3}$ &~{\bf 5.19$\times$10$^{-2}$} & \\
\hline
 $\Lambda_{QCD}^{(4)}$        &~1.59$\times$10$^{-4}$ &~-1.99$\times$10$^{-4}$& 1.94$\times$10$^{-4}$ &~8.07$\times$10$^{-4}$ &~{\bf 7.53$\times$10$^{-4}$} \\
\hline \hline
\end{tabular}

\caption{\label{tab:cov}{\sf Our results for the covariance matrix
of the N$^3$LO nonsinglet QCD fit for Pad\'e [1/1] and [0/2]  at
$Q_0^2 = 4~\GeV^2$ by using MINUIT\cite{MINUIT}.}}
\end{table*}



Figure~(\ref{fig:2}) illustrates our fit results for
$xu_v(x,Q^2_0)$, $xd_v(x,Q^2_0)$ at $Q_0^2 = 4 \GeV^2$ up to
N$^3$LO and for Pad\'e [1/1] with correlated errors. In this
figure our results for N$^3$LO  compared with results obtained
from \cite{Khorramian:2008yh} at LO, NLO, and N$^2$LO QCD
analysis. The shaded areas represent the fully correlated one
$\sigma$ statistical error bands.


In Fig.~(\ref{fig:3}) we show the evolution of the valence quark
distributions $xu_v(x,Q^2)$ and $xd_v(x,Q^2)$ from $Q^2 = 1
\GeV^2$ to $Q^2 = 10^4 \GeV^2$ in the region  $x\in [10^{-4},1]$
at N$^3$LO. In this figure we also compared our results with the
nonsinglet QCD analysis from \cite{Blumlein:2006be}. With rising
values of $Q^2$ the distributions flatten at large values of $x$
and rise at low values.

Another way to test the N$^3$LO fit results is comparison of low
order moments of the distributions $u_v(x,Q^2), d_v(x,Q^2),$ and
$u_v(x,Q^2)- d_v(x,Q^2)$. In Table~\ref{tab:table6} we present the
lowest nontrivial moments of these distributions at $Q^2 = Q_0^2$
in N$^3$LO and compare to the respective moments obtained for the
parameterizations \cite{Blumlein:2006be}.

 We should note that the unknown parameters
 are correlated and almost depend on the method of the QCD fits.
 We believe that the source of the small
difference between the results of our analysis and reported
results in \cite{Blumlein:2006be} is the kind of the different
method of the QCD analysis. We used the Jacobi polynomial method
as an expansion method to do QCD fits but they used the exact
inverse Mellin technique to obtain some unknown parameters. We
also found that the results of Pad\'e [1/1] and [0/2] in 4-loop
level are almost the same.

\renewcommand{\arraystretch}{0.9}
\begin{table*}
\begin{tabular}{|c|c|c|c|c|}
\hline \hline  $f$ & N & BBG~\cite{Blumlein:2006be} & N$^{3}$LO Pad\'e[1/1] & N$^{3}$LO Pad\'e[0/2] \\
\hline
$u_{v}$ & $2$ & $0.3006\pm 0.0031$ & $0.30757\pm 0.0026$ & $0.30806\pm 0.0028$\\
& $3$ & $0.0877\pm 0.0012$ & $0.08771\pm 0.0011$& $0.08781\pm 0.0012$ \\
& $4$ & $0.0335\pm 0.0006$ & $0.03320\pm 0.0006$& $0.03323\pm 0.0006$ \\
\hline
$d_{v}$ & $2$ & $0.1252\pm 0.0027$ & $0.12450\pm 0.0024$& $0.12495\pm 0.0025$ \\
& $3$ & $0.0318\pm 0.0009$ & $0.03040\pm 0.0008$& $0.03012\pm 0.0008$ \\
& $4$ & $0.0106\pm 0.0004$ & $0.00992\pm 0.0004$& $0.00993\pm 0.0005$ \\
\hline
$u_{v}-d_{v}$ & $2$ & $0.1754\pm 0.0041$ & $0.18305\pm 0.0036$& $0.18310\pm 0.0038$ \\
& $3$ & $0.0559\pm 0.0015$ & $0.05767\pm 0.0013$ & $0.05769\pm 0.0014$ \\
& $4$ & $0.0229\pm 0.0007$ & $0.02329\pm 0.0007$& $0.02329\pm
0.0007$
\\ \hline \hline
\end{tabular}
\normalsize
%
\vspace{2mm} \caption{\label{tab:table6}{\sf Comparison of low
order moments from our nonsinglet N$^3$LO QCD analysis at $Q_0^2 =
4~\mbox{GeV}^2$ with the N$^3$LO analysis from
Ref.~\cite{Blumlein:2006be}. }}
\end{table*}

To perform higher twist QCD analysis of the nonsinglet world data
in N$^3$LO, we consider the $Q^2 \geq 4 \GeV^2, 4<W^2 < 12.5
\GeV^2$ cuts. The number of data points in the above range for
proton and deuteron is 279 and  278, respectively.
 The extracted distributions for $h(x)$ in N$^3$LO are depicted in
Fig.(\ref{fig:4}) for the nonsinglet case considering scattering
off the proton and deuteron target. According to our results the
coefficient $h(x)$ grows towards large $x$.To compare, we also
present the reported results of the early N$^2$LO analysis
\cite{Khorramian:2008yh} in Fig.(\ref{fig:4}). Also in this figure
HT contributions have the tendency to decrease form N$^2$LO to
N$^3$LO. This effect was observed for the first time in the case
of fits of $F_3$ DIS $\nu N$ data in \cite{Kataev:1997nc} and then
studied in more detail in \cite{Kataev:1999bp,Kataev:2001kk}.

This similar effect was also observed in the fits of $F_2$ charge
lepton-nucleon DIS data
\cite{Yang:1999xg,Blumlein:2006be,Gluck:2006yz,Blumlein:2008kz}.
% Note that the results for $h(x)$ in LO are not
%presented in the BBG model \cite{Blumlein:2006be,Blumlein:2008kz}.
In Ref. \cite{Gluck:2006yz}, the functional form for $h(x)$ is
chosen by
\begin{equation}
h(x) = a\left( \frac{x^b}{1-x} -c\right)\,, \label{eq.hx}
\end{equation}
and it is possible to compare $h(x)$ results in N$^2$LO and
N$^3$LO . In Table \ref{tab:hx} we present our results for $a, b,
c$ in the above equation.


\begin{table*}
\begin{tabular}{cccc}
\hline\hline & $a$ & $b$ & $c$ \\ \hline Proton & $~~1.015~~$ &
$~~3.928~~$ & $~~-0.193~~$ \\ \hline Deuteron & $4.481$ & $7.759$
&
$~~-0.064~~$ \\
\hline
\end{tabular}
\vspace{2mm} \caption{\label{tab:hx}{\sf Our results for h(x)
function according to Eq.~(\ref{eq.hx}) and for N$^3$LO.}}
\end{table*}

As seen from Fig.(\ref{fig:4})  $h(x)$ is widely independent of
the target comparing the results for deeply inelastic scattering
off protons and deuterons.


\section{Discussion}
A study \cite{Khorramian:2008yh} of the available world data on
deep-inelastic lepton-proton and lepton-deuteron scattering
provided a determination of the valence quark parton densities and
$\alpha_s$ in wide ranges of the Bjorken scaling variable $x$ and
$Q^2$ up to 3-loop. In the nonsinglet case, where heavy flavor
effects are negligibly small, the analysis can be extended to
4-loop level, i.e. to QCD in N$^3$LO perturbative expansion.



 The analysis was performed
using the Jacobi polynomials method to determine the parameters of
the problem in a fit to the data. A new aspect in comparison with
previous analysis is that we determine the parton densities and
the QCD scale up to N$^3$LO by using  the Jacobi polynomial
expansion method and using Pad\'e approximations. The benefit of
this approach is the possibility to determine nonsinglet parton
distributions analytically and not numerically. In Ref.
\cite{Program:summary} we arrange the MATHEMATICA program to
extract $xu_v(x,Q^2)$ and $xd_v(x,Q^2)$ up to the 4-loops.

In this analysis we adopt the $\bar{d}-\bar{u}$ distribution at
$Q_0^2=4$ GeV$^2$ from
Refs.~\cite{Gluck:2006yz,Blumlein:2006be,ref19,Blumlein:2004pr},
which gives a good description of the Drell-Yan dimuon production
data \cite{E866}. The nonsinglet regime is manifesting itself at
$x\geq 0.1$ as the rule. In this regime, when we changed the sea
distribution from the other groups,  the value of $\chi^2$,
valence distributions, $\Lambda$ and $\alpha_s$ varied, but only
slightly. For example, we used the $\bar{d}-\bar{u}$ distribution
from
\cite{Gluck:2007ck,Martin:2007bv,JimenezDelgado:2008hf,Martin:2009iq}
and we found that the value of  $\chi^2$ varies by about 3\%
 and $\Lambda$ by about 1\%-2\%.


In the QCD analysis we parameterized the strong coupling constant
$\alpha_s$ in terms of four massless flavors determining
$\Lambda_{\rm QCD}$. Up to N$^3$LO results fitting the data, are
\begin{eqnarray}
\Lambda_{\rm QCD}^{(4)} &=& 213.2 \pm 28
\; \mbox{MeV},~~{\tt LO}, \nonumber \\
\Lambda_{\rm QCD}^{(4)} &=& 263.8 \pm 30 \;
\mbox{MeV},~~{\tt NLO}, \nonumber \\
\Lambda_{\rm QCD}^{(4)} &=& 239.9 \pm 27 \;
\mbox{MeV},~~{\tt N^2LO},\nonumber \\
\Lambda_{\rm QCD}^{(4)} &=& 241.4 \pm 29 \; \mbox{MeV},~~{\tt
N^3LO}.
\end{eqnarray}
These results can be expressed in terms of $\alpha_s(M_Z^2)$:
\begin{eqnarray}
\alpha_s(M_Z^2) &=& 0.1281 \pm 0.0028,~~{\tt LO},  \nonumber \\
\alpha_s(M_Z^2) &=& 0.1149 \pm 0.0021,~~{\tt NLO}, \nonumber \\
\alpha_s(M_Z^2) &=& 0.1131 \pm 0.0019,~~{\tt N^2LO}, \nonumber \\
\alpha_s(M_Z^2) &=& 0.1139 \pm 0.0020,~~{\tt N^3LO}.
\end{eqnarray}
Note that in above results we use the matching between $n_f$ and
$n_{f+1}$ flavor couplings calculated in Ref.
\cite{Chetyrkin:1997sg}. We adopt this prescription to be able to
compare our results with other measurement of $\Lambda_{\rm QCD}$.

 The $\alpha_s(M_Z^2)$ values can be
compared with results from other QCD analysis of inclusive
deep--inelastic scattering data in N$^2$LO
\normalsize
\begin{center}
\begin{tabular}{rcllll}
 & & & & & \\
A02   \cite{A02}:      & & $\alpha_s(M_Z^2)$=0.1143 & $\pm$0.0014 &   \\
GRS  \cite{Gluck:2006yz}:&  & $\alpha_s(M_Z^2)$=0.111 &  &
\\
 MRST03 \cite{MRST03}:     & &$\alpha_s(M_Z^2)$= 0.1153 & $\pm$0.0020 &    \\
 SY01(ep)  \cite{SY}:  & & $\alpha_s(M_Z^2)$=0.1166 & $\pm$0.0013 &                \\
 SY01($\nu$N) \cite{SY}:& & $\alpha_s(M_Z^2)$=0.1153 & $\pm$0.0063 &                \\
A06 \cite{Alekhin:2006zm}: & &  $\alpha_s(M_Z^2)$=0.1128 & $\pm$
0.0015 &
\\
BBG \cite{Blumlein:2006be}:  & & $\alpha_s(M_Z^2)$=0.1134  & \footnotesize{$\begin{array}{c} +0.0019 \\
-0.0021 \end{array}$} &\\
BM07   \cite{Brooks:2006wh}:      & & $\alpha_s(M_Z^2)$=0.1189 &
$\pm$0.0019 &
\\
KPS00($\nu N$) \cite{Kataev:1999bp}:      & &
$\alpha_s(M_Z^2)$=0.118 &$\pm$ 0.002 ~($stat$)$\pm$ 0.005~($syst$) \\
&&& $\pm$
0.003~($theory$) &  \\
KPS03($\nu N$) \cite{Kataev:2001kk}:      & &
$\alpha_s(M_Z^2)$=0.119 &$\pm$ 0.002 ~($stat$)$\pm$ 0.005~($syst$)\\
&&&  $\pm$ 0.002 ~($threshold$)
$^{+0.004}_{-0.002}$~ ($scale$) &  \\
KT08 \cite{Khorramian:2008yh}:  & & $\alpha_s(M_Z^2)$=0.1131  & \footnotesize{$\begin{array}{c} \pm 0.0019  \end{array}$} &\\
\end{tabular}
\end{center}
\normalsize The N$^3$LO values of $\alpha_s(M_Z^2)$ can also be
compared with results from other QCD analysis \noindent
\normalsize
\begin{center}
\begin{tabular}{rcllll}
 & & & & & \\
BBG \cite{Blumlein:2006be}:  & & $\alpha_s(M_Z^2)$=0.1134  & \footnotesize{$\begin{array}{c} +0.0019 \\
-0.0021 \end{array}$} &\\
\end{tabular}
\end{center}
\normalsize and with the value of the world average $0.1189 \pm
0.0010$ \cite{Bethke:2006ac} and also the current world average
%-------------------------------------------------------------------------
\begin{eqnarray}
 \alpha_s(M_Z^2) = 0.1184 \pm 0.0007~,
\end{eqnarray}
%-----------------
which has been extracted in \cite{Bethke:2009jm} very recently. It
seems that our results confirm that the value of $\alpha_s(M_Z^2)$
from DIS turns out to be sizably below the world average. In this
case, it would be useful to find out which set of data is mainly
responsible for the low value of $\alpha_s(M_Z^2)$. We will try to
see which subset makes $\alpha_s(M_Z^2)$ particularly small in a
future work.

We hope our results of QCD analysis of structure functions in
terms of Jacobi polynomials
 could be able to describe more complicated hadron structure functions. We also
hope to be able to consider massive quark contributions by using
the structure function expansion in terms of the Jacobi
polynomials.

\newpage
\begin{acknowledgments}
We are especially grateful to G. Altarelli  for guidance and
critical remarks. A.N.K. is grateful to F. Olness and J.
Bl\"umlein for useful discussions and constructive comments. H.K.
is grateful to S.~Moch for his guidance and discussion.  We would
like to thank Z. Karamloo and M. Ghominejad for reading the
manuscript of this paper. A.N.K is grateful to TH-PH division at
CERN for their hospitality while he visited there and  amended
this paper. A.N.K. thanks Semnan University for partial financial
support of this project. We acknowledge the School of Particles
and Accelerators, Institute for Research in Fundamental Sciences
(IPM) for financially supporting this project.
\end{acknowledgments}


\newpage
\begin{figure}
\vspace{2 cm}
\includegraphics[angle=0, width=8.0cm]{alphaN3LO}% Here is how to import EPS art
\caption{\label{fig:0}{\sf The strong running of $\alpha_s(Q^2)$,
according to Eq.~\ref{as-exp1}, in 4- loop approximation and for
different values of $\Lambda$. }}
\end{figure}


\newpage
\begin{figure}
\vspace{2 cm}
\includegraphics[angle=0, width=12.0cm]{F2-P-D-NS}% Here is how to import EPS art
\caption{\label{fig:1}{\sf The structure functions $F_2^p$,
$F_2^d$, and $F_2^{NS}$  as a function of $Q^2$ in intervals of
$x$. Shown are the Pad\'e [1/1] QCD fits in N$^3$LO (solid line)
and the contributions from target mass corrections (dashed line)
and higher twist (dash--dotted line). The vertical dashed line
indicates the regions with $W^2 > 12.5$ GeV$^2$.
%The vertical dashed line
%indicate the regions with $W^2 > 12.5~{\GeV^2}$.
}}
\end{figure}

\newpage
\begin{figure}
\vspace{2 cm}
\begin{center}
%\includegraphics[angle=0, width=5.0cm]{xuvxdvQ04.eps}
\includegraphics[angle=0, width=9.5cm]{xuvxdv-Q04.eps}
\end{center}
\vspace{0 cm} {\sf \caption{\label{fig:2} {\sf The parton
densities $xu_v$ and $xd_v$ up to 4-loop (Pad\'e [1/1]) at the
input scale $Q_0^2 = 4.0~{\rm{\GeV^2}}$ (solid line) compared with
results obtained from N$^2$LO analysis (dashed-- line), NLO
analysis (dash--dotted line), and LO analysis(dash--dott--dotted
line)~\cite{Khorramian:2008yh}. The shaded areas represent the
fully correlated one $\sigma$ statistical error bands.}}}
\end{figure}

\newpage
\begin{figure}
\vspace{2 cm}
\includegraphics[angle=0, width=16 cm]{xuvxdvallQ.eps}
\begin{center}
\end{center}
{\sf \caption{\label{fig:3} {\sf The parton densities $xu_v$ and
$xd_v$ at N$^3$LO evolved up to $Q^2 = 10000~{ \GeV^2}$ (solid
lines) compared with results obtained by
%A05 (dashed
%line)~\cite{Alekhin:2005gq},
BBG (dashed line)~\cite{Blumlein:2006be}
%, and MRST (dashed--dotted-dottedline)~\cite{Martin:2007bv,MRST04}
.}}}
\end{figure}
\newpage

\begin{figure}
\begin{center}
\includegraphics[angle=0, width=13.0cm]{ht2new.eps}
\end{center}
{\sf \caption{\label{fig:4} The higher twist coefficient $h(x)$
for the proton and deuteron data as a function of $x$ in N$^3$LO
(solid line) compared with results obtained by N$^2$LO (dashed
line) \cite{Khorramian:2008yh}
% and BBG (dashed--dotted line) \cite{Blumlein:2008kz}
. }}
\end{figure}
\normalsize




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Bibliography
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage

\bibliography{apssamp}% Produces the bibliography via BibTeX.
\begin{thebibliography}{99}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\cite{Altarelli:2009za}
\bibitem{Altarelli:2009za}
  G.~Altarelli,
  %``DIS 2009 Concluding Talk: Outlook and Perspective,''
  arXiv:0907.1751 [hep-ph].
  %%CITATION = ARXIV:0907.1751;%%


%\cite{CERN:2009}
\bibitem{CERN:2009}
http://indico.cern.ch/conferenceOtherViews.py?view=standard\&confId=43022.


%\cite{Martin:2009iq}
\bibitem{Martin:2009iq}
  A.~D.~Martin, W.~J.~Stirling, R.~S.~Thorne and G.~Watt,
  %``Parton distributions for the LHC,''
  arXiv:0901.0002 [hep-ph].
  %%CITATION = ARXIV:0901.0002;%%

%\cite{Thorne:2009ky}
%\bibitem{Thorne:2009ky}
%  R.~S.~Thorne, A.~D.~Martin, W.~J.~Stirling and G.~Watt,
  %``Status of MRST/MSTW PDF sets,''
%  arXiv:0907.2387 [hep-ph].
  %%CITATION = ARXIV:0907.2387;%%

%\cite{Thorne:2009me}
\bibitem{Thorne:2009me}
  R.~S.~Thorne,
  %``Summary of the HERA-LHC workshop,''
  Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 191}, 295 (2009)
  [arXiv:0902.4820 [hep-ph]].
  %%CITATION = NUPHZ,191,295;%%

%\cite{Nadolsky:2008zw}
\bibitem{Nadolsky:2008zw}
  P.~M.~Nadolsky {\it et al.},
  %``Implications of CTEQ global analysis for collider observables,''
  Phys.\ Rev.\  D {\bf 78}, 013004 (2008)
  [arXiv:0802.0007 [hep-ph]].
  %%CITATION = PHRVA,D78,013004;%%




%\cite{Khorramian:2008yh}
\bibitem{Khorramian:2008yh}
  A.~N.~Khorramian and S.~A.~Tehrani,
  %``NNLO QCD contributions to the flavor nonsinglet sector of ${\bf
  %{F_2(x,Q^2)}}$,''
  Phys.\ Rev.\  D {\bf 78}, 074019 (2008)
  [arXiv:0805.3063 [hep-ph]].
  %%CITATION = PHRVA,D78,074019;%%

%\cite{Altarelli:2008fd}
%\bibitem{Altarelli:2008fd}
%  G.~Altarelli,
  %``Status of the Standard Model at the LHC Start,''
 % arXiv:0804.4147 [hep-ph].
  %%CITATION = ARXIV:0804.4147;%%

%\cite{Altarelli:2008aj}
%\bibitem{Altarelli:2008aj}
%  G.~Altarelli, R.~D.~Ball and S.~Forte,
  %``Small x Resummation with Quarks: Deep-Inelastic Scattering,''
 % Nucl.\ Phys.\  B {\bf 799}, 199 (2008)
 % [arXiv:0802.0032 [hep-ph]].
  %%CITATION = NUPHA,B799,199;%%


%\cite{Altarelli:2008xp}
%\bibitem{Altarelli:2008xp}
%  G.~Altarelli, R.~D.~Ball and S.~Forte,
  %``Structure Function Resummation in small-x QCD,''
%  PoS {\bf RADCOR2007}, 028 (2007)
%  [arXiv:0802.0968 [hep-ph]].
  %%CITATION = POSCI,RADCOR2007,028;%%

%\cite{Barker}
\bibitem{Barker}
I. S. Barker, B. R. Martin and G. Shaw,
\newblock{\it Z. Phys.} {\bf C19} (1983)  147;
\\I. S. Barker and B. R. Martin,
\newblock{\it Z. Phys.} {\bf C24} (1984) 255;
\\S. P. Kurlovich, A. V. Sidorov and N. B. Skachkov,
\newblock JINR Report E2-89-655, Dubna, 1989.


\bibitem{parisi}
G. Parisi and N. Sourlas,
\newblock{\it Nucl. Phys.} {\bf B151} (1979) 421;
\\ I. S. Barker, C. B. Langensiepen and  G. Shaw,
\newblock{\it Nucl. Phys.} {\bf B186} (1981) 61.



%\cite{Krivokhizhin:1987rz}

\bibitem{Krivokhizhin:1987rz}
  V.~G.~Krivokhizhin, S.~P.~Kurlovich, V.~V.~Sanadze, I.~A.~Savin, A.~V.~Sidorov and N.~B.~Skachkov,
%  ``QCD ANALYSIS OF SINGLET STRUCTURE FUNCTIONS USING JACOBI POLYNOMIALS: THE
 % DESCRIPTION OF THE METHOD,''
  Z.\ Phys.\ C {\bf 36} (1987) 51.
  %%CITATION = ZEPYA,C36,51;%%

%\cite{Krivokhizhin:1990ct}

\bibitem{Krivokhizhin:1990ct}
  V.~G.~Krivokhizhin {\it et al.},
%  ``Next-to-leading order QCD analysis of structure functions with the help of
%  Jacobi polynomials,''
  Z.\ Phys.\ C {\bf 48}, 347 (1990).
  %%CITATION = ZEPYA,C48,347;%%

%\cite{Chyla:1986eb}

\bibitem{Chyla:1986eb}
  J.~Chyla and J.~Rames,
%  ``On Methods Of Analyzing Scaling Violation In Deep Inelastic Scattering,''
  Z.\ Phys.\  C {\bf 31} (1986) 151.
  %%CITATION = ZEPYA,C31,151;%%



%\cite{Barker:1980wu}

\bibitem{Barker:1980wu}
  I.~S.~Barker, C.~S.~Langensiepen and G.~Shaw,
 % ``General Parametrization Of Scale Breaking,''
  Nucl.\ Phys.\  B {\bf 186} (1981) 61.
  %%CITATION = NUPHA,B186,61;%%


%\cite{Kataev:1997nc}

\bibitem{Kataev:1997nc}
  A.~L.~Kataev, A.~V.~Kotikov, G.~Parente and A.~V.~Sidorov,
%   ``Next-to-next-to-leading order QCD analysis of the revised CCFR data for
 % xF3 structure function,''
  Phys.\ Lett.\ B {\bf 417}, (1998) 374
  [arXiv:hep-ph/9706534].
  %%CITATION = HEP-PH 9706534;%%

%\cite{Kataev:1998ce}

\bibitem{Kataev:1998ce}
  A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
 %  ``The {QCD} analysis of the CCFR data for xF3: Higher twists and
 % alpha(s)(M(Z)) extractions at the N$^2$LO and beyond,''
  arXiv:hep-ph/9809500.
  %%CITATION = HEP-PH 9809500;%%

%\cite{Alekhin:1998df}

%\bibitem{Alekhin:1998df}
%  S.~I.~Alekhin and A.~L.~Kataev,
%   ``The NLO DGLAP extraction of alpha(s) and higher twist terms from CCFR  xF3
 % and F2 structure functions data for nu N DIS,''
%  Phys.\ Lett.\ B {\bf 452}, (1999) 402
%  [arXiv:hep-ph/9812348].
  %%CITATION = HEP-PH 9812348;%%



%\cite{Kataev:1999bp}

\bibitem{Kataev:1999bp}
  A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
 %  ``Higher twists and alpha(s)(M(Z)) extractions from the N$^2$LO {QCD} analysis
 % of the CCFR data for the xF3 structure function,''
  Nucl.\ Phys.\ B {\bf 573}, (2000) 405
  [arXiv:hep-ph/9905310].
  %%CITATION = HEP-PH 9905310;%%




%\cite{Kataev:2001kk}
\bibitem{Kataev:2001kk}
  A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
%   ``Fixation of theoretical ambiguities in the improved fits to the xF3  CCFR
%  data at the next-to-next-to-leading order and beyond,''
  Phys.\ Part.\ Nucl.\  {\bf 34},  (2003) 20
   [arXiv:hep-ph/0106221]; \\A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
%``N$^3$LO fits to xF3 data: alpha(s) vs 1/Q**2 contributions,''
Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 116} (2003) 105
[arXiv:hep-ph/0211151].



%\cite{Khorramian:2006wg}
%\bibitem{Khorramian:2006wg}
%  A.~N.~Khorramian and S.~Atashbar Tehrani,
  %``The N$^2$LO nonsinglet QCD analysis of parton distributions based on
  %Bernstein polynomials,''
 % JHEP {\bf 0703}, 051 (2007)
 % [arXiv:hep-ph/0610136].
  %%CITATION = JHEPA,0703,051;%%

%\cite{Khorramian:2007zz}
\bibitem{Khorramian:2007zz}
  A.~N.~Khorramian, S.~Atashbar~Tehrani and M.~Ghominejad,
  %``The nonsinglet QCD analysis of parton distributions up to N$^2$LO based on
  %Jacobi polynomials,''
  Acta Phys.\ Polon.\  B {\bf 38}, 3551 (2007).
  %%CITATION = APPOA,B38,3551;%%

%\cite{Khorramian:2008zz}
\bibitem{Khorramian:2008zz}
  A.~N.~Khorramian and S.~A.~Tehrani,
  %``Nonsinglet QCD analysis in the NNLO approximation,''
  J.\ Phys.\ Conf.\ Ser.\  {\bf 110}, 022022 (2008).
  %%CITATION = 00462,110,022022;%%

%\cite{Khorramian:2008zza}
\bibitem{Khorramian:2008zza}
  A.~N.~Khorramian and S.~A.~Tehrani,
  %``nonsinglet Sector Of Proton And Deuteron Structure Functions,''
  AIP Conf.\ Proc.\  {\bf 1006} (2008) 118.
  %%CITATION = APCPC,1006,118;%%



%\cite{AtashbarTehrani:2009zz}
\bibitem{AtashbarTehrani:2009zz}
  S.~Atashbar Tehrani and A.~N.~Khorramian,
  %``3-Loops QCD Analysis For nonsinglet Sector Of F2,''
  Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 186}, 58 (2009).
  %%CITATION = NUPHZ,186,58;%%



%\cite{Khorramian:2009zz1}
\bibitem{Khorramian:2009zz1}
A.~N.~Khorramian, S.~Atashbar Tehrani, H. Khanpour and S. Taheri
Monfared, Hyperfine Interactions  {\bf 194}, 337 (2009).

%\cite{Khorramian:2009zz2}
\bibitem{Khorramian:2009zz2}
A.~N.~Khorramian, S.~Atashbar Tehrani, M. Soleymaninia. and S.
Batebi, Hyperfine Interactions  {\bf 194}, 341 (2009).

%\cite{Khorramian:2009zz3}
\bibitem{Khorramian:2009zz3}
S.~Atashbar Tehrani and A.~N.~Khorramian, Hyperfine Interactions
{\bf 194}, 331 (2009).

%\cite{Khorramian:2009zz4}
\bibitem{Khorramian:2009zz4}
S.~Atashbar Tehrani and A.~N.~Khorramian, Applied Mathematics \&
Information Sciences (2009),~367-373.


%\cite{Leader:1997kw}
\bibitem{Leader:1997kw}
  E.~Leader, A.~V.~Sidorov and D.~B.~Stamenov,
  %``NLO QCD analysis of polarized deep inelastic scattering,''
  Int.\ J.\ Mod.\ Phys.\  A {\bf 13}, 5573 (1998)
  [arXiv:hep-ph/9708335].
  %%CITATION = IMPAE,A13,5573;%%

%\cite{Atashbar Tehrani:2007be}
\bibitem{Atashbar Tehrani:2007be}
  S.~Atashbar Tehrani and A.~N.~Khorramian,
  %``The Jacobi Polynomials QCD analysis for the polarized structure function,''
  JHEP {\bf 0707}, 048 (2007)
  [arXiv:0705.2647 [hep-ph]].
  %%CITATION = JHEPA,0707,048;%%

%\cite{Khorramian:2007gu}
\bibitem{Khorramian:2007gu}
  A.~N.~Khorramian and S.~Atashbar~Tehrani,
  %``Simple model for QCD analysis of the proton helicity structure,''
  arXiv:0712.2373 [hep-ph].
  %%CITATION = ARXIV:0712.2373;%%


%\cite{Khorramian:2007zza}
\bibitem{Khorramian:2007zza}
  A.~N.~Khorramian and S.~Atashbar Tehrani,
  %``The Jacobi polynomials QCD analysis for proton spin structure,''
  AIP Conf.\ Proc.\  {\bf 915}, 420 (2007).
  %%CITATION = APCPC,915,420;%%

%\cite{Mirjalili:2007ep}

\bibitem{Mirjalili:2007ep}
  A.~Mirjalili, A.~N.~Khorramian and S.~Atashbar-Tehrani,
  %``QCD analysis of flavor-asymmetry for polarized sea quark distributions,''
  Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 164}, 38 (2007).
  %%CITATION = NUPHZ,164,38;%%

%\cite{Mirjalili:2006hf}
\bibitem{Mirjalili:2006hf}
  A.~Mirjalili, S.~Atashbar Tehrani and A.~N.~Khorramian,
  %``The role of polarized valons in the flavor symmetry breaking of nucleon
  %sea,''
  Int.\ J.\ Mod.\ Phys.\  A {\bf 21}, 4599 (2006)
  [arXiv:hep-ph/0608224].
  %%CITATION = IMPAE,A21,4599;%%


%[2]
\bibitem{BCDMS}
%\cite{Benvenuti:1989fm}
%\bibit{\it et al.}em{Benvenuti:1989fm}
A.C.~Benvenuti {\it et al.} [BCDMS Collaboration],
%``A High Statistics Measurement Of The Deuteron Structure Functions F2 (X,
%Q**2) And R From Deep Inelastic Muon Scattering At High Q**2,''
Phys.\ Lett.\ B {\bf 237} (1990) 592;
%%CITATION = PHLTA,B237,592;%%
\\
A.C.~Benvenuti {\it et al.} [BCDMS Collaboration], Phys. Lett.
{\bf B223} (1989) 485; Phys. Lett. {\bf B237} (1990) 592.
\\
%\cite{Benvenuti:1989gs}
%\bibitem{Benvenuti:1989gs}
A.C.~Benvenuti {\it et al.} [BCDMS Collaboration],
%``A Comparison Of The Structure Functions F2 Of The Proton And The Neutron From
%Deep Inelastic Muon Scattering At High Q**2,''
Phys.\ Lett.\ B {\bf 237} (1990) 599.
%%CITATION = PHLTA,B237,599;%%
%----------------------------------------------------------------------
%
%[3]
\bibitem{SLAC}
%\cite{Whitlow:1991uw}
%%\bibitem{Whitlow:1991uw}
L.~W.~Whitlow, E.~M.~Riordan, S.~Dasu, S.~Rock and A.~Bodek, %
% ``Precise measurements of the proton and deuteron structure functions from a
%global analysis of the SLAC deep inelastic electron scattering
%cross-sections,''
Phys.\ Lett.\ B {\bf 282} (1992). %
%%CITATION = PHLTA,B282,475;%%
%----------------------------------------------------------------------
%
%[4]
\bibitem{NMC}
%\cite{Arneodo:1996qe}
%\bibitem{Arneodo:1996qe}
M.~Arneodo {\it et al.}  [New Muon Collaboration],
%``Measurement of the proton and deuteron structure functions, F2(p) and  F2(d),
%and of the ratio sigma(L)/sigma(T),''
Nucl.\ Phys.\ B {\bf 483} (1997) 3 [arXiv:hep-ph/9610231].
%%CITATION = HEP-PH 9610231;%%
%----------------------------------------------------------------------
%
%[5]
\bibitem{H1}
%\cite{Adloff:2000qk}
%\bibitem{Adloff:2000qk}
C.~Adloff {\it et al.}  [H1 Collaboration],
%``Deep-inelastic inclusive e p scattering at low x and a determination of
%alpha(s),''
Eur.\ Phys.\ J.\ C {\bf 21} (2001) 33
[arXiv:hep-ex/0012053];\\
%%CITATION = HEP-EX 0012053;%%
%\cite{Adloff:2003uh}
%\bibitem{Adloff:2003uh}
C.~Adloff {\it et al.}  [H1 Collaboration],
%``Measurement and QCD analysis of neutral and charged current cross  sections
%at HERA,''
Eur.\ Phys.\ J.\ C {\bf 30} (2003) 1 [arXiv:hep-ex/0304003].
%%CITATION = HEP-EX 0304003;%%
%----------------------------------------------------------------------
%
%[6]
\bibitem{ZEUS}
%\cite{Breitweg:1998dz}
%\bibitem{Breitweg:1998dz}
J.~Breitweg {\it et al.}  [ZEUS Collaboration],
%``ZEUS results on the measurement and phenomenology of F2 at low x and  low
%Q**2,''
Eur.\ Phys.\ J.\ C {\bf 7} (1999) 609
[arXiv:hep-ex/9809005];\\
%%CITATION = HEP-EX 9809005;%%
%\cite{Chekanov:2001qu}
%\bibitem{Chekanov:2001qu}
S.~Chekanov {\it et al.}  [ZEUS Collaboration],
%``Measurement of the neutral current cross section and F2 structure  function
%for deep inelastic e+ p scattering at HERA,''
Eur.\ Phys.\ J.\ C {\bf 21} (2001) 443 [arXiv:hep-ex/0105090].
%%CITATION = HEP-EX 0105090;%%
%----------------------------------------------------------------------
%

\bibitem{Pade1}
  M.A. Samuel, J. Ellis and M. Karliner,
  Phys.\ Rev.\ Lett.\ 74 (1995) 4380
\bibitem{Pade2}
  J. Ellis, E. Gardi, M. Karliner and M.A. Samuel,
  Phys.~Lett.\ B366 (1996) 268
\bibitem{Pade3}
  J. Ellis, E. Gardi, M. Karliner and M.A. Samuel,
  Phys.\ Rev.\ D54 (1996) 6986

\bibitem{Baker}
G.A. Baker, Jr. {\em Essentials of Pad\'e Approximants},
 Academic Press, 1975.

\bibitem{BenderOrszag}C.M. Bender and S.A. Orszag,
{\sl Advanced Mathematical Methods for Scientists and Engineers},
McGraw-Hill, 1978.


%[11]
\bibitem{MRST03}
%\cite{Martin:2003tt}
%\bibitem{Martin:2003tt}
A.~D.~Martin, R.~G.~Roberts, W.~J.~Stirling and R.~S.~Thorne,
%``MRST partons and uncertainties,''
arXiv:hep-ph/0307262.
%%CITATION = HEP-PH 0307262;%%
%----------------------------------------------------------------------
%
%[12]
%\cite{Alekhin:2002fv}
\bibitem{A02}
%\bibitem{Alekhin:2002fv}
S.~Alekhin,
%``Parton distributions from deep-inelastic scattering data,''
Phys.\ Rev.\ D {\bf 68} (2003) 014002 [arXiv:hep-ph/0211096].
%%CITATION = HEP-PH 0211096;%%

%----------------------------------------------------------------------
%
%[14]
\bibitem{SY}
%\cite{Santiago:2001mh}
%\bibitem{Santiago:2001mh}
J.~Santiago and F.~J.~Yndurain,
%``Improved calculation of F2 in electroproduction and xF3 in neutrino
%scattering to NNLO and determination of alpha(s),''
Nucl.\ Phys.\ B {\bf 611} (2001) 447 [arXiv:hep-ph/0102247].
%%CITATION = HEP-PH 0102247;%%
%----------------------------------------------------------------------
%
%[15]
\bibitem{KAT}
%\cite{Kataev:2002rj}
%\bibitem{Kataev:2002rj}
A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
%``N**3LO fits to xF3 data: alpha(s) vs 1/Q**2 contributions,''
Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 116} (2003) 105
[arXiv:hep-ph/0211151];\\
%%CITATION = HEP-PH 0211151;%%
%\newpage
%\cite{Kataev:2001kk}
%\bibitem{Kataev:2001kk}
A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
%``Fixation of theoretical ambiguities in the improved fits to the xF3  CCFR
%data at the next-to-next-to-leading order and beyond,''
Phys.\ Part.\ Nucl.\  {\bf 34} (2003) 20 [Fiz.\ Elem.\ Chast.\
Atom.\ Yadra {\bf 34} (2003) 43] [arXiv:hep-ph/0106221].
%%CITATION = HEP-PH 0106221;%%
%----------------------------------------------------------------------
%
%[16]
\bibitem{SYS}
%\cite{Santiago:1999pr}
%\bibitem{Santiago:1999pr}
J.~Santiago and F.~J.~Yndurain,
%``Calculation of electroproduction to NNLO and precision determination of
%alpha(s),''
Nucl.\ Phys.\ B {\bf 563} (1999) 45 [arXiv:hep-ph/9904344].
%%CITATION = HEP-PH 9904344;%%
%----------------------------------------------------------------------
%
%[17]
\bibitem{BBG04}
%\cite{Blumlein:2004ip}
%\bibitem{Blumlein:2004ip}
J.~Bl\"umlein, H.~B\"ottcher and A.~Guffanti,
%``Nonsinglet QCD analysis of the structure function F2 in 3-loops,''
Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 135} (2004) 152
[arXiv:hep-ph/0407089];
%\cite{Blumlein:2006ws}
%\bibitem{Blumlein:2006ws}
%J.~Blumlein, H.~Bottcher and A.~Guffanti,
%``NNLO analysis of unpolarized DIS structure functions,''
arXiv:hep-ph/0606309.
%%CITATION = HEP-PH 0606309;%%
%----------------------------------------------------------------------
%
%[18]
%\cite{Gluck:2006yz}
\bibitem{Gluck:2006yz}
M.~Gl\"uck, E.~Reya and C.~Schuck,
%``Nonsinglet QCD analysis of F2(x,Q**2) up to NNLO,''
arXiv:hep-ph/0604116.
%%CITATION = HEP-PH 0604116;%%


%\cite{Blumlein:2006be}
\bibitem{Blumlein:2006be}
  J.~Blumlein, H.~Bottcher and A.~Guffanti,
  %``Nonsinglet QCD analysis of deep inelastic world data at O(alpha(s)**3),''
  Nucl.\ Phys.\  B {\bf 774}, 182 (2007)
  [arXiv:hep-ph/0607200].
  %%CITATION = NUPHA,B774,182;%%

%\cite{ABKM:2009aa}
\bibitem{ABKM:2009aa}
  S. Alekhin, J.~Blumlein, S. Klein,  S. Moch
  %``...,''
  arXiv:0908.2766 [hep-ph] and references  therein.
  %%CITATION = ARXIV:0907.1751;%%

\bibitem{ref3} W.L.~van Neerven, A.~Vogt, {\it Nucl.~Phys.}
               {\bf B568}, 263 (2000)
\bibitem{ref4} W.L.~van Neerven, A.~Vogt, {\it Nucl.~Phys.}
               {\bf B588}, 345 (2000) and arXiv:hep-ph/0006154
               (corrected)
\bibitem{ref5} J.~Bl\"umlein, A.~Vogt,
               {\it Phys.~Rev.} {\bf D58}, 014020 (1998).
               %\cite{Gluck:2006pm}
\bibitem{Gluck:2006pm}
  M.~Gluck, C.~Pisano and E.~Reya,
  %``The curvature of F2(p)(x,Q**2) as a probe of the range of validity of
  %perturbative QCD evolutions in the small-x region,''
  Eur.\ Phys.\ J.\  C {\bf 50}, 29 (2007)
  [arXiv:hep-ph/0610060].
  %%CITATION = EPHJA,C50,29;%%

%\cite{Vermaseren:2005qc}
\bibitem{Vermaseren:2005qc}
J.~A.~M.~Vermaseren, A.~Vogt and S.~Moch,
%``The third-order QCD corrections to deep-inelastic scattering by photon
%exchange,''
Nucl.\ Phys.\ B {\bf 724} (2005) 3 [arXiv:hep-ph/0504242].
%%CITATION = HEP-PH 0504242;%%

%\cite{Barone:1992ej}
\bibitem{Barone:1992ej}
  V.~Barone, M.~Genovese, N.~N.~Nikolaev, E.~Predazzi and B.~G.~Zakharov,
  %``Structure functions of bound nucleons: From the EMC effect to nuclear
  %shadowing,''
  Z.\ Phys.\  C {\bf 58}, 541 (1993).
  %%CITATION = ZEPYA,C58,541;%%

  %\cite{Kwiecinski:1987tb}
  \bibitem{Kwiecinski:1987tb}
  J.~Kwiecinski and B.~M.~Badelek,
  %``ANALYSIS OF THE ELECTROPRODUCTION STRUCTURE FUNCTIONS IN THE LOW Q**2
  %REGION COMBINING THE VECTOR MESON DOMINANCE AND THE PARTON MODEL WITH
  %POSSIBLE SCALING VIOLATION,''
  Z.\ Phys.\  C {\bf 43}, 251 (1989).
  %%CITATION = ZEPYA,C43,251;%%



%\cite{Amaudruz:1991at}
\bibitem{Amaudruz:1991at}
  P.~Amaudruz {\it et al.}  [New Muon Collaboration],
  %``The Gottfried Sum From The Ratio F2(N) / F2(P),''
  Phys.\ Rev.\ Lett.\  {\bf 66}, 2712 (1991);
  %%CITATION = PRLTA,66,2712;%%
M.~Arneodo {\it et al.}  [New Muon Collaboration],
  %``A Reevaluation of the Gottfried sum,''
  Phys.\ Rev.\  D {\bf 50}, 1 (1994);
  %\cite{Arneodo:1996kd}
%\bibitem{Arneodo:1996kd}
  M.~Arneodo {\it et al.}  [New Muon Collaboration],
  %``Accurate measurement of F2(d)/F2(p) and R(d)-R(p),''
  Nucl.\ Phys.\  B {\bf 487}, 3 (1997)
  [arXiv:hep-ex/9611022].
  %%CITATION = NUPHA,B487,3;%%

\bibitem{ref19} A.D.~Martin  et al., {\it Eur.~Phys.~J.} {\bf C23}
               (2002) 73.
%\cite{Blumlein:2004pr}
\bibitem{Blumlein:2004pr} J.~Bl\"umlein, H,~B\"ottcher, and A.~Guffanti,
               {\it Nucl.~Phys.~B} (Proc.~Suppl.)
               {\bf 135} (2004) 152.

\bibitem{E866} R.S.~Towell et al., E866 Collab.,
               {\it Phys.~Rev.} {\bf D64} (2001) 052002.




\bibitem{FP}
%\cite{Furmanski:1981cw}
%\bibitem{Furmanski:1981cw}
W.~Furmanski and R.~Petronzio,
%``Lepton - Hadron Processes Beyond Leading Order In Quantum Chromodynamics,''
Z.\ Phys.\ C {\bf 11} (1982) 293.
%%CITATION = ZEPYA,C11,293;%%


\bibitem{NS2}
%\cite{vanNeerven:1991nn}
%\bibitem{vanNeerven:1991nn}
W.~L.~van Neerven and E.~B.~Zijlstra,
%``Order alpha-s**2 contributions to the deep inelastic Wilson coefficient,''
Phys.\ Lett.\ B {\bf 272} (1991) 127;\\
%%CITATION = PHLTA,B272,127;%%
%\cite{Zijlstra:1992qd}
%\bibitem{Zijlstra:1992qd}
E.~B.~Zijlstra and W.~L.~van Neerven,
%``Order alpha-s**2 QCD corrections to the deep inelastic proton structure
%functions F2 and F(L),''
Nucl.\ Phys.\ B {\bf 383} (1992) 525.
%%CITATION = NUPHA,B383,525;%%



%\cite{Tarasov:1980au}
\bibitem{Tarasov:1980au}
  O.~V.~Tarasov, A.~A.~Vladimirov and A.~Y.~Zharkov,
  %``The Gell-Mann-Low Function Of QCD In The Three Loop Approximation,''
  Phys.\ Lett.\  B {\bf 93}, 429 (1980).
  %%CITATION = PHLTA,B93,429;%%


%\cite{Larin:1993tp}
\bibitem{Larin:1993tp}
  S.~A.~Larin and J.~A.~M.~Vermaseren,
  %``The Three Loop QCD Beta Function And Anomalous Dimensions,''
  Phys.\ Lett.\  B {\bf 303}, 334 (1993)
  [arXiv:hep-ph/9302208].
  %%CITATION = PHLTA,B303,334;%%

%\cite{Vogt:2004ns}
\bibitem{Vogt:2004ns}
  A.~Vogt,
  %``Efficient evolution of unpolarized and polarized parton distributions  with
  %QCD-PEGASUS,''
  Comput.\ Phys.\ Commun.\  {\bf 170}, 65 (2005)
  [arXiv:hep-ph/0408244].
  %%CITATION = CPHCB,170,65;%%

%\cite{Chetyrkin:1997sg}
\bibitem{Chetyrkin:1997sg}
K.~G.~Chetyrkin, B.~A.~Kniehl and M.~Steinhauser,
%``Strong coupling constant with flavor thresholds at four loops in the
%MS-bar scheme,''
, Phys.\ Rev.\ Lett.\  {\bf 79} (1997) 2184
[arXiv:hep-ph/9706430]; S.~Bethke,
%``Determination of the QCD coupling alpha(s)''
, J.\ Phys.\ G {\bf 26} (2000) R27 [arXiv:hep-ex/0004021];
W.~Bernreuther and W.~Wetzel,
  %``Decoupling of heavy quarks in the minimal subtraction scheme,''
  Nucl.\ Phys.\  B {\bf 197}, 228 (1982)
  [Erratum-ibid.\  B {\bf 513}, 758 (1998)].


%[33]
\bibitem{BAR}
%\cite{bibitem{BAR}
%\bibitem{Bardeen:1978yd}
W.~A.~Bardeen, A.~J.~Buras, D.~W.~Duke and T.~Muta,
%``Deep Inelastic Scattering Beyond The Leading Order In Asymptotically Free
%Gauge Theories,''
Phys.\ Rev.\ D {\bf 18} (1978) 3998.
%%CITATION = PHRVA,D18,3998;%%
%----------------------------------------------------------------------

%[34]
\bibitem{R1998}
%\cite{Abe:1998ym}
%\bibitem{Abe:1998ym}
K.~Abe {\it et al.}  [E143 Collaboration],
%``Measurements of R = sigma(L)/sigma(T) for 0.03 < x < 0.1 and fit to  world
%data,''
Phys.\ Lett.\ B {\bf 452} (1999) 194 [arXiv:hep-ex/9808028].
%%CITATION = HEP-EX 9808028;%%


%\cite{Parisi:1978jv}
\bibitem{Parisi:1978jv}
  G.~Parisi and N.~Sourlas,
  %``A Simple Parametrization Of The Q**2 Dependence Of The Quark Distributions
  %In QCD,''
  Nucl.\ Phys.\  B {\bf 151}, 421 (1979).
  %%CITATION = NUPHA,B151,421;%%






\bibitem{Stump:2001gu}
  D.~Stump {\it et al.},
  %``Uncertainties of predictions from parton distribution functions. I:  The
  %Lagrange multiplier method,''
  Phys.\ Rev.\  D {\bf 65}, 014012 (2002)
  [arXiv:hep-ph/0101051].
  %%CITATION = PHRVA,D65,014012;%%


\bibitem{MINUIT}
F.~James, CERN Program Library, Long Writeup D506 (MINUIT).


%\cite{Georgi:1976ve}
\bibitem{Georgi:1976ve}
H.~Georgi and H.~D.~Politzer,
%``Freedom At Moderate Energies: Masses In Color Dynamics,''
Phys.\ Rev.\ D {\bf 14} (1976) 1829.
%%CITATION = PHRVA,D14,1829;%%





%\cite{Alekhin:2005gq}
%\bibitem{Alekhin:2005gq}
%  S.~Alekhin,
  %``Parton distribution functions from the precise N$^2$LO QCD fit,''
%  JETP Lett.\  {\bf 82}, 628 (2005)
%  [Pisma Zh.\ Eksp.\ Teor.\ Fiz.\  {\bf 82}, 710 (2005)]
%  [arXiv:hep-ph/0508248].
  %%CITATION = ZFPRA,82,710;%%


%  \bibitem{MRST04}
%\cite{Martin:2004ir}
%\bibitem{Martin:2004ir}
%A.~D.~Martin, R.~G.~Roberts, W.~J.~Stirling and R.~S.~Thorne,
%``Physical gluons and high-E(T) jets,''
%Phys.\ Lett.\ B {\bf 604} (2004) 61 [arXiv:hep-ph/0410230].
%%CITATION = HEP-PH 0410230;%%



%\cite{Yang:1999xg}
\bibitem{Yang:1999xg}
  U.~K.~Yang and A.~Bodek,
  %``Studies of higher twist and higher order effects in NLO and N$^2$LO QCD
  %analysis of lepton nucleon scattering data on F2 and R =
  %sigma(L)/sigma(T),''
  Eur.\ Phys.\ J.\  C {\bf 13}, 241 (2000)
  [arXiv:hep-ex/9908058].
  %%CITATION = EPHJA,C13,241;%%

%\cite{Blumlein:2008kz}
\bibitem{Blumlein:2008kz}
  J.~Blumlein and H.~Bottcher,
  %``Higher Twist Contributions to the Structure Functions F_2^p(x,Q^2) and
  %F_2^d(x,Q^2) at Large x and Higher Orders,''
  Phys.\ Lett.\  B {\bf 662}, 336 (2008)
  [arXiv:0802.0408 [hep-ph]].
  %%CITATION = PHLTA,B662,336;%%


%\cite{Program:summary}
\bibitem{Program:summary}
Program summary URL: http://particles.ipm.ir/links/QCD.htm.


%\cite{Gluck:2007ck}
\bibitem{Gluck:2007ck}
  M.~Gluck, P.~Jimenez-Delgado and E.~Reya,
  %``Dynamical parton distributions of the nucleon and very small-x physics,''
  Eur.\ Phys.\ J.\  C {\bf 53}, 355 (2008)
  [arXiv:0709.0614 [hep-ph]].
  %%CITATION = EPHJA,C53,355;%%

%\cite{Martin:2007bv}
\bibitem{Martin:2007bv}
  A.~D.~Martin, W.~J.~Stirling, R.~S.~Thorne and G.~Watt,
  %``Update of parton distributions at NNLO,''
  Phys.\ Lett.\  B {\bf 652}, 292 (2007)
  [arXiv:0706.0459 [hep-ph]].
  %%CITATION = PHLTA,B652,292;%%

%\cite{JimenezDelgado:2008hf}
\bibitem{JimenezDelgado:2008hf}
  P.~Jimenez-Delgado and E.~Reya,
  %``Dynamical NNLO parton distributions,''
  Phys.\ Rev.\  D {\bf 79}, 074023 (2009)
  [arXiv:0810.4274 [hep-ph]].
  %%CITATION = PHRVA,D79,074023;%%


%\cite{Alekhin:2006zm}
\bibitem{Alekhin:2006zm}
  S.~Alekhin, K.~Melnikov and F.~Petriello,
  %``Fixed target Drell-Yan data and N$^2$LO QCD fits of parton distribution
  %functions,''
  Phys.\ Rev.\  D {\bf 74}, 054033 (2006)
  [arXiv:hep-ph/0606237].
  %%CITATION = PHRVA,D74,054033;%%


%\cite{Brooks:2006wh}
\bibitem{Brooks:2006wh}
  P.~M.~Brooks and C.~J.~Maxwell,
  %``Improved analysis of moments of F3 in neutrino-nucleon scattering using
  %the Bernstein polynomial method,''
  Nucl.\ Phys.\  B {\bf 780}, 76 (2007)
  [arXiv:hep-ph/0610137].
  %%CITATION = NUPHA,B780,76;%%

%\cite{Bethke:2006ac}
\bibitem{Bethke:2006ac}
  S.~Bethke,
  %``Experimental tests of asymptotic freedom,''
  Prog.\ Part.\ Nucl.\ Phys.\  {\bf 58}, 351 (2007)
  [arXiv:hep-ex/0606035].
  %%CITATION = PPNPD,58,351;%%

%\cite{Bethke:2009jm}
\bibitem{Bethke:2009jm}
  S.~Bethke,
  %``The 2009 Wolrd Average of $\alpha_s (M_Z)$,''
  arXiv:0908.1135 [hep-ph].
  %%CITATION = ARXIV:0908.1135;%%




\end{thebibliography}






\end{document}

%
% ****** End of file apssamp.tex ******

0909.2665.extracted_bib
\begin{thebibliography}{99}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\cite{Altarelli:2009za}
\bibitem{Altarelli:2009za}
  G.~Altarelli,
  %``DIS 2009 Concluding Talk: Outlook and Perspective,''
  arXiv:0907.1751 [hep-ph].
  %%CITATION = ARXIV:0907.1751;%%


%\cite{CERN:2009}
\bibitem{CERN:2009}
http://indico.cern.ch/conferenceOtherViews.py?view=standard\&confId=43022.


%\cite{Martin:2009iq}
\bibitem{Martin:2009iq}
  A.~D.~Martin, W.~J.~Stirling, R.~S.~Thorne and G.~Watt,
  %``Parton distributions for the LHC,''
  arXiv:0901.0002 [hep-ph].
  %%CITATION = ARXIV:0901.0002;%%

%\cite{Thorne:2009ky}
%\bibitem{Thorne:2009ky}
%  R.~S.~Thorne, A.~D.~Martin, W.~J.~Stirling and G.~Watt,
  %``Status of MRST/MSTW PDF sets,''
%  arXiv:0907.2387 [hep-ph].
  %%CITATION = ARXIV:0907.2387;%%

%\cite{Thorne:2009me}
\bibitem{Thorne:2009me}
  R.~S.~Thorne,
  %``Summary of the HERA-LHC workshop,''
  Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 191}, 295 (2009)
  [arXiv:0902.4820 [hep-ph]].
  %%CITATION = NUPHZ,191,295;%%

%\cite{Nadolsky:2008zw}
\bibitem{Nadolsky:2008zw}
  P.~M.~Nadolsky {\it et al.},
  %``Implications of CTEQ global analysis for collider observables,''
  Phys.\ Rev.\  D {\bf 78}, 013004 (2008)
  [arXiv:0802.0007 [hep-ph]].
  %%CITATION = PHRVA,D78,013004;%%




%\cite{Khorramian:2008yh}
\bibitem{Khorramian:2008yh}
  A.~N.~Khorramian and S.~A.~Tehrani,
  %``NNLO QCD contributions to the flavor nonsinglet sector of ${\bf
  %{F_2(x,Q^2)}}$,''
  Phys.\ Rev.\  D {\bf 78}, 074019 (2008)
  [arXiv:0805.3063 [hep-ph]].
  %%CITATION = PHRVA,D78,074019;%%

%\cite{Altarelli:2008fd}
%\bibitem{Altarelli:2008fd}
%  G.~Altarelli,
  %``Status of the Standard Model at the LHC Start,''
 % arXiv:0804.4147 [hep-ph].
  %%CITATION = ARXIV:0804.4147;%%

%\cite{Altarelli:2008aj}
%\bibitem{Altarelli:2008aj}
%  G.~Altarelli, R.~D.~Ball and S.~Forte,
  %``Small x Resummation with Quarks: Deep-Inelastic Scattering,''
 % Nucl.\ Phys.\  B {\bf 799}, 199 (2008)
 % [arXiv:0802.0032 [hep-ph]].
  %%CITATION = NUPHA,B799,199;%%


%\cite{Altarelli:2008xp}
%\bibitem{Altarelli:2008xp}
%  G.~Altarelli, R.~D.~Ball and S.~Forte,
  %``Structure Function Resummation in small-x QCD,''
%  PoS {\bf RADCOR2007}, 028 (2007)
%  [arXiv:0802.0968 [hep-ph]].
  %%CITATION = POSCI,RADCOR2007,028;%%

%\cite{Barker}
\bibitem{Barker}
I. S. Barker, B. R. Martin and G. Shaw,
\newblock{\it Z. Phys.} {\bf C19} (1983)  147;
\\I. S. Barker and B. R. Martin,
\newblock{\it Z. Phys.} {\bf C24} (1984) 255;
\\S. P. Kurlovich, A. V. Sidorov and N. B. Skachkov,
\newblock JINR Report E2-89-655, Dubna, 1989.


\bibitem{parisi}
G. Parisi and N. Sourlas,
\newblock{\it Nucl. Phys.} {\bf B151} (1979) 421;
\\ I. S. Barker, C. B. Langensiepen and  G. Shaw,
\newblock{\it Nucl. Phys.} {\bf B186} (1981) 61.



%\cite{Krivokhizhin:1987rz}

\bibitem{Krivokhizhin:1987rz}
  V.~G.~Krivokhizhin, S.~P.~Kurlovich, V.~V.~Sanadze, I.~A.~Savin, A.~V.~Sidorov and N.~B.~Skachkov,
%  ``QCD ANALYSIS OF SINGLET STRUCTURE FUNCTIONS USING JACOBI POLYNOMIALS: THE
 % DESCRIPTION OF THE METHOD,''
  Z.\ Phys.\ C {\bf 36} (1987) 51.
  %%CITATION = ZEPYA,C36,51;%%

%\cite{Krivokhizhin:1990ct}

\bibitem{Krivokhizhin:1990ct}
  V.~G.~Krivokhizhin {\it et al.},
%  ``Next-to-leading order QCD analysis of structure functions with the help of
%  Jacobi polynomials,''
  Z.\ Phys.\ C {\bf 48}, 347 (1990).
  %%CITATION = ZEPYA,C48,347;%%

%\cite{Chyla:1986eb}

\bibitem{Chyla:1986eb}
  J.~Chyla and J.~Rames,
%  ``On Methods Of Analyzing Scaling Violation In Deep Inelastic Scattering,''
  Z.\ Phys.\  C {\bf 31} (1986) 151.
  %%CITATION = ZEPYA,C31,151;%%



%\cite{Barker:1980wu}

\bibitem{Barker:1980wu}
  I.~S.~Barker, C.~S.~Langensiepen and G.~Shaw,
 % ``General Parametrization Of Scale Breaking,''
  Nucl.\ Phys.\  B {\bf 186} (1981) 61.
  %%CITATION = NUPHA,B186,61;%%


%\cite{Kataev:1997nc}

\bibitem{Kataev:1997nc}
  A.~L.~Kataev, A.~V.~Kotikov, G.~Parente and A.~V.~Sidorov,
%   ``Next-to-next-to-leading order QCD analysis of the revised CCFR data for
 % xF3 structure function,''
  Phys.\ Lett.\ B {\bf 417}, (1998) 374
  [arXiv:hep-ph/9706534].
  %%CITATION = HEP-PH 9706534;%%

%\cite{Kataev:1998ce}

\bibitem{Kataev:1998ce}
  A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
 %  ``The {QCD} analysis of the CCFR data for xF3: Higher twists and
 % alpha(s)(M(Z)) extractions at the N$^2$LO and beyond,''
  arXiv:hep-ph/9809500.
  %%CITATION = HEP-PH 9809500;%%

%\cite{Alekhin:1998df}

%\bibitem{Alekhin:1998df}
%  S.~I.~Alekhin and A.~L.~Kataev,
%   ``The NLO DGLAP extraction of alpha(s) and higher twist terms from CCFR  xF3
 % and F2 structure functions data for nu N DIS,''
%  Phys.\ Lett.\ B {\bf 452}, (1999) 402
%  [arXiv:hep-ph/9812348].
  %%CITATION = HEP-PH 9812348;%%



%\cite{Kataev:1999bp}

\bibitem{Kataev:1999bp}
  A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
 %  ``Higher twists and alpha(s)(M(Z)) extractions from the N$^2$LO {QCD} analysis
 % of the CCFR data for the xF3 structure function,''
  Nucl.\ Phys.\ B {\bf 573}, (2000) 405
  [arXiv:hep-ph/9905310].
  %%CITATION = HEP-PH 9905310;%%




%\cite{Kataev:2001kk}
\bibitem{Kataev:2001kk}
  A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
%   ``Fixation of theoretical ambiguities in the improved fits to the xF3  CCFR
%  data at the next-to-next-to-leading order and beyond,''
  Phys.\ Part.\ Nucl.\  {\bf 34},  (2003) 20
   [arXiv:hep-ph/0106221]; \\A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
%``N$^3$LO fits to xF3 data: alpha(s) vs 1/Q**2 contributions,''
Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 116} (2003) 105
[arXiv:hep-ph/0211151].



%\cite{Khorramian:2006wg}
%\bibitem{Khorramian:2006wg}
%  A.~N.~Khorramian and S.~Atashbar Tehrani,
  %``The N$^2$LO nonsinglet QCD analysis of parton distributions based on
  %Bernstein polynomials,''
 % JHEP {\bf 0703}, 051 (2007)
 % [arXiv:hep-ph/0610136].
  %%CITATION = JHEPA,0703,051;%%

%\cite{Khorramian:2007zz}
\bibitem{Khorramian:2007zz}
  A.~N.~Khorramian, S.~Atashbar~Tehrani and M.~Ghominejad,
  %``The nonsinglet QCD analysis of parton distributions up to N$^2$LO based on
  %Jacobi polynomials,''
  Acta Phys.\ Polon.\  B {\bf 38}, 3551 (2007).
  %%CITATION = APPOA,B38,3551;%%

%\cite{Khorramian:2008zz}
\bibitem{Khorramian:2008zz}
  A.~N.~Khorramian and S.~A.~Tehrani,
  %``Nonsinglet QCD analysis in the NNLO approximation,''
  J.\ Phys.\ Conf.\ Ser.\  {\bf 110}, 022022 (2008).
  %%CITATION = 00462,110,022022;%%

%\cite{Khorramian:2008zza}
\bibitem{Khorramian:2008zza}
  A.~N.~Khorramian and S.~A.~Tehrani,
  %``nonsinglet Sector Of Proton And Deuteron Structure Functions,''
  AIP Conf.\ Proc.\  {\bf 1006} (2008) 118.
  %%CITATION = APCPC,1006,118;%%



%\cite{AtashbarTehrani:2009zz}
\bibitem{AtashbarTehrani:2009zz}
  S.~Atashbar Tehrani and A.~N.~Khorramian,
  %``3-Loops QCD Analysis For nonsinglet Sector Of F2,''
  Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 186}, 58 (2009).
  %%CITATION = NUPHZ,186,58;%%



%\cite{Khorramian:2009zz1}
\bibitem{Khorramian:2009zz1}
A.~N.~Khorramian, S.~Atashbar Tehrani, H. Khanpour and S. Taheri
Monfared, Hyperfine Interactions  {\bf 194}, 337 (2009).

%\cite{Khorramian:2009zz2}
\bibitem{Khorramian:2009zz2}
A.~N.~Khorramian, S.~Atashbar Tehrani, M. Soleymaninia. and S.
Batebi, Hyperfine Interactions  {\bf 194}, 341 (2009).

%\cite{Khorramian:2009zz3}
\bibitem{Khorramian:2009zz3}
S.~Atashbar Tehrani and A.~N.~Khorramian, Hyperfine Interactions
{\bf 194}, 331 (2009).

%\cite{Khorramian:2009zz4}
\bibitem{Khorramian:2009zz4}
S.~Atashbar Tehrani and A.~N.~Khorramian, Applied Mathematics \&
Information Sciences (2009),~367-373.


%\cite{Leader:1997kw}
\bibitem{Leader:1997kw}
  E.~Leader, A.~V.~Sidorov and D.~B.~Stamenov,
  %``NLO QCD analysis of polarized deep inelastic scattering,''
  Int.\ J.\ Mod.\ Phys.\  A {\bf 13}, 5573 (1998)
  [arXiv:hep-ph/9708335].
  %%CITATION = IMPAE,A13,5573;%%

%\cite{Atashbar Tehrani:2007be}
\bibitem{Atashbar Tehrani:2007be}
  S.~Atashbar Tehrani and A.~N.~Khorramian,
  %``The Jacobi Polynomials QCD analysis for the polarized structure function,''
  JHEP {\bf 0707}, 048 (2007)
  [arXiv:0705.2647 [hep-ph]].
  %%CITATION = JHEPA,0707,048;%%

%\cite{Khorramian:2007gu}
\bibitem{Khorramian:2007gu}
  A.~N.~Khorramian and S.~Atashbar~Tehrani,
  %``Simple model for QCD analysis of the proton helicity structure,''
  arXiv:0712.2373 [hep-ph].
  %%CITATION = ARXIV:0712.2373;%%


%\cite{Khorramian:2007zza}
\bibitem{Khorramian:2007zza}
  A.~N.~Khorramian and S.~Atashbar Tehrani,
  %``The Jacobi polynomials QCD analysis for proton spin structure,''
  AIP Conf.\ Proc.\  {\bf 915}, 420 (2007).
  %%CITATION = APCPC,915,420;%%

%\cite{Mirjalili:2007ep}

\bibitem{Mirjalili:2007ep}
  A.~Mirjalili, A.~N.~Khorramian and S.~Atashbar-Tehrani,
  %``QCD analysis of flavor-asymmetry for polarized sea quark distributions,''
  Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 164}, 38 (2007).
  %%CITATION = NUPHZ,164,38;%%

%\cite{Mirjalili:2006hf}
\bibitem{Mirjalili:2006hf}
  A.~Mirjalili, S.~Atashbar Tehrani and A.~N.~Khorramian,
  %``The role of polarized valons in the flavor symmetry breaking of nucleon
  %sea,''
  Int.\ J.\ Mod.\ Phys.\  A {\bf 21}, 4599 (2006)
  [arXiv:hep-ph/0608224].
  %%CITATION = IMPAE,A21,4599;%%


%[2]
\bibitem{BCDMS}
%\cite{Benvenuti:1989fm}
%\bibit{\it et al.}em{Benvenuti:1989fm}
A.C.~Benvenuti {\it et al.} [BCDMS Collaboration],
%``A High Statistics Measurement Of The Deuteron Structure Functions F2 (X,
%Q**2) And R From Deep Inelastic Muon Scattering At High Q**2,''
Phys.\ Lett.\ B {\bf 237} (1990) 592;
%%CITATION = PHLTA,B237,592;%%
\\
A.C.~Benvenuti {\it et al.} [BCDMS Collaboration], Phys. Lett.
{\bf B223} (1989) 485; Phys. Lett. {\bf B237} (1990) 592.
\\
%\cite{Benvenuti:1989gs}
%\bibitem{Benvenuti:1989gs}
A.C.~Benvenuti {\it et al.} [BCDMS Collaboration],
%``A Comparison Of The Structure Functions F2 Of The Proton And The Neutron From
%Deep Inelastic Muon Scattering At High Q**2,''
Phys.\ Lett.\ B {\bf 237} (1990) 599.
%%CITATION = PHLTA,B237,599;%%
%----------------------------------------------------------------------
%
%[3]
\bibitem{SLAC}
%\cite{Whitlow:1991uw}
%%\bibitem{Whitlow:1991uw}
L.~W.~Whitlow, E.~M.~Riordan, S.~Dasu, S.~Rock and A.~Bodek, %
% ``Precise measurements of the proton and deuteron structure functions from a
%global analysis of the SLAC deep inelastic electron scattering
%cross-sections,''
Phys.\ Lett.\ B {\bf 282} (1992). %
%%CITATION = PHLTA,B282,475;%%
%----------------------------------------------------------------------
%
%[4]
\bibitem{NMC}
%\cite{Arneodo:1996qe}
%\bibitem{Arneodo:1996qe}
M.~Arneodo {\it et al.}  [New Muon Collaboration],
%``Measurement of the proton and deuteron structure functions, F2(p) and  F2(d),
%and of the ratio sigma(L)/sigma(T),''
Nucl.\ Phys.\ B {\bf 483} (1997) 3 [arXiv:hep-ph/9610231].
%%CITATION = HEP-PH 9610231;%%
%----------------------------------------------------------------------
%
%[5]
\bibitem{H1}
%\cite{Adloff:2000qk}
%\bibitem{Adloff:2000qk}
C.~Adloff {\it et al.}  [H1 Collaboration],
%``Deep-inelastic inclusive e p scattering at low x and a determination of
%alpha(s),''
Eur.\ Phys.\ J.\ C {\bf 21} (2001) 33
[arXiv:hep-ex/0012053];\\
%%CITATION = HEP-EX 0012053;%%
%\cite{Adloff:2003uh}
%\bibitem{Adloff:2003uh}
C.~Adloff {\it et al.}  [H1 Collaboration],
%``Measurement and QCD analysis of neutral and charged current cross  sections
%at HERA,''
Eur.\ Phys.\ J.\ C {\bf 30} (2003) 1 [arXiv:hep-ex/0304003].
%%CITATION = HEP-EX 0304003;%%
%----------------------------------------------------------------------
%
%[6]
\bibitem{ZEUS}
%\cite{Breitweg:1998dz}
%\bibitem{Breitweg:1998dz}
J.~Breitweg {\it et al.}  [ZEUS Collaboration],
%``ZEUS results on the measurement and phenomenology of F2 at low x and  low
%Q**2,''
Eur.\ Phys.\ J.\ C {\bf 7} (1999) 609
[arXiv:hep-ex/9809005];\\
%%CITATION = HEP-EX 9809005;%%
%\cite{Chekanov:2001qu}
%\bibitem{Chekanov:2001qu}
S.~Chekanov {\it et al.}  [ZEUS Collaboration],
%``Measurement of the neutral current cross section and F2 structure  function
%for deep inelastic e+ p scattering at HERA,''
Eur.\ Phys.\ J.\ C {\bf 21} (2001) 443 [arXiv:hep-ex/0105090].
%%CITATION = HEP-EX 0105090;%%
%----------------------------------------------------------------------
%

\bibitem{Pade1}
  M.A. Samuel, J. Ellis and M. Karliner,
  Phys.\ Rev.\ Lett.\ 74 (1995) 4380
\bibitem{Pade2}
  J. Ellis, E. Gardi, M. Karliner and M.A. Samuel,
  Phys.~Lett.\ B366 (1996) 268
\bibitem{Pade3}
  J. Ellis, E. Gardi, M. Karliner and M.A. Samuel,
  Phys.\ Rev.\ D54 (1996) 6986

\bibitem{Baker}
G.A. Baker, Jr. {\em Essentials of Pad\'e Approximants},
 Academic Press, 1975.

\bibitem{BenderOrszag}C.M. Bender and S.A. Orszag,
{\sl Advanced Mathematical Methods for Scientists and Engineers},
McGraw-Hill, 1978.


%[11]
\bibitem{MRST03}
%\cite{Martin:2003tt}
%\bibitem{Martin:2003tt}
A.~D.~Martin, R.~G.~Roberts, W.~J.~Stirling and R.~S.~Thorne,
%``MRST partons and uncertainties,''
arXiv:hep-ph/0307262.
%%CITATION = HEP-PH 0307262;%%
%----------------------------------------------------------------------
%
%[12]
%\cite{Alekhin:2002fv}
\bibitem{A02}
%\bibitem{Alekhin:2002fv}
S.~Alekhin,
%``Parton distributions from deep-inelastic scattering data,''
Phys.\ Rev.\ D {\bf 68} (2003) 014002 [arXiv:hep-ph/0211096].
%%CITATION = HEP-PH 0211096;%%

%----------------------------------------------------------------------
%
%[14]
\bibitem{SY}
%\cite{Santiago:2001mh}
%\bibitem{Santiago:2001mh}
J.~Santiago and F.~J.~Yndurain,
%``Improved calculation of F2 in electroproduction and xF3 in neutrino
%scattering to NNLO and determination of alpha(s),''
Nucl.\ Phys.\ B {\bf 611} (2001) 447 [arXiv:hep-ph/0102247].
%%CITATION = HEP-PH 0102247;%%
%----------------------------------------------------------------------
%
%[15]
\bibitem{KAT}
%\cite{Kataev:2002rj}
%\bibitem{Kataev:2002rj}
A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
%``N**3LO fits to xF3 data: alpha(s) vs 1/Q**2 contributions,''
Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 116} (2003) 105
[arXiv:hep-ph/0211151];\\
%%CITATION = HEP-PH 0211151;%%
%\newpage
%\cite{Kataev:2001kk}
%\bibitem{Kataev:2001kk}
A.~L.~Kataev, G.~Parente and A.~V.~Sidorov,
%``Fixation of theoretical ambiguities in the improved fits to the xF3  CCFR
%data at the next-to-next-to-leading order and beyond,''
Phys.\ Part.\ Nucl.\  {\bf 34} (2003) 20 [Fiz.\ Elem.\ Chast.\
Atom.\ Yadra {\bf 34} (2003) 43] [arXiv:hep-ph/0106221].
%%CITATION = HEP-PH 0106221;%%
%----------------------------------------------------------------------
%
%[16]
\bibitem{SYS}
%\cite{Santiago:1999pr}
%\bibitem{Santiago:1999pr}
J.~Santiago and F.~J.~Yndurain,
%``Calculation of electroproduction to NNLO and precision determination of
%alpha(s),''
Nucl.\ Phys.\ B {\bf 563} (1999) 45 [arXiv:hep-ph/9904344].
%%CITATION = HEP-PH 9904344;%%
%----------------------------------------------------------------------
%
%[17]
\bibitem{BBG04}
%\cite{Blumlein:2004ip}
%\bibitem{Blumlein:2004ip}
J.~Bl\"umlein, H.~B\"ottcher and A.~Guffanti,
%``Nonsinglet QCD analysis of the structure function F2 in 3-loops,''
Nucl.\ Phys.\ Proc.\ Suppl.\  {\bf 135} (2004) 152
[arXiv:hep-ph/0407089];
%\cite{Blumlein:2006ws}
%\bibitem{Blumlein:2006ws}
%J.~Blumlein, H.~Bottcher and A.~Guffanti,
%``NNLO analysis of unpolarized DIS structure functions,''
arXiv:hep-ph/0606309.
%%CITATION = HEP-PH 0606309;%%
%----------------------------------------------------------------------
%
%[18]
%\cite{Gluck:2006yz}
\bibitem{Gluck:2006yz}
M.~Gl\"uck, E.~Reya and C.~Schuck,
%``Nonsinglet QCD analysis of F2(x,Q**2) up to NNLO,''
arXiv:hep-ph/0604116.
%%CITATION = HEP-PH 0604116;%%


%\cite{Blumlein:2006be}
\bibitem{Blumlein:2006be}
  J.~Blumlein, H.~Bottcher and A.~Guffanti,
  %``Nonsinglet QCD analysis of deep inelastic world data at O(alpha(s)**3),''
  Nucl.\ Phys.\  B {\bf 774}, 182 (2007)
  [arXiv:hep-ph/0607200].
  %%CITATION = NUPHA,B774,182;%%

%\cite{ABKM:2009aa}
\bibitem{ABKM:2009aa}
  S. Alekhin, J.~Blumlein, S. Klein,  S. Moch
  %``...,''
  arXiv:0908.2766 [hep-ph] and references  therein.
  %%CITATION = ARXIV:0907.1751;%%

\bibitem{ref3} W.L.~van Neerven, A.~Vogt, {\it Nucl.~Phys.}
               {\bf B568}, 263 (2000)
\bibitem{ref4} W.L.~van Neerven, A.~Vogt, {\it Nucl.~Phys.}
               {\bf B588}, 345 (2000) and arXiv:hep-ph/0006154
               (corrected)
\bibitem{ref5} J.~Bl\"umlein, A.~Vogt,
               {\it Phys.~Rev.} {\bf D58}, 014020 (1998).
               %\cite{Gluck:2006pm}
\bibitem{Gluck:2006pm}
  M.~Gluck, C.~Pisano and E.~Reya,
  %``The curvature of F2(p)(x,Q**2) as a probe of the range of validity of
  %perturbative QCD evolutions in the small-x region,''
  Eur.\ Phys.\ J.\  C {\bf 50}, 29 (2007)
  [arXiv:hep-ph/0610060].
  %%CITATION = EPHJA,C50,29;%%

%\cite{Vermaseren:2005qc}
\bibitem{Vermaseren:2005qc}
J.~A.~M.~Vermaseren, A.~Vogt and S.~Moch,
%``The third-order QCD corrections to deep-inelastic scattering by photon
%exchange,''
Nucl.\ Phys.\ B {\bf 724} (2005) 3 [arXiv:hep-ph/0504242].
%%CITATION = HEP-PH 0504242;%%

%\cite{Barone:1992ej}
\bibitem{Barone:1992ej}
  V.~Barone, M.~Genovese, N.~N.~Nikolaev, E.~Predazzi and B.~G.~Zakharov,
  %``Structure functions of bound nucleons: From the EMC effect to nuclear
  %shadowing,''
  Z.\ Phys.\  C {\bf 58}, 541 (1993).
  %%CITATION = ZEPYA,C58,541;%%

  %\cite{Kwiecinski:1987tb}
  \bibitem{Kwiecinski:1987tb}
  J.~Kwiecinski and B.~M.~Badelek,
  %``ANALYSIS OF THE ELECTROPRODUCTION STRUCTURE FUNCTIONS IN THE LOW Q**2
  %REGION COMBINING THE VECTOR MESON DOMINANCE AND THE PARTON MODEL WITH
  %POSSIBLE SCALING VIOLATION,''
  Z.\ Phys.\  C {\bf 43}, 251 (1989).
  %%CITATION = ZEPYA,C43,251;%%



%\cite{Amaudruz:1991at}
\bibitem{Amaudruz:1991at}
  P.~Amaudruz {\it et al.}  [New Muon Collaboration],
  %``The Gottfried Sum From The Ratio F2(N) / F2(P),''
  Phys.\ Rev.\ Lett.\  {\bf 66}, 2712 (1991);
  %%CITATION = PRLTA,66,2712;%%
M.~Arneodo {\it et al.}  [New Muon Collaboration],
  %``A Reevaluation of the Gottfried sum,''
  Phys.\ Rev.\  D {\bf 50}, 1 (1994);
  %\cite{Arneodo:1996kd}
%\bibitem{Arneodo:1996kd}
  M.~Arneodo {\it et al.}  [New Muon Collaboration],
  %``Accurate measurement of F2(d)/F2(p) and R(d)-R(p),''
  Nucl.\ Phys.\  B {\bf 487}, 3 (1997)
  [arXiv:hep-ex/9611022].
  %%CITATION = NUPHA,B487,3;%%

\bibitem{ref19} A.D.~Martin  et al., {\it Eur.~Phys.~J.} {\bf C23}
               (2002) 73.
%\cite{Blumlein:2004pr}
\bibitem{Blumlein:2004pr} J.~Bl\"umlein, H,~B\"ottcher, and A.~Guffanti,
               {\it Nucl.~Phys.~B} (Proc.~Suppl.)
               {\bf 135} (2004) 152.

\bibitem{E866} R.S.~Towell et al., E866 Collab.,
               {\it Phys.~Rev.} {\bf D64} (2001) 052002.




\bibitem{FP}
%\cite{Furmanski:1981cw}
%\bibitem{Furmanski:1981cw}
W.~Furmanski and R.~Petronzio,
%``Lepton - Hadron Processes Beyond Leading Order In Quantum Chromodynamics,''
Z.\ Phys.\ C {\bf 11} (1982) 293.
%%CITATION = ZEPYA,C11,293;%%


\bibitem{NS2}
%\cite{vanNeerven:1991nn}
%\bibitem{vanNeerven:1991nn}
W.~L.~van Neerven and E.~B.~Zijlstra,
%``Order alpha-s**2 contributions to the deep inelastic Wilson coefficient,''
Phys.\ Lett.\ B {\bf 272} (1991) 127;\\
%%CITATION = PHLTA,B272,127;%%
%\cite{Zijlstra:1992qd}
%\bibitem{Zijlstra:1992qd}
E.~B.~Zijlstra and W.~L.~van Neerven,
%``Order alpha-s**2 QCD corrections to the deep inelastic proton structure
%functions F2 and F(L),''
Nucl.\ Phys.\ B {\bf 383} (1992) 525.
%%CITATION = NUPHA,B383,525;%%



%\cite{Tarasov:1980au}
\bibitem{Tarasov:1980au}
  O.~V.~Tarasov, A.~A.~Vladimirov and A.~Y.~Zharkov,
  %``The Gell-Mann-Low Function Of QCD In The Three Loop Approximation,''
  Phys.\ Lett.\  B {\bf 93}, 429 (1980).
  %%CITATION = PHLTA,B93,429;%%


%\cite{Larin:1993tp}
\bibitem{Larin:1993tp}
  S.~A.~Larin and J.~A.~M.~Vermaseren,
  %``The Three Loop QCD Beta Function And Anomalous Dimensions,''
  Phys.\ Lett.\  B {\bf 303}, 334 (1993)
  [arXiv:hep-ph/9302208].
  %%CITATION = PHLTA,B303,334;%%

%\cite{Vogt:2004ns}
\bibitem{Vogt:2004ns}
  A.~Vogt,
  %``Efficient evolution of unpolarized and polarized parton distributions  with
  %QCD-PEGASUS,''
  Comput.\ Phys.\ Commun.\  {\bf 170}, 65 (2005)
  [arXiv:hep-ph/0408244].
  %%CITATION = CPHCB,170,65;%%

%\cite{Chetyrkin:1997sg}
\bibitem{Chetyrkin:1997sg}
K.~G.~Chetyrkin, B.~A.~Kniehl and M.~Steinhauser,
%``Strong coupling constant with flavor thresholds at four loops in the
%MS-bar scheme,''
, Phys.\ Rev.\ Lett.\  {\bf 79} (1997) 2184
[arXiv:hep-ph/9706430]; S.~Bethke,
%``Determination of the QCD coupling alpha(s)''
, J.\ Phys.\ G {\bf 26} (2000) R27 [arXiv:hep-ex/0004021];
W.~Bernreuther and W.~Wetzel,
  %``Decoupling of heavy quarks in the minimal subtraction scheme,''
  Nucl.\ Phys.\  B {\bf 197}, 228 (1982)
  [Erratum-ibid.\  B {\bf 513}, 758 (1998)].


%[33]
\bibitem{BAR}
%\cite{bibitem{BAR}
%\bibitem{Bardeen:1978yd}
W.~A.~Bardeen, A.~J.~Buras, D.~W.~Duke and T.~Muta,
%``Deep Inelastic Scattering Beyond The Leading Order In Asymptotically Free
%Gauge Theories,''
Phys.\ Rev.\ D {\bf 18} (1978) 3998.
%%CITATION = PHRVA,D18,3998;%%
%----------------------------------------------------------------------

%[34]
\bibitem{R1998}
%\cite{Abe:1998ym}
%\bibitem{Abe:1998ym}
K.~Abe {\it et al.}  [E143 Collaboration],
%``Measurements of R = sigma(L)/sigma(T) for 0.03 < x < 0.1 and fit to  world
%data,''
Phys.\ Lett.\ B {\bf 452} (1999) 194 [arXiv:hep-ex/9808028].
%%CITATION = HEP-EX 9808028;%%


%\cite{Parisi:1978jv}
\bibitem{Parisi:1978jv}
  G.~Parisi and N.~Sourlas,
  %``A Simple Parametrization Of The Q**2 Dependence Of The Quark Distributions
  %In QCD,''
  Nucl.\ Phys.\  B {\bf 151}, 421 (1979).
  %%CITATION = NUPHA,B151,421;%%






\bibitem{Stump:2001gu}
  D.~Stump {\it et al.},
  %``Uncertainties of predictions from parton distribution functions. I:  The
  %Lagrange multiplier method,''
  Phys.\ Rev.\  D {\bf 65}, 014012 (2002)
  [arXiv:hep-ph/0101051].
  %%CITATION = PHRVA,D65,014012;%%


\bibitem{MINUIT}
F.~James, CERN Program Library, Long Writeup D506 (MINUIT).


%\cite{Georgi:1976ve}
\bibitem{Georgi:1976ve}
H.~Georgi and H.~D.~Politzer,
%``Freedom At Moderate Energies: Masses In Color Dynamics,''
Phys.\ Rev.\ D {\bf 14} (1976) 1829.
%%CITATION = PHRVA,D14,1829;%%





%\cite{Alekhin:2005gq}
%\bibitem{Alekhin:2005gq}
%  S.~Alekhin,
  %``Parton distribution functions from the precise N$^2$LO QCD fit,''
%  JETP Lett.\  {\bf 82}, 628 (2005)
%  [Pisma Zh.\ Eksp.\ Teor.\ Fiz.\  {\bf 82}, 710 (2005)]
%  [arXiv:hep-ph/0508248].
  %%CITATION = ZFPRA,82,710;%%


%  \bibitem{MRST04}
%\cite{Martin:2004ir}
%\bibitem{Martin:2004ir}
%A.~D.~Martin, R.~G.~Roberts, W.~J.~Stirling and R.~S.~Thorne,
%``Physical gluons and high-E(T) jets,''
%Phys.\ Lett.\ B {\bf 604} (2004) 61 [arXiv:hep-ph/0410230].
%%CITATION = HEP-PH 0410230;%%



%\cite{Yang:1999xg}
\bibitem{Yang:1999xg}
  U.~K.~Yang and A.~Bodek,
  %``Studies of higher twist and higher order effects in NLO and N$^2$LO QCD
  %analysis of lepton nucleon scattering data on F2 and R =
  %sigma(L)/sigma(T),''
  Eur.\ Phys.\ J.\  C {\bf 13}, 241 (2000)
  [arXiv:hep-ex/9908058].
  %%CITATION = EPHJA,C13,241;%%

%\cite{Blumlein:2008kz}
\bibitem{Blumlein:2008kz}
  J.~Blumlein and H.~Bottcher,
  %``Higher Twist Contributions to the Structure Functions F_2^p(x,Q^2) and
  %F_2^d(x,Q^2) at Large x and Higher Orders,''
  Phys.\ Lett.\  B {\bf 662}, 336 (2008)
  [arXiv:0802.0408 [hep-ph]].
  %%CITATION = PHLTA,B662,336;%%


%\cite{Program:summary}
\bibitem{Program:summary}
Program summary URL: http://particles.ipm.ir/links/QCD.htm.


%\cite{Gluck:2007ck}
\bibitem{Gluck:2007ck}
  M.~Gluck, P.~Jimenez-Delgado and E.~Reya,
  %``Dynamical parton distributions of the nucleon and very small-x physics,''
  Eur.\ Phys.\ J.\  C {\bf 53}, 355 (2008)
  [arXiv:0709.0614 [hep-ph]].
  %%CITATION = EPHJA,C53,355;%%

%\cite{Martin:2007bv}
\bibitem{Martin:2007bv}
  A.~D.~Martin, W.~J.~Stirling, R.~S.~Thorne and G.~Watt,
  %``Update of parton distributions at NNLO,''
  Phys.\ Lett.\  B {\bf 652}, 292 (2007)
  [arXiv:0706.0459 [hep-ph]].
  %%CITATION = PHLTA,B652,292;%%

%\cite{JimenezDelgado:2008hf}
\bibitem{JimenezDelgado:2008hf}
  P.~Jimenez-Delgado and E.~Reya,
  %``Dynamical NNLO parton distributions,''
  Phys.\ Rev.\  D {\bf 79}, 074023 (2009)
  [arXiv:0810.4274 [hep-ph]].
  %%CITATION = PHRVA,D79,074023;%%


%\cite{Alekhin:2006zm}
\bibitem{Alekhin:2006zm}
  S.~Alekhin, K.~Melnikov and F.~Petriello,
  %``Fixed target Drell-Yan data and N$^2$LO QCD fits of parton distribution
  %functions,''
  Phys.\ Rev.\  D {\bf 74}, 054033 (2006)
  [arXiv:hep-ph/0606237].
  %%CITATION = PHRVA,D74,054033;%%


%\cite{Brooks:2006wh}
\bibitem{Brooks:2006wh}
  P.~M.~Brooks and C.~J.~Maxwell,
  %``Improved analysis of moments of F3 in neutrino-nucleon scattering using
  %the Bernstein polynomial method,''
  Nucl.\ Phys.\  B {\bf 780}, 76 (2007)
  [arXiv:hep-ph/0610137].
  %%CITATION = NUPHA,B780,76;%%

%\cite{Bethke:2006ac}
\bibitem{Bethke:2006ac}
  S.~Bethke,
  %``Experimental tests of asymptotic freedom,''
  Prog.\ Part.\ Nucl.\ Phys.\  {\bf 58}, 351 (2007)
  [arXiv:hep-ex/0606035].
  %%CITATION = PPNPD,58,351;%%

%\cite{Bethke:2009jm}
\bibitem{Bethke:2009jm}
  S.~Bethke,
  %``The 2009 Wolrd Average of $\alpha_s (M_Z)$,''
  arXiv:0908.1135 [hep-ph].
  %%CITATION = ARXIV:0908.1135;%%




\end{thebibliography}