QCD analysis of non-singlet structure functions at NNLO accuracy, based on the Laplace transform

TL;DR

This paper employs Laplace transform techniques to analyze non-singlet quark distributions and nucleon structure functions at NNLO accuracy, incorporating target mass and higher twist corrections, and compares results with contemporary PDFs and experimental data.

Contribution

It introduces a novel application of Laplace transforms combined with Jacobi polynomial expansion for NNLO QCD analysis of non-singlet structure functions, including target mass and higher twist effects.

Findings

Results agree well with recent experimental data.

Comparison with modern PDF sets shows consistency.

Method provides a precise framework for structure function analysis.

Abstract

In this work, using the Laplace transformation technique we present our results for non-singlet quark distributions as well as nucleon structure function $F_2(x,Q^2)$ in unpolarized case at next-to-next-to-leading order (NNLO) QCD accuracy. We shall particularly compare our results for the sets of valence-quark parton distribution functions with the contemporary collaborations like CT14, CT18, MMHT14, MKAM16 and NNPDF. To construct the nucleon structure function we employ the expansion of Jacobi polynomials which is a suitable transform to convert the results of non-singlet structure function from the Laplace $s$-space to Bjorken $x$-space. We shall also consider the contributions of target mass correction as well as the higher twist effects at large-$x$ region for the proton and deuteron structure functions. Our results for the unpolarized quark distribution functions and nucleon structure functions are in good agreement with recent theoretical models and available experimental data.

Tables

Table 1. Table 1: Parameter values of NNLO non-singlet QCD fit at Q = 0 2 2 GeV 2 {}_{0}^{2}=2\,{\rm GeV}^{2} .

NNLO non-singlet QCD fit
$u_v$	${\alpha }_u$	0.70639 $\pm$ 0.0205
	${\beta }_u$	3.5318 $\pm$ 0.0183
	${\gamma }_u$	1.000
	${\eta }_u$	1.1400
$d_v$	${\alpha }_d$	0.71113 $\pm$ 0.0173
	${\beta }_d$	4.2114 $\pm$ 0.0955
	${\gamma }_d$	1.5999
	${\eta }_d$	4.2899
${\alpha }_s^{N_f=4}\left(Q_0^2\right)$		0.3632 $\pm$ 0.0149
${\chi }^2/n.d.f$		510.113/563 = 0.906

Table 2. Table 2: Parameter values of NLO non-singlet QCD fit at Q = 0 2 2 GeV 2 {}_{0}^{2}=2\,{\rm GeV}^{2}   Khanpour:2016uxh .

NLO non-singlet QCD fit
$u_v$	${\alpha }_u$	0.7108 $\pm$ 0.1295
	${\beta }_u$	3.3595 $\pm$ 0.027
	${\gamma }_u$	0.2979
	${\eta }_u$	1.3440
$d_v$	${\alpha }_d$	0.9467 $\pm$ 0.0261
	${\beta }_d$	2.8468 $\pm$ 0.3130
	${\gamma }_d$	1.1004
	${\eta }_d$	-1.1330
${\alpha }_s^{N_f=4}\left(Q_0^2\right)$		0.3521 $\pm$ 0.0139
${\chi }^2/n.d.f$		521.303/563 = 0.92

Table 3. Table 3: Data sets used in our QCD analysis for (a) F 2 p ​ ( x , Q 2 ) superscript subscript 𝐹 2 𝑝 𝑥 superscript 𝑄 2 F_{2}^{p}(x,Q^{2}) , (b) F 2 d ​ ( x , Q 2 ) superscript subscript 𝐹 2 𝑑 𝑥 superscript 𝑄 2 F_{2}^{d}(x,Q^{2}) , and (c) F 2 NS ​ ( x , Q 2 ) superscript subscript 𝐹 2 NS 𝑥 superscript 𝑄 2 F_{2}^{\rm NS}(x,Q^{2}) . The name of different data sets as well as the range of x 𝑥 x and Q 2 are given in three first columns. The normalization shifts are also listed in the last column.

Experiment	$x$	Q${}_{}{}^2\left({\mathrm{GeV}}^2\right)$	$F_2^p$	$F_2^p\mathrm{cuts}$	$𝒩$
BCDMS (100)	0.35–0.75	11.75–75.00	51	29	0.996805
BCDMS (120)	0.35–0.75	13.25–75.00	59	32	0.996805
BCDMS (200)	0.35–0.75	32.50–137.50	50	28	0.997833
BCDMS (280)	0.35–0.75	43.00–230.00	49	26	1.002131
NMC (comb)	0.35–0.50	7.00–65.00	15	14	1.000158
SLAC (comb)	0.30–0.62	7.30–21.39	57	57	1.000714
H1 (hQ2)	0.40–0.65	200–30000	26	26	1.001000
ZEUS (hQ2)	0.40–0.65	650–30000	15	15	0.999929
H1 (comb)	0.40–0.65	90–30000	145	21	0.999947
proton			467	248

Table 4. Table 4: The correlation matrix elements for the 4 + 1 free parameters at NNLO fit.

	${\alpha }_{u_v}$	${\beta }_{u_v}$	${\alpha }_{d_v}$	${\beta }_{d_v}$	${\Lambda }_{\overline{\mathrm{MS}}}^{\left(4\right)}$
${\alpha }_{u_v}$	1
${\beta }_{u_v}$	0.0957708	1
${\alpha }_{d_v}$	-0.0500717	-0.294017	1
${\beta }_{d_v}$	-0.0315322	-0.23711	0.910278	1
${\Lambda }_{\overline{\mathrm{MS}}}^{\left(2\right)}$	0.0445004	-0.585629	-0.0247452	-0.127968	1

Table 5. Table 5: The values of HT parameters in Eq. ( 32 ) extracted through the NNLO fit at Q = 0 2 2 GeV 2 {}_{0}^{2}=2\,{\rm GeV}^{2} .

NNLO
$h\left(x\right)$	$\alpha =1.07092$	$\beta =0.842737$	$\gamma =1.12905$

Equations

$$\displaystyle\frac{4\pi}{\alpha_{s}(Q^{2})}\frac{\partial q_{\text{NS}}}{\partial\ln Q^{2}}(x,Q^{2})=q_{\text{NS}}\otimes\Big{(}p_{NS}^{\text{LO}}+$$

$$\displaystyle\frac{\alpha_{s}(Q^{2})}{4\pi}p_{NS}^{\text{NLO}}+(\frac{\alpha_{s}(Q^{2})}{4\pi})^{2}p_{NS}^{\text{NNLO}}\Big{)}(x,Q^{2})\;.$$

$$p_{NS}({\alpha }_s(Q^2))$$

$$\displaystyle\frac{\partial\hat{F}_{\text{ NS}}}{\partial\tau}(\nu,\tau)=\int_{0}^{\nu}\hat{F}_{\text{ NS}}(w,\tau)e^{-(\nu-w)}\,dw\Big{(}p_{NS}^{\text{ LO}}(\nu-w)+$$

$$\displaystyle\frac{\alpha_{s}(\tau)}{4\pi}p_{NS}^{\text{ NLO}}(\nu-w)+(\frac{\alpha_{s}(\tau)}{4\pi})^{2}p_{NS}^{\text{ NNLO}}(\nu-w)\Big{)}\,.$$

$$\frac{\partial f_{\text{ NS}}}{\partial \tau }(s,\tau )$$

$$f_{\text{ NS}}(s,\tau )=e^{\tau {\Phi }_{\text{NS}}(s)}{f}_{\text{ NS}}^0(s).$$

$${\Phi }_{\text{NS}}(s)\equiv {\Phi }_{\text{ NS}}^{\text{LO}}(s)+\frac{{\tau }_2}{\tau }{\Phi }_{\text{NS}}^{\text{ NLO}}(s)+\frac{{\tau }_3}{\tau }{\Phi }_{\text{NS}}^{\text{NNLO}}(s).$$

$${\tau }_2\equiv \frac{1}{(4\pi {)}^2}{\int }_{Q_0^2}^{Q^2}{\alpha }_s^2(Q^{\prime 2})d\ln Q^{\prime 2},$$

$${\tau }_3\equiv \frac{1}{(4\pi {)}^3}{\int }_{Q_0^2}^{Q^2}{\alpha }_s^3(Q^{\prime 2})d\ln Q^{\prime 2},$$

$${\Phi }_f^{\text{LO}}=4-\frac{8}{3}\left(\frac{1}{s+1}+\frac{1}{s+2}+2\left({\gamma }_E+\psi (s+1)\right)\right),$$

$${\mathcal{M}}_2^p(s,\tau )$$

$${\mathcal{M}}_2^d(s,\tau )$$

$${\mathcal{M}}_2^{\mathrm{NS}}(s,\tau )$$

$$x(\bar{d}-\bar{u})(x,Q_0^2)=37.0x^{2.2}(1-x{)}^{19.2}(1+2.1\sqrt{x}).$$

$$x{u}_v={\mathcal{N}}_u{x}^{{\alpha }_{u_v}}(1-x{)}^{{\beta }_{u_v}}(1+{\gamma }_{u_v}{x}^{0.5}+{\eta }_{u_v}x),$$

$$x{d}_v={\mathcal{N}}_d{x}^{{\alpha }_{d_v}}(1-x{)}^{{\beta }_{d_v}}(1+{\gamma }_{d_v}{x}^{0.5}+{\eta }_{d_v}x),$$

$${\mathcal{N}}_u$$

$${\mathcal{N}}_d$$

$$u_v(s)=\mathcal{L}[e^{-\nu }{u}_v(e^{-\nu });s],$$

$$d_v(s)=\mathcal{L}[e^{-\nu }{d}_v(e^{-\nu });s],$$

$$u_v(s)=2(B[{\alpha }_{u_v}+s,{\beta }_{u_v}+1]+{\eta }_{u_v}B[{\alpha }_{u_v}+s+1,{\beta }_{u_v}+1]$$

$$+{\gamma }_{u_v}B[{\alpha }_{u_v}+s+0.5,{\beta }_{u_v}+1])/(B[{\alpha }_{u_v},{\beta }_{u_v}+1]+$$

$${\eta }_{u_v}B[{\alpha }_{u_v}+1,{\beta }_{u_v}+1]+{\gamma }_{u_v}B[{\alpha }_{u_v}+0.5,{\beta }_{u_v}+1]),$$

$$d_v(s)=(B[{\alpha }_{d_v}+s,{\beta }_{d_v}+1]+{\eta }_{d_v}B[{\alpha }_{d_v}+s+1,{\beta }_{d_v}+1]$$

$$+{\gamma }_{d_v}B[{\alpha }_{d_v}+s+0.5,{\beta }_{d_v}+1])/(B[{\alpha }_{d_v},{\beta }_{d_v}+1]+$$

$${\eta }_{d_v}B[{\alpha }_{d_v}+1,{\beta }_{d_v}+1]+{\gamma }_{d_v}B[{\alpha }_{d_v}+0.5,{\beta }_{d_v}+1]).$$

$$F_2^{p,d,NS}(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 )\mathcal{L}[F_2^{p,d,NS},j+1].$$

$${\Theta }_n^{\alpha ,\beta }(x)=\sum _{j=0}^n{c}_j^{(n)}(\alpha ,\beta )x^j,$$

$${\int }_0^{1}dx{x}^{\beta }(1-x{)}^{\alpha }{\Theta }_k^{\alpha ,\beta }(x){\Theta }_l^{\alpha ,\beta }(x)={\delta }_{k,l}.$$

$${\chi }_{\mathrm{global}}^2=\sum _n{w}_n{\chi }_n^2,$$

$${\chi }_n^2={\left(\frac{1-{\mathcal{N}}_n}{\Delta {\mathcal{N}}_n}\right)}^2+\sum _i{\left(\frac{{\mathcal{N}}_n{F}_{2,i}^{data}-F_{2,i}^{theor}}{{\mathcal{N}}_n\Delta F_{2,i}^{data}}\right)}^2.$$

$$\Delta {\chi }^2\equiv {\chi }_{\text{global}}^2-{\chi }_{\text{m}in}^2=\sum _{\text{ i, j=1}}^p(a_i-a_i^0)H_{\mathrm{ij}}(a_j-a_j^0).$$

Full text

I Introduction

As is well-known, our understanding of Quantum Chromodynamics (QCD) as well as the nucleon structure does profoundly depend on the deep inelastic scattering (DIS) processes. In this regard, the theory of QCD contains the necessary ingredients for the scale evolution of nucleon structure functions. On this base the proton, as a specific state of nucleons, can be described in perturbative QCD (pQCD) framework in terms of parton distribution functions (PDFs) which are known as the nonperturbative part of evolution process. The PDFs are typically related to the partons, i.e., gluon and quarks, and refer to the probability to find partons carrying away a specific fraction of proton’s momentum. In other words, the nonperturbative inputs are the PDFs at the initial energy scale which can be evolved to higher scales within the pQCD framework. Finally, they denote the probability density to find a parton which is carrying a specific fraction of proton’s momentum at a desired energy scale.

As usual, these nonperturbative PDFs can be determined by fitting the available experimental data including the DIS processes and the ones from hadron colliders Ball:2017nwa ; Bourrely:2015kla ; Harland-Lang:2014zoa ; Hou:2017khm ; Alekhin:2017kpj ; Khanpour:2017slc ; Goharipour:2018yov ; Soleymaninia:2018uiv ; Nejad:2015oca ; Nejad:2015far ; MoosaviNejad:2012ju ; Shoeibi:2017zha . To get numerical solutions for the evolved non-singlet PDFs, there are various methods such as the Brute-force Cabibbo:1978ez ; Miyama:1995bd ; Hirai:1997gb , the Laguerre transformation technique Toldra:2001yz ; Furmanski:1981ja ; Blumlein:1989pd ; Kumano:1992vd ; Kobayashi:1994hy ; GhasempourNesheli:2015tva and the Mellin-transform Gluck:1989ze ; Graudenz:1995sk ; Blumlein:1997em ; Blumlein:2000hw ; Stratmann:2001pb , and for the analytical solutions we can refer to the Jacobi polynomials model  Khorramian:2008yh ; Khorramian:2009xz as well as Laplace transformation technique Khanpour:2016uxh ; MoosaviNejad:2016ebo ; Boroun:2015hiy . The evolution of PDFs in $x$-space can be also done using the common programs like HOPPET hop , QCDNUM qcd , and APFEL apf .

In this paper, using the Laplace transformation technique we compute the non-singlet quark distributions as well as the unpolarized nucleon structure function at next-to-next-to-leading order (NNLO) of QCD analysis. To extract the NNLO PDFs we apply the deep-inelastic data for non-singlet QCD analysis employing the Laplace transformation at NNLO accuracy. Our results for the valence-quark PDFs are also compared with the famous Collaborations like CT14, CT18, MMHT14, MKAM16 and NNPDF. The expansion of Jacobi polynomials is also employed to construct the nucleon structure function. This is a convenient approach to convert the results of non-singlet structure function from the Laplace $s$-space to the Bjorken $x$-space. We will also consider the corrections due to the target mass as well as the higher twist effect for the proton and deuteron structure functions. These effects improve the quality of fit for structure functions at low energy scales.

The structure of our paper is organized as follows: In Sec. II, we present an analytical solution for NNLO non-singlet quark density based on the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations in Laplace $s$-space. Sec. III is devoted to the Jacobi polynomial technique which yields the non-singlet nucleon structure function in $x$-space. Our global analysis of valence-quark densities is presented in Sec. IV where we describe our procedure for the QCD fit of non-singlet $F_{2}$ structure function data. Corrections due to the target mass and the Higher twist effect are discussed in Sec. V. Our results and discussions are listed in Sec. VI.

As a supplement to this work, our analytical results for the Willson coefficient functions as well as the splitting functions in the Laplace $s$-space at NLO and NNLO accuracies are presented in appendix $A$.

II NNLO Non-singlet solution in Laplace space

In this work, our main interest is to investigate the proton spin-independent structure function, $F_{2}^{p}(x,Q^{2})$, in non-singlet case at the NNLO QCD analysis, specifically, in large values of $x$. Some analytical solutions of the DGLAP equations based on the Laplace transform technique have been recently presented, see for example Refs. Block:2010du ; Block:2011xb ; Block:2010fk ; Block:2009en ; Block:2010ti ; Zarei:2015jvh ; Taghavi-Shahri:2016ktz ; Boroun:2015cta ; Boroun:2014dka which contain remarkable success from phenomenological point of view. In Ref. Khanpour:2016uxh , using the mentioned technique authors have presented their spin-independent analysis of the structure function $F_{2}^{p}(x,Q^{2})$ at NLO accuracy. In Refs. MoosaviNejad:2016ebo and Sheibani:2018gxt , the application of Laplace transform technique to the analysis of charged-current structure functions $xF_{3}(x,Q^{2})$ as well as the EMC effects are studied. The spin-dependent structure functions $xg_{1}^{p}(x,Q^{2})$ are analyzed via the same technique in Refs. AtashbarTehrani:2013qea ; Salajegheh:2018hfs at NLO and NNLO accuracies.

For the non-singlet sector of quark densities, presented by $q_{\rm NS}(x,Q^{2})$, one can write the following DGLAP evolution equations at the NNLO expansion

[TABLE]

Here, the symbol $\otimes$ denotes the convolution integral and $\alpha_{s}(Q^{2})$ stands for the renormalized strong coupling constant. Moreover, the non-singlet Altarelli-Parisi splitting kernels up to three-loops corrections are presented by $p_{NS}^{LO}(\alpha_{s}(Q^{2}))$, $p_{NS}^{NLO}(\alpha_{s}(Q^{2}))$ and $p_{NS}^{NNLO}(\alpha_{s}(Q^{2}))$, respectively. Considering these kernels, the required splitting function has the following expansion

[TABLE]

A brief description to extract the analytical solution of the valence quark distribution function through the DGLAP evolution equations based on the Laplace transform technique is now at hand. Considering the explicit form of convolution integral in Eq. (II), in addition to the $x$-variable an extra $z$-variable would be also appeared. Taking the variable changes as $\nu\equiv\ln(1/x)$ and $w\equiv\ln(1/z)$, the evolution equation (II) is expressed in terms of the convolution integrals as

[TABLE]

As is seen, the $q_{\text{NS}}(x,Q^{2})$-function in Eq. (II) is now presented by $\hat{F}_{\text{NS}}$, including new variables $\nu$ and $\tau$. The $\tau$-variable is utilizing to present the $Q^{2}$-dependence of Eq. (II) and is entirely given by $\tau(Q^{2},Q_{0}^{2})=(1/4\pi)\int_{Q_{0}^{2}}^{Q^{2}}\alpha_{s}(Q^{\prime 2})d\ln Q^{\prime 2}$.

The Laplace transform on $\hat{F}_{NS}(\nu,\tau)$ now leads to $f_{NS}(s,\tau)\equiv{\cal L}[\hat{F}_{NS}(\nu,\tau);s]$. On the other hand, it is known that the Laplace transform of convolution factors will yield the ordinary product of the Laplace transform of those factors Block:2010du ; Block:2011xb ; MoosaviNejad:2016ebo . Therefore, by imposing the Laplace transform on Eq. (II), the result would be the ordinary first order differential equation with respect to the $\tau$-variable in Laplace space $s$. The result for the non-singlet distributions $f_{\rm NS}(s,\tau)$ reads

[TABLE]

where, $\Phi_{\text{NS}}^{\text{i}}(\text{i=LO, NLO, NNLO})$ represents the Laplace transform of the splitting kernels at the desired order. The analytical results for the splitting kernels in $x$-space are given in Refs. Curci:1980uw ; Moch:2004pa . From Eq. (4), one can achieve a solution involving a simple form as

[TABLE]

In this equation the splitting function, $\Phi_{\text{NS}}(s)$, can be presented in terms of new expansion parameters as

[TABLE]

The variables $\tau_{2}$ and $\tau_{3}$ in Eq. (6) are given by

[TABLE]

and

[TABLE]

which are $Q^{2}$-dependent.

For the splitting function $\Phi_{\text{NS}}^{\text{LO}}(s)$ in the Laplace $s$-space, the result is typically given by

[TABLE]

where, $\gamma_{E}=0.577216\cdots$ is the Euler constant and $\psi(s)=d\ln\Gamma(s)/ds$ is the digamma function. The corresponding results for the NLO and NNLO splitting kernels are too lengthy so they are presented in appendix A.

Considering Eq. (5) and using Eq. (6) the solution of non-singlet sector of quark densities can be obtained in the Laplace $s$-space. Applying the inverse Laplace transform in two stages one achieves the quark densities in Bjorken $x$-space as a function of $Q^{2}$. Details of this procedure can be found in Ref. Khanpour:2016uxh . Back to Eq. (5), we see that the solution of non-singlet parton densities would be obtained provided that the $f^{0}_{\text{ NS}}(s)$ at the initial scale $Q_{0}$ is known. In this regards, one needs to extract the parton densities at the initial scale $Q_{0}$ which can be done by a global fit over whole available data. One of the reliable techniques to extract the initial parton densities is the Jacobi polynomials approach which will be illustrated in the following sections.

III Jacobi polynomials and non-singlet structure function

In this section, we apply the Jacobi polynomials to convert the results for the non-singlet structure function from the Laplace $s$-space to the known Bjorken $x$-space. Based on this method, we will be able to perform a global fit over all available experimental data to extract free parameters in the proposed form of the parton densities at the initial scale $Q_{0}^{2}$. This procedure will be explained later. In this approach, the results which will be obtained for the parton densities in Laplace $s$-space can be considered as the moment of densities. Therefore, the theoretical perspectives on Jacobi polynomial approach is to extract the non-singlet structure function $F_{2}^{NS}(x,Q^{2})$ from the NNLO analytical solution of non-singlet DGLAP equations in Laplace $s$-space. It should be noted that, the solution of evolved structure functions at any value of $x$ and Q2 could be attained through the method described. This is a crucial feature for the phenomenological implications.

Now, the results in Laplace $s$-space for the non-singlet structure function $F_{2}^{NS}(x,Q^{2})$ are considered as the Laplace $s$-space moments ${\cal M}^{NS}(s,Q^{2})$. These provide the possibility to include DIS data for performing a QCD analysis up to NNLO accuracy. The DIS data contains a wide range of the transferred momentum from $Q_{0}^{2}\gtrsim 2GeV^{2}$ to $Q^{2}\sim 30000GeV^{2}$ where the Jacobi polynomials technique works, reasonably.

To construct the moments of proton structure function, ${\cal M}^{p}_{2}(s,\tau)$, at NNLO accuracy the following combination of parton densities at the valence region $x\geq 0.3$ in Laplace $s$-space (for the non-singlet sector) are required

[TABLE]

where, $C_{2NS}^{(1)}$ and $C_{2NS}^{(2)}$ are the Wilson coefficients at NLO and NNLO accuracies in the Laplace $s$-space, respectively. Their analytical expressions are presented in appendix $A$.

For the non-singlet sector of deuteron structure function, ${\cal M}_{2}^{d}$ (with $d=(p+n)/2$), in Laplace $s$-space the relevant combination of parton densities at NNLO expansion can be written as

[TABLE]

It is now possible to consider the difference of proton and deuteron structure functions which is important to analysis the data in the region $x\leq 0.3$. It reads

[TABLE]

For $x<0.3$, the effect of sea quark densities , $\bar{u}-\bar{d}$ in above equation, can not be neglected. In our calculations, we take from JR14 Jimenez-Delgado:2014twa this distribution at the initial scale Q${}_{0}^{2}$ = 2 GeV2 as,

[TABLE]

For practical purposes, the combinations of $d-\bar{d}$ and $u-\bar{u}$ are also considered as the proton valence densities. These are denoted as $d_{v}$ and $u_{v}$, respectively. The following valence distributions are employed in our analysis at the input scale Q${}_{0}^{2}$ = 2 GeV2 Jimenez-Delgado:2014twa ,

[TABLE]

[TABLE]

where, ${\cal N}_{u}$ and ${\cal N}_{d}$ are the normalization factors which are given by

[TABLE]

where, the B-function stands for the Euler beta function.

Using a variable change as $x=e^{-\nu}$ and considering the Laplace transformation via the following relations

[TABLE]

[TABLE]

the valence distributions in Eqs. (14) and (15) can be presented in the Laplace $s$-space as

[TABLE]

and

[TABLE]

To extract the unknown parameters in the valance distribution, one needs to do a global fit based on the Jacobi polynomial approach. Details of this approach can be found in Refs. Khorramian:2008yh ; Khorramian:2009xz . Here, we review this method briefly. Using this approach, the structure functions can be reconstructed as

[TABLE]

In the above equation, the Jacobi polynomials $\Theta_{n}^{\alpha,\beta}(x)$ have the following expansion

[TABLE]

where, $c_{j}^{(n)}(\alpha,\beta)$ are the coefficients which can be written in terms of the Euler Gamma function. The parameters $\alpha$ and $\beta$ are fixed and taken to be 3 and 0.7, respectively Khorramian:2008yh ; Khorramian:2009xz . As a final point, by considering the weight function $x^{\beta}(1-x)^{\alpha}$ the orthogonality condition for the Jacobi polynomials is given by

[TABLE]

The moment of structure function in Laplace $s$-space, given by last term in Eq. (III), can be written in terms of Eq. (III). The global fit over all available experimental data for the structure functions will provide us the unknown parameters of valance densities. Technical details of the global fit are investigated in the following section.

Numerical results for the free parameters in the $u_{v}$ and $d_{v}$ densities at NNLO and NLO accuracies are listed in Tables. 1 and 2, respectively. They are resulted from non-singlet QCD fit at $Q^{2}_{0}=2\;GeV^{2}$, as would be described in Sec. IV. Although, the uncertainty bands of PDFs are related to the experimental uncertainties but they are also affected by increasing the theoretical precision. As is seen from Figs. 1 and 2 , the error bands for NNLO results are smaller than the NLO one which indicates an enhancement in accuracy of analysis.

IV Global analyses of valence-quarks densities

IV.1 Different features of data sets

The valence PDFs can be determined by fitting to a global data set over 572 data points including a variety range of scattering processes at high energies. In Table. 3, the data sets used in our analysis are listed. The DIS data from BCDMS Benvenuti:1989fm ; Benvenuti:1989rh ; Benvenuti:1989gs and SLAC Whitlow:1991uw as well as NMC Arneodo:1996qe ; Arneodo:1995cq experiments make the required sets. The flavor separation of PDFs at large $x$ is facilitated by these data sets. In our analysis, the DIS data from H1 H1 and ZEUS ZEUS Collaborations are also employed. Additionally, the combined measurements of H1 and ZEUS Collaborations at HERA for the inclusive $e^{\pm}p$ scattering cross sections are applied as the new data sets Aaron:2009aa . For the valence quarks in the region $x\geq 0.3$, we use the data set given for $F_{2}^{p}(x,Q^{2})$ and $F_{2}^{d}(x,Q^{2})$ while for the region $x<0.3$ the data are related to $F_{2}^{\rm NS}(x,Q^{2})=2(F_{2}^{p}(x,Q^{2})-F_{2}^{d}(x,Q^{2}))$.

In order to eliminate the higher twist (HT) effects, it is needed to take different cuts on data sets before analysis. The DIS data cuts are considered for Q${}^{2}>4$ GeV2 and on the hadronic mass $W^{2}>12.5$ GeV2 which perform the required converge in the concerned kinematic region. For the BCDMS and NMC data sets, additional cuts as $y>0.35$ and $Q^{2}>8$ GeV2 should be applied, respectively. Considering the required cuts, we listed in Table. 3 both the required DIS data and the number of data points for each experiment in the fitting process. In the 5th column of Table. 3, taking the additional cuts the number of reduced data points are presented. Consequently, the number of data points are reduced from 467 to 248 for $F_{2}^{p}(x,Q^{2})$ so that for $F_{2}^{d}(x,Q^{2})$ this reduction is from 232 to 159. At last, for $F_{2}^{\rm NS}(x,Q^{2})$ the number of data points are reduced from 208 to 165.

IV.2 $\chi^{2}$-Minimization approach

The best values of fit parameters at NNLO accuracy are extracted by minimizing the $\chi^{2}$ with respect to four unknown parameters in valence-quark distributions (14) and (15), along with $\Lambda_{\rm{\overline{MS}}}^{(4)}$ as the QCD cutoff parameter. As is shown in Tables. 1 and 2, we fix four parameters in the fit from the beginning, i.e. $\gamma_{u},\eta_{u},\gamma_{d}$ and $\eta_{d}$, so there are totally five free parameters remaining to be extracted from the fit process.

The global goodness-of-fit can be done using the following definition for the effective $\chi^{2}$ cor

[TABLE]

where, $n$ labels the different experiments and

[TABLE]

In the above relation, $F_{2,i}^{data}$, $\Delta F_{2,i}^{data}$ and $F_{2,i}^{theor}$ stand for the data value, the data uncertainty including a combination of statistical and systematic errors, and the theoretical value for the $i^{\mathrm{th}}$ data point in the $n^{\mathrm{th}}$ experiment, respectively. Moreover, ${\Delta{\cal N}_{n}}$ in the first parenthesis is the experimental normalization uncertainty and for the data of $n^{\mathrm{th}}$ experiment an overall normalization factor is given by ${\cal N}_{n}$. This factor is extracted through the fit procedure and when it is determined through the first fit we fix its value in next stages. The possible weighting factor denotes by $w_{n}$ which is usually taken to be 1. In our analysis, a relative normalization shift between the different data sets, i.e., ${\cal N}_{n}$, is used for which the normalization uncertainties ${\Delta{\cal N}_{n}}$ are given by the experiments.

In our analysis, the total number of data points are $n^{data}=572$ and there are also five unknown parameters in the fit. The best parametrization of valence-quark densities are obtained using the CERN program library MINUIT CERN-Minuit . The results are listed in Tables. 1 and 2.

IV.3 Determining the input uncertainties

In this subsection, we describe our method to determine the uncertainties of valence-quark PDFs as well as the error propagation from experimental data points. In practice, the uncertainties in global PDF analysis can be achieved by procedures which are well-defined to propagate through the experimental uncertainties on the fitted data points to the PDF uncertainties. For this purpose, we apply the Hessian method (or error matrix approach) Pumplin:2001ct which is based linearly on error propagation and the suitable production of PDF eigenvector. The Hessian method was firstly used in analysis by MRST04 Martin:2003sk , MSTW08 Martin:2009iq , MRST03 Martin:2002aw , CETQ group Dulat:2015mca and its up-to-dated version by CT18 Hou:2019efy , and also in our previous works AtashbarTehrani:2012xh ; Khanpour:2016pph ; Monfared:2011xf ; Khanpour:2012tk . By running the CERN program library MINUIT CERN-Minuit , an error analysis can be done using the Hessian matrix. A simple and efficient method to calculate the uncertainties of parton densities can be obtained by applying this method which is essentially related to diagonalize the covariance matrix. In the Hessian approach, the main assumption is to perform a quadratic expansion of $\chi^{2}_{\text{ global}}$ with respect to the fitted parameters $a_{i}$ about its global minimum. It reads

[TABLE]

In this equation, the elements of the Hessian matrix are denoted by $H_{\rm ij}$ and the number of parameters in the global fit is presented by $p$.

IV.4 Error calculations

The method presented in Refs. Martin:2002aw ; Martin:2003sk ; Pumplin:2001ct ; Martin:2009iq ; AtashbarTehrani:2012xh ; Khanpour:2016pph can be followed to determine the error of calculations. The eigenvectors and eigenvalues of the covariance (or Hessian) matrix are the basic mathematical tools. For practical purposes we need, at first, to have a set of appropriate fit parameters considered in the valence quark densities. They are usually labeled by $a_{i}(s_{\text{min}})$ and required to minimize the $\chi^{2}_{\text{ global}}$. Considering the parton sets as $s_{k}^{\pm}$, it is possible to have an expansion in terms of the eigenvectors and eigenvalues for the variation of parameters around the global minimum, i.e.

[TABLE]

Here, $a_{i}(s_{k}^{\pm})$ and $a_{i}(s_{\text{min}})$ are the same as the abbreviated $a_{i}$ and $a_{i}^{0}$ symbols in Eq. (26). In the above equation, the k${}^{\text{th}}$-eigenvalue and the i${}^{\text{th}}$-component of the orthonormal eigenvectors of Hessian matrix are denoted by $\lambda_{k}$ and $v_{\text{ik}}$, respectively. By adjusting the $t$ variable proportionally, it is possible to set $t=T$ in the quadratic approximation and therefore to get $T^{2}=\Delta\chi^{2}$. The quadratic approximation of Eq. (26) can be tested by considering the dependence of $\Delta\chi^{2}$ on the eigenvector directions for some selected samples. In Fig. 3, $\Delta\chi^{2}$ is plotted as a function of $t$ for all eigenvectors at the NNLO accuracy.

In our NNLO global analysis, the covariance matrix elements for five free parameters are given in Table. 4. The estimation for uncertainties of parton densities which is generally denoted by $F$ can be determined, using the following relation introduced in Martin:2002aw :

[TABLE]

where, $C_{ij}(a)=(H^{-1})_{ij}$ is the covariance matrix and $\Delta\chi^{2}$ is the allowed variation in $\chi^{2}$, as was specified previously. By suitable choice of $\Delta\chi^{2}=1$, which corresponds to one sigma confidence level, and by considering the Hessian (or equivalently covariance) matrix it is now possible to calculate the errors on parton densities.

Now, the valence quark densities can be calculated at higher-Q2 (i.e. $Q^{2}>Q_{0}^{2}$) using the DGLAP evolution equations in Laplace $s$-space. Using the analytical solution, based on the Laplace transform technique as well as Jacobi polynomials approach, the global QCD analysis can be done. In Table. 1, we listed the numerical values of PDFs parameters which are related to the QCD fit for the NNLO non-singlet sector at the input scale Q${}_{0}^{2}=2$ GeV2. After the first minimization, some parameters take values with very small error so that one can ignore their errors and fix them to their central values through the fit procedure. By fixing some parameters, the fit is repeated which could be done easier and faster. Hence some parameters are presented in Tables. 1 and 2 without any error. We can see from Table. 1 that the central values of PDF parameters are rather stable. From the fit, one can find the numerical value for $\alpha_{\rm s}^{\rm N_{f}=4}(Q_{0}^{2})$ as presented in Tables. 1 and 2. This implies the numerical value for QCD cutoff quantity. Using this quantity and considering the threshold matching condition for running coupling constant match , we get $\alpha_{s}(M_{Z}^{2})=0.1147\pm 0.0009$ at the mass scale of $Z$-boson. In Fig. 1, we depicted the $xu_{v}$- and $xd_{v}$-distributions at the energy scale $Q_{0}^{2}=2$ GeV2. This also contains the PDF uncertainties with required confidence level corresponding to $\Delta\chi^{2}=1$ in Eq. (26). In this figure, we have also added the results from other models such MKAM16 MoosaviNejad:2016ebo , NNPDF Ball:2012cx , MMHT14 Harland-Lang:2014zoa , CT14 Dulat:2015mca and up-to-dated results from CT18 PDFs Hou:2019efy to have a more qualitative comparison. To confirm the considerable progress of PDFs at NNLO accuracy, the results of NLO PDFs for up and down valence densities have been also depicted in Fig. 2 at the scale $Q_{0}^{2}=2$ GeV2. The numerical values of fitted parameters at NLO accuracy are listed in Table. 2. As is seen from Tables. 1 and 2, the value of $\chi^{2}/{\rm n.d.f}$ is reduced from 0.920 to 0.906 by increasing the accuracy of corrections from NLO to NNLO. Moreover, in Fig. 2 the small uncertainty raised in our NNLO analysis confirms the validity of NNLO analysis in comparison with the NLO one. Of course, the higher order correction is not the only reason for decreasing the errors and additionally the type of fitting procedure would be also effective on getting the small error bands. This point will be described in the following section.

V Higher twist and target mass corrections

To determine the moment of structure function, it is assumed that in the limit of $Q^{2}\rightarrow\infty$ the mass of target hadron approaches to zero. This assumption would be violated at the intermediate and low values of $Q^{2}$. In these regions, the moment is significantly affected by the power correction of order ${\cal O}(m_{N}^{2}/Q^{2})$ where $m_{N}$ denotes the nucleon mass Kataev:1996vu ; Kataev:1997nc . This correction is known as the target mass correction (TMC). Taking into account the TMC effect leads to change the moment of non-singlet structure function as it follows  Khorramian:2008yh ; Khorramian:2009xz ; Khanpour:2016uxh ; Georgi:1976ve ; Gluck:2006yz ; Steffens:2012jx ,

[TABLE]

In the relevant region where $x<0.8$, the higher powers of $(m_{\rm N}^{2}/Q^{2})^{n}$ (for $n\geqslant 2$) are negligible Georgi:1976ve . By substituting Eq. (V) into Eq. (III), one obtains

[TABLE]

where, ${\cal M}_{2,{\rm TMC}}^{k}(j+1,Q^{2})$ is the moment of concerned structure function in the Laplace $s$-space. This term contains the TMC effect and is given by Eq. (V) where the corrections due to higher powers (for $n\geqslant 2$) are neglected.

The effect of higher twist (HT) correction is also significant at large values of $x$ and moderate $Q^{2}$ Abt:2016vjh ; Jimenez-Delgado:2013boa ; Leader:2006xc ; Nath:2016phi ; Wei:2016far where the TMC effect is also considerable. Consequently, inclusion of higher twists has its own importance to analysis the parton densities.

The required kinematic cuts, containing the HT effect for $F_{2}^{p}(x,Q^{2})$ and $F_{2}^{d}(x,Q^{2})$, are as $Q^{2}\geq 4$ GeV2 and $W^{2}\geq 12.5$ GeV2 where $W^{2}$ denotes the hadronic invariant mass. The cuts can be extended to the kinematic region $4<W^{2}<12.5{\rm GeV}^{2}$. Therefore, an extrapolation is needed to this region in our QCD fit. It should be noted that, the $W^{2}$ is related to the $Q^{2}$-variable through $W^{2}=(1/x-1)Q^{2}+m^{2}_{n}$, where $m_{n}$ is the nucleon mass. From this relation, it can be seen that the required kinematic cut is such that for $W^{2}\geq 12.5$ GeV2, the $Q^{2}$ is also increasing. Consequently, the TMC and HT effects can be ignored in this region. Therefore, as is seen from Figs. 4 and 5, three plots labeled as QCD, QCD+TMC and QCD+TMC+HT show the same behavior.

Concerning the higher twist effect it should be noted that, the parametrization form of higher twist contributions are practically independent of the leading twist one. This parametrization is typically given by polynomial functions of $x$. For the DIS data where the power corrections in the concerned region can not be neglected, the corrections are defined by an ansatz which is completely motivated from phenomenological points of view and it reads Blumlein:2008kz ; Blumlein:2006be ; Martin:1998np ; Martin:2003sk

[TABLE]

where $F_{2}^{\rm TMC}(x,Q^{2})$ is given by Eq. (V). The operation ${\cal O}_{\rm TMC}[\cdots]$ in the above equation refers to the target mass correction while the twist-2 contribution is considered for the desired structure function. It should be noted that, the $h(x,Q^{2})$-coefficient is determined at individual intervals of $x$ and Q2 but finally the averaged result over Q2 is taken into account. Therefore, the $x$-dependence of $h(x)$ as the higher twist contribution could be defined as Blumlein:2008kz

[TABLE]

To achieve a flexible higher twist contribution during the data analysis, the above choice of $h(x,Q^{2})$ would be effective. As was mentioned, the cuts $Q^{2}\geq 4\,{\rm GeV}^{2}$ and $4<W^{2}<12.5\,{\rm GeV}^{2}$ should be considered for the higher twist QCD analysis of the non-singlet data. The value of $\Lambda_{\rm QCD}$ and the free parameters of valence PDFs and function $h(x)$ could be extracted through the fit procedure to all available data. Our strategy is to perform a fit for the PDF parameters at large values of $W^{2}$ and then to impose the TMC effect to the structure functions. Following that a HT-term is added to the calculation while the PDFs parameters are fixed. The free parameters of HT-term are determined during the second stage of fit procedure using the lower values of $W^{2}$. Our reason for this strategy of fitting is to indicate the ability of Laplace transformation technique at the NNLO expansion to extract the non-singlet structure function. The fitted result for the HT term is listed in Table. 5.

VI Results and Discussion

In this work, we discussed about the results arising from the PDF fits. In this regard, we analyzed the global fitting on parton densities in Sec. IV. They are mostly related to Figs. 1, 2 and 3 and also to the corresponding Tables. 1, 2 and 4. In Fig. 3, we indicated the results of $\Delta\chi^{2}$ with respect to the $t$-variable considering all five eigenvectors. As can be seen, a very good description of data is provided by considering our theoretical predictions which are based on the analytical solutions resulted from the Laplace transform technique and Jacobi polynomials approach.

The rest of results are related to the TMC and HT effects which have been described in Sec. V, in detail. In this section, Fig. 4 involves the data for inclusive proton structure functions ($F_{2}^{p}$) from BCDMS and SLAC experiments. A comparison between available experimental data and our NNLO fit, as a function of Q2 at approximately constant values of $x$, has been also done. The effects of TMC and HT corrections have been also included so that a considerable agreement between the experimental data and the theoretical predictions is achieved especially at low energy scales. Briefly, our strategy includes two steps. In the first step, we performed the required analysis to get the best values for the PDF parameters. At the second step, the TMC and HT effects have been included. They affected the results, specially at low energy scales. Final results, which are shown in Figs. 4 and 5, indicate good agreements with the experimental data for the available range of energy scales. A more detailed comparison between the theoretical predictions of our NNLO fit and the data for the deuteron structure function ($F^{d}_{2}$) has been depicted in Fig. 5. The results have contained again the effect of TMC and HT corrections.

In Fig. 6, a comparison between the structure function for $(d-p)$, i.e. $F^{\rm NS}_{2}$, and the data from BCDMS and NMC experiments has been shown. For better presentation, the data have been scaled by $c=0.2$, $c=0.4$ and $c=0.6$ for $x=0.275$, $x=0.225$ and $x=0.18$, respectively.

Here, we just remind that the agreement between the theoretical predictions and the experimental data becomes significantly better when one employs the effect of TMC and HT corrections. This point can be seen clearly in Figs. 4 and 5. In this regard, we have also got an excellent agreement between the theoretical prediction for the structure functions $F^{p}_{2}$ and $F^{d}_{2}$ and the available experimental data which included the TMC and HT effects. The data included different ranges of Q2 and $x$-values.

The technique of Laplace transformation can be also employed in analyzing the Drell-Yan data which contain antiquark distributions. This type of distribution has not been yet considered in our recent analysis. It would be also valuable if we can improve the method of our fitting by including the TMC calculation and the parameters of HT term together with the PDF parameters, simultaneously. These projects are in progress and remain as our future works.

Acknowledgment

S. A. T is thankful to the School of Particles and Accelerators, Institute for Research in Fundamental Sciences (IPM) for financial support of this project. S. M. M. N and A. M are also grateful the Yazd university to provide the required facilities to do this project.

Appendix A

In this appendix, we present the Laplace transforms for the NLO and NNLO splitting functions as well as the Wilson coefficient functions. They are related to the non-singlet sectors of structure function used in Eq. (6) and Eqs. (III)-(III). In deriving the NLO and NNLO Wilson coefficients $C_{2NS}^{(1)}$ and $C_{2NS}^{(2)}$ in Laplace $s$-space, we use the corresponding results at $x$-space given in Refs. Moch:1999eb ; Vermaseren:2005qc . The NLO and NNLO splitting functions in $s$-space, i.e. $\Phi_{(NS)}^{\rm NLO}$ and $\Phi_{(NS)}^{\rm NNLO}$, are calculated from the corresponding results in $x$-space which are given in Refs. Curci:1980uw ; Moch:2004pa . The obtained results for splitting functions and the Wilson coefficients in $s$-space are as following:

[TABLE]

[TABLE]

[TABLE]

[TABLE]