Dynamical parton distributions from DGLAP equations with nonlinear corrections

Determination of proton parton distribution functions is present under the dynamical parton model assumption by applying DGLAP equations with GLR-MQ-ZRS corrections. We provide two data sets, referred as IMParton16, which are from two different nonperturbative inputs. One is the naive input of three valence quarks and the other is the input of three valence quarks with flavor-asymmetric sea components. Basically, both data sets are compatible with the experimental measurements at high scale ($Q^2>2$ GeV$^2$). Furthermore, our analysis shows that the input with flavor-asymmetric sea components better reproduce the structure functions at high $Q^2$. Generally, the obtained parton distribution functions, especially the gluon distribution function, are the good options of inputs for simulations of high energy scattering processes. The analysis is performed under the fixed-flavor number scheme for $n_f=$ 3, 4, 5. Both data sets start from very low scales around 0.07 GeV$^2$, where the nonperturbative input is directly connected to the simple picture of quark model. These results may shed some lights on the origin of the parton distributions observed at high $Q^2$.


I. INTRODUCTION
Hadrons are the complex systems consisting of quarks and gluons, which makes a long and continuous way to precisely understand the hadron structure. Thanks to the collinear factorization theorem [1][2][3] in quantum chromodynamics (QCD), the calculation of high energy hadron collision becomes much straightforward. The calculation is the product of the calculable hard process and the incalculable soft part which is absorbed into the parton distribution functions (PDFs). Although incalculable so far, parton distribution functions are universal coefficients which can be determined by the experiments conducted worldwide. Moreover there are some models [4][5][6][7] and Lattice QCD (LQCD) calculations [8][9][10] which try to predict/match the PDFs of proton. PDFs in wide kinematic ranges of Q 2 and x is an important tool to give some theoretical predictions of high energy hadron collisions and simulations of expected interesting physics in modern colliders or JLab experiments of high luminosity.
Determination of PDFs of proton attracts a lot of interests on both theoretical and experimental sides. To Date, the most reliable and precise PDFs data comes from the global QCD analysis of experimental data. There have been a lot of efforts and progresses achieved on this issue [11][12][13][14][15][16][17][18]. In the global * email: rwang@impcas.ac.cn † email: xchen@impcas.ac.cn analysis, firstly, the initial parton distributions at low scale Q 2 0 ∼ 1 GeV 2 , commonly called the nonperturbative input, is parameterized using complicated functions with many parameters. Given the nonperturbative input, the PDFs at high Q 2 are predicted by using DGLAP equations from QCD theory. Secondly, the nonperturbative input is determined by comparing the theoretical predictions to the experimental data measured at high scale. This procedure is usually chosen to be the least square regression method. Finally, PDFs in a wide kinematic range is given with the obtained optimized nonperturbative input. Although a lot of progresses have been made, the gluon distribution at small x is still poorly estimated, which has large uncertainties [19,20]. Even worse, the gluon distributions from different collaborations exhibit large differences. Gluon distribution needs to be more quantitative in terms with a number of physics issues relating to the behavior of it [21][22][23].
PDFs at low resolution scale is always confusing since it is in the nonperturbative QCD region. However it is related to the nucleon structure information measured at high resolution scale. Therefore the nonperturbative input gives some valuable information of the nucleon. Besides the powerful predictions of the QCD theory, other fundamental rules of hadron physics should also be reflected in the nonperturbative input. How does the PDFs relate to the simple picture of the proton made up of three quarks? In the dynamical parton model [11,12,[24][25][26], the input contains only valence quarks, valence-like light seas and valence-like gluon, which is consistent with the dressed constitute quark model. All sea quarks and gluons at small x are dynamically produced. In the dynamical parton approach, the gluon and sea quark distributions are excellently constrained by the experimental data, since there are no parametrizations for input dynamical parton distributions. Parton radiation is the dynamical origin of sea quarks and gluons inside the proton. It is also worthwhile to point out that the valence-like input and PDFs generated from it are positive. In some analysis of MSTW [27][28][29], the negative gluon density distributions are allowed for the nonperturbative input in order to fit the small-x behavior observed at high scale. MMHT2014 PDFs [18] supersede the MSTW2008 PDFs, which make some changes to the parameterisation to stop negative distribution in some region of x, and include LHC, updated Tevatron data and the HERA combined H1 and ZEUS data on the total and charm structure functions in the global analysis. The dynamical parton model is developed and extended to even low scale around Q 2 0 ∼ 0.06 GeV 2 in our previous works [30,31]. The naive nonperturbative input [30] with merely three valence quarks are realized, which is the simplest input for the nucleon. In the later research [31], we composed a nonperturbative input which consists of three valence quarks and flavor-asymmetric sea components, and extracted the flavor-asymmetric sea components from various experimental data measured at high Q 2 . The flavor-asymmetric sea components here refer to the sea quark distributions generated not from the QCD evolution but from the complicated nonperturbative QCD mechanisms. In terms of the interpretation of the nonperturbative input, the extended dynamical parton model gives the clearest physics picture. This work is mainly based on our previous works [30,31]. The extended dynamical parton model is taken in the analysis.
DGLAP equations [32][33][34] based on parton model and perturbative QCD theory successfully and quantitatively interpret the Q 2 -dependence of PDFs. It is so successful that most of the PDFs are extracted by using the DGLAP equations up to now. And the common way of improving the accuracy of the determined PDFs is to apply the higher order calculations of DGLAP equations. However there are many QCD-based evolution equations and corrections to DGLAP equations [35][36][37][38][39] being worked out. It is worthwhile to apply new evolution equations in the global analysis. There are some pioneering works [30,31,40] trying to reach this aim. In this work, DGLAP equations with GLR-MQ-ZRS corrections are taken to do the global analysis.
The main purpose of this study is to give purely dynamical gluon distributions (g(x, Q 2 0 ) = 0), which is expected to be more reliable at small x. The second purpose is to con-nect the quark model picture of proton to the QCD description at high energy scale. The aim is to resolve the origin of sea quarks and gluons at high resolution scale. The third purpose is to understand the QCD dynamics of parton radiation and parton recombination. We want to quantify the strength of GLR-MQ-ZRS corrections by determining the value of parton correlation length R.
The organization of this paper is as follows. Section II lists the experimental data we used in the analysis. Section III discusses the QCD evolution equations, which is the most important tool to evaluate the PDFs. The nonperturbative input inspired by quark model and other nonperturbative effects are discussed in Sec. IV. The other details of the QCD analysis are explained in Sec. V. Section VI shows the results of the global fits and the comparisons of the obtained PDFs to experimental measurements and other widely used PDF data sets. Section VII introduces the IMParton package which gives the interface of the obtained PDFs. Finally, a simple discussion and summary is given in Sec. VIII.

II. EXPERIMENTAL DATA
The deeply inelastic scattering (DIS) of charged leptons on nucleon has been the powerful tool to study nucleon structure for a long time. The quark structure of matter is clearly acquired by decades of measurements starting from the late 1960s with the lepton probes interacting mainly through the electromagnetic force. The DIS data of leptons is so important that we include only the DIS data in this work. The structure function F 2 (x, Q 2 ) data used in this analysis are taken from SLAC [41], BCDMS [42], NMC [43], E665 [44] and HERA (H1 and ZEUS) [17,45] collaborations.
In order to make sure the data is in the deep inelastic region, and to eliminate the contributions of nucleon resonances, two kinematic requirements shown in Eq. (1) are performed to select the experimental data.
For the neutral-current DIS, the contribution of the Z-boson exchange can not be neglected at high Q 2 . Therefore we compose another kinematic cut to reduce the influence of the Z-boson exchange contribution, which is shown in Eq. (2). The Z-boson exchange contribution is of the order ∼ 1% at Q 2 = 1000 GeV 2 .
SLAC was the first to perform the fixed-target DIS experiments. The SLAC data we used is from the reanalysis of a series of eight electron inclusive scattering experiments conducted between 1970 and 1985. The reanalysis procedure implement some improved treatments of radiative correction and the value of R = σ L /σ T . The minus four-momentum transfer squared Q 2 of SLAC experiments are not big (Q 2 ≤ 30 GeV 2 ), and the x is mainly at large x (0.06 ≤ x ≤ 0.9) because of the relative low beam energy. The target mass correction (TMC) should not be ignored for the SLAC data, because of the low Q 2 and large x. In this work, the formula of TMC [46,47] is taken as, with r = 1 + 4x 2 M 2 /Q 2 , and ξ the Nachtmann variable defined as 2x/(1 + 1 + 4x 2 M 2 /Q 2 ). Compared to the later experiments, the uncertainties of the structure functions and the absolute normalization of SLAC data are big.   The precise measurements of the structure function F 2 was followed by the experiments at CERN, Fermilab, and HERA at DESY. Both BCDMS and NMC data are collected from the muon-proton DIS with CERN SPS muon beam but with radically different detectors. The BCDMS data are taken at beam energy of 100, 120, 200 and 280 GeV, and the NMC data are taken at beam energy of 90, 120, 200 and 280 GeV. The absolute normalization for the NMC data was based on an empirical data model motivated basically by leading order QCD calculations. Therefore we should fit the NMC normalization factors for each incident beam energy. The H1 and ZEUS data at HERA span a wide kinematic region of both Q 2 and x. The small x information of the structure function primarily comes from the HERA data. The HERA data we used is the combined analysis of H1 and ZEUS experiments. The normalization uncertainty in this data is 0.5%. A complementary set of the inclusive HERA data was obtained by the H1 Collaboration in the run with a reduced collision energy. These data are particularly sensitive to the structure function F L and thereby to the small-x shape of the gluon distribution.
Finally, the kinematic coverage of the charged leptonproton DIS data is shown in Fig. 1. The kinematic of all the data covers 3 magnitudes in both x and Q 2 . Since the SLAC and the NMC data distribute from relatively low Q 2 , the target mass corrections are applied when comparing theoretical calculations to these data. All the normalization factors of the experimental data are fitted in the analysis except for the combined data of H1 and ZEUS, as the normalization uncertainties are not small for other data.

III. NONLINEAR CORRECTIONS TO DGLAP EQUATIONS
DGLAP equations [32][33][34] is the important and widely used tool to describe the Q 2 dependence of quark and gluon densities. The equations are derived from the perturbative QCD theory using the quark-parton model instead of the rigorous renormalization group equations, which offers a illuminating interpretation of the scaling violation and the picture of parton evolution with the Q 2 . The DGLAP equations are written as, in which P q i q j , P q i g , P gq j and P gg are the parton splitting functions [34]. The prominent characteristic of the solution of the equations is the rising sea quark and gluon densities toward small x. The QCD radiatively generated parton distributions at small x and at high Q 2 are tested extensively by the measurements of hard processes at modern accelerators.
The most important correction to DGLAP evolution is the parton recombination effect. The theoretical prediction of this effect is initiated by Gribov, Levin and Ryskin (GLR) [35], and followed by Mueller, Qiu (MQ) [36], Zhu, Ruan and Shen (ZRS) [37][38][39] with concrete and different methods. The number densities of partons increase rapidly at small x. At some small x, the number density become so large that the quanta of partons overlap spatially. One simple criterion to estimate this saturation region is x f g (x, Q 2 ) ≥ Q 2 R 2 p , with R p the proton radius. Therefore the parton-parton interaction effect becomes essential at small x, and it expected to stop the increase of the cross sections near their unitarity limit. In ZRS's work, the time-ordered perturbative theory (TOPT) is used instead of the AGK cutting rules [37]. The corrections to DGLAP equations are calculated in the leading logarithmic (Q 2 ) approximation, and extended to the whole x region, which satisfy the momentum conservation rule [38].
In this analysis, DGLAP equations with GLR-MQ-ZRS corrections are used to evaluate the PDFs of proton. The GLR-MQ-ZRS corrections are very important to slow down the parton splitting at low Q 2 < 1 GeV 2 . Up to date, ZRS have derived all the recombination functions for gluon-gluon, quark-gluon and quark-quark processes [37][38][39]. Our previous work finds that the gluon-gluon recombination effect is dominant [31], since the gluon density is significantly larger than the quark density at small x. Therefore, we use the simplified form of DGLAP equations with GLR-MQ-ZRS corrections, which is written as, for the flavor non-singlet quark distributions, for the dynamical sea quark distributions, and for the gluon distribution, in which the factor 1/(4πR 2 ) is from the normalization of the two-parton densities, and R is the correlation length of the two interacting partons. In most cases, R is supposed to be smaller than the hadron radius [36]. Note that the integral terms as 1/2 x in above equations should removed when x is larger than 1/2. Σ in Eq. (7) is defined as The splitting functions of the linear terms are given by DGLAP equations, and the recombination functions of the nonlinear terms are written as [38],

IV. QUARK MODEL AND NONPERTURBATIVE INPUTS
Quark model achieved a remarkable success in explaining the hadron spectropy and some dynamical behaviors of high energy reactions with hadrons involved. Quark model uncovers the internal symmetry of hadrons. Moreover, it implies that the hadrons are composite particles containing two or three quarks. According to the quark model assumption, the sea quarks and gluons of proton at high Q 2 are radiatively produced from three valence quarks. There are some model calculations of the initial valence quark distributions at some low Q 2 0 from MIT bag model [4,5], Nambu-Jona-Lasinio model [6] and maximum entropy [7] estimation.
Inspired by the quark model, an ideal assumption is that proton consists of only three colored quarks at some low scale Q 0 . This assumption results in the naive non-perturbative input -three valence quarks input. At the input scale, the sea quark and gluon distributions are all zero. This thought is widely studied soon after the advent of QCD theory [26,48,49]. The initial scale of the naive nonperturbative input is lower than 1 GeV 2 , since gluons already take comparable part of the proton energy at Q 2 = 1 GeV 2 . To properly evolve the naive nonperturbative input should be considered at such low Q 2 . Partons overlap more often at low Q 2 because of the big size at low resolution scale. In our analysis, the recombination corrections are implemented.
In the dynamical PDF model, all sea quarks and gluons at small x are generated by the QCD evolution processes. Global QCD analysis based on the dynamical PDF model [11,12,[24][25][26] reproduced the experimental data at high Q 2 with high precision using the input of three dominated valence quarks and valence-like components which are of small quantities. Partons produced by the QCD evolution are called the dynamical partons. The input scale for the valence-like input is aroud 0.3 GeV 2 [11,12] and the evolution of the valence-like input is performed with DGLAP equations. In our works, the dynamical PDF model is developed and extended to even low Q 2 [30,31]. The naive nonperturbative input is realized in our approach. The input of valence quarks with flavor-asymmetric sea components is also investigated and found to be a rather better nonperturbative input. The flavor-asymmetric sea components here refer to the intrinsic sea quarks in the light front theory [50,51] or the connected sea quarks in LQCD [52][53][54], or the cloud sea in the π cloud model [55][56][57]. Although there are different theories for the flavor-asymmetric sea components, the flavor-asymmetric sea components are generated by the nonperturbative mechanisms. These types of sea quarks are completely different from the dynamical sea quarks. In this analysis, the evolutions of the flavor-asymmetric sea components obey the equation for the non-singlet quark distributions.
In this work, we try to use two different inputs. One is the naive nonperturbative input and the other is the three valence quarks adding a few flavor-asymmetric sea components. For convenience, three valence quarks input is called input A, and the one with flavor-asymmetric sea components is called input B in this paper. Accordingly, PDFs from inputs A and B are called data set A and data set B respectively. The simplest function form to approximate valence quark distribution is the time-honored canonical parametrization f (x) = Ax B (1 − x) C , which is found to well depict the valence distribution at large x. Therefore the parameterization of the naive input is written as, with zero sea quark distributions and zero gluon distribution. One proton has two up valence quarks and one down valence quark. Therefore we have the valence sum rules for the nonperturbtive inputs, For the naive input, the valence quarks take all the momentum of proton. Hence, we have the momentum sum rule for valence quarks in the naive input, With above constraints, there are only three free parameters left for the parametrizations of the naive input. The naive in-put (Eq. (9)) is the simplest nonperturbative input for proton, which simplifies the nucleon structure greatly. For input B, the parametrizations of valence quarks and the valence sum rules are the same. For simplicity, the parameterizations of the flavor-asymmetric sea components in input B are given by, This parameterizations easily predict thed −ū difference. The dynamical sea quark and gluon distributions are all zero for input B. With the flavor-asymmetric sea components, the momentum sum rule for input B is modified as follows, In order to determine the quantity of the flavor-asymmetric sea components with accuracy, the following constraint Eq. (14) from E866 experiment [58] is taken in this analysis.
Therefore, there are only 7 free parameters left for the parametrization of input B. For better discussion on the quantity of flavor-asymmetric sea, we define δ the momentum fraction of the flavor-asymmetric sea components, One last thing about the nonperturbative input is the input scale Q 0 . According to the naive nonperturbative input, the momentum fraction taken by valence quarks is one. By using QCD evolution for the second moments (momentum) of the valence quark distributions [59] and the measured moments of the valence quark distributions at a higher Q 2 [12], we get the specific starting scale Q 0 = 0.253 GeV for LO evolution (with Λ QCD = 0.204 GeV for f = 3 flavors). This energy scale is very close to the starting scale for bag model PDFs which is 0.26 GeV [60]. In all, the initial scale Q 0 depends on the running coupling constant and the experimental measurements at high Q 2 . We are sure that the initial scale Q 0 for the naive input is close to the pole (Λ QCD ) of coupling constant. In this analysis, the initial scale Q 0 is viewed as a free parameter which can be determined by experimental data.

V. QCD ANALYSIS
The running coupling constant α s and the quark masses are the fundamental parameters of perturbative QCD. In fact these parameters can be determined by the DIS data at high Q 2 . However these fundamental parameters are already determined by a lot of experiments. Hence there is no need to let these parameters to be free. The running coupling constant we choose is in which β 0 = 11 − 2n f /3 and Λ 3,4,5,6 LO = 204, 175, 132, 66.5 MeV [12]. For the α s matchings, we take m c = 1.4 GeV, m b = 4.5 GeV, m t = 175 GeV.
The fixed flavor number scheme (FFNS) is used to deal with heavy quarks in this work. In this approach, the heavy quarks (c, b and t) will not be considered as massless partons within the nucleon. The number of active flavors n f in the DGLAP evolution and the corresponding Wilson coefficients is fixed at n f = 3 (only u, d and s light quarks). The heavy quark flavors are entirely produced perturbatively from the initial light quarks and gluons. The FFNS predictions agree with the DIS data with excellence [12,24]. In this analysis, only charm quark distribution is given, since bottom and top distributions are trivial. The charm quark distribution comes mainly from the gluon distribution through the photon-gluon fusion subprocesses as γ * g → cc, γ * g → ccg and γ * q(q) → ccq(q) [61,62]. The LO contribution of charm quarks to the structure function [61,62] is calculated in this analysis.
The flavor-dependence of sea quarks is an interesting finding in the nucleon structure study [63]. As discussed in Sec. III, the flavor-asymmetric sea componentsū AS andd AS result in thed −ū difference naturally. As found in experiments [64,65] and predicted by the LQCD [54], the strange quark distribution is lower than the up or down quark distribution. In order to reflect the suppression of strange quark distribution, the suppression ratio is applied ass = R(ū DS +d DS )/2 with R = 0.8 [31,54].ū DS +d DS here denotes the dynamical sea quarks. In this approach, the strange quarks are all dynamical sea quarks without any intrinsic components.
The least square method is used to determine the optimal parameterized nonperturbative input. Using DGLAP evolution with recombination corrections, the χ 2 function is calculated by the formula, where N e is the number of data points in experiment e, D i is a data in a experiment, T i is the predicted value from QCD evolution, and σ i is the total uncertainty combing both statistical and systematic errors.

VI. RESULTS
Two separate fits are performed for input A, which consists only three valence quarks. One of them is the fit to all x range (Fit 1) and the other is to fit the data excluding the region of 2 × 10 −3 < x < 0.15 (Fit 2). The results of the fits are listed in Table I. The obtained input valence quark distributions from Fit 1 and Fit 2 are expressed as The initial scale Q 0 and the parton correlation length R for parton recombination are shown in Table II. The obtained R values are smaller than the proton radius, which are consistent with the previous studies [30,36]. In order to justify the importance of parton-parton recombination corrections, we also performed a global fit using DGLAP equations without GLR-MQ-ZRS corrections to the experimental data in the range of x ≤ 2 × 10 −3 or x ≥ 0.15, as a baseline. The obtained χ 2 /N is 25349/1319 = 19.3, and the input scale is Q 0 = 361 MeV. The quality of the fit is bad if we use DGLAP equations without parton-parton recombination corrections, because parton splitting process only generates very steep and high parton distributions at small x [30]. Parton-parton recombination corrections can not be neglected if the evolution of PDFs starts from very low resolution scale. The obtained χ 2 /N is big for input A, especially in the case of Fit 1. Basically, the predicted F 2 structure function gives the similar shape as that measured in experiments, which are shown in Fig. 2(a). However it fails in depicting the experimental data in details around x = 0.02. The experimental data are obviously higher than Fit 2 in the intermediate x region, which is demonstrated clearly in Fig. 2(b). It is interesting to find that the PDFs generated from three valence quarks input miss a peak-like component in the transition region from valence-domain to sea-domain. Three valence quarks input needs to be modified and developed. This discrepancy is expected to be removed by the intrinsic light quarks or cloud sea quarks or connected quarks. In order to get reliable valence quark distributions, the experimental data in the region of 2 × 10 −3 < x < 0.15 should be excluded in the global fit, since the discrepancy around x = 0.02 distorts the optimal three valence quarks input from the analysis. This is the reason why we performed Fit 2 to input A. Fit 2 is in excellence agreement with the experimental data at both large x and small x, which are shown in Fig.  2. Quarks at small x ( 10 −3 ) are mainly the dynamical sea quarks. Generally, our obtained valence quark distributions and the dynamical sea quark distributions are consistent with the experimental observables.
For input B, we performed a fit to the data in all x range (Fit 3). The quality of the fit improves greatly compared to input A, which is shown in Table I. The additional flavorasymmetric sea components are important to remove the discrepancy around x = 0.02. The obtained input is shown in Eq. (19). So far, we have introduced the simplest parametrization for flavor-asymmetric sea components. We argue that more complex parametrization will further improve the result. The total momentum δ carried by flavor-asymmetric sea components at the input scale is obtained to be 0.1. The obtained parameters Q 0 and R are shown in Table II, which is close to that of Fit 1 and Fit 2. The determined input scales are close to the simple theoretical estimation 0.253 GeV as discussed in Sec. IV. The obtained normalization factor of SLAC data is 1.007. The obtained normalization factors of NMC data at beam energy of 90, 120, 200, 280 GeV are 1.07, 1.08, 1.07, and 1.04 respectively. The normalization factors of NMC data by ABM11 global analysis [15] are also large than one. The obtained normalization factor of E665 data is 1.09. The obtained normalization factor of H1 data is 1.02. The obtained normalization factors of BCDMS data at beam energy of 100, 120, 200, 280 GeV are 1.02, 1.01, 1.007, and 1.01 respectively.
The predictions of x-dependence of structure functions at different Q 2 are shown in Figs. 3, 4 and 5 with the experimental data. Our obtained PDFs agree well with the experimental measurements in a wide kinematical range at high resolution scale. The evolutions of F 2 structure function with Q 2 and the comparisons with the experimental data are shown in Fig.  6. Parton distribution functions generated from the the valence quarks and the flavor-asymmetric sea components at the nonperturbative region are consistent with the experimental measurements at high Q 2 > 2 GeV 2 in whole x region with the application of GLR-MQ-ZRS corrections to the standard DGLAP evolution. The experimental data favors some intrinsic components in the nonperturbative input besides three valence quarks.     tribution is close to that of GRV98 and MSTW08, but it is higher than that of CT10. One thing we need to point out is that our gluon distributions are purely dynamically produced in the QCD evolution. We argue that this gluon distribution is more reliable since no arbitrary parametrization of input gluon distribution is involved.  [12], MSTW08 [13] and CT10 [14] at Q 2 = 20 GeV 2 .
Our predicted difference betweend andū are shown in Fig  9. Since the up and down dynamical sea quarks are produced from the gluon splitting, their distributions are the same. The flavor asymmetry between up and down sea quarks are merely from the flavor-asymmetric sea components in this approach. The parametrization of the flavor-asymmetric sea components in this work basically can reproduce the observedd-ū difference observed in Drell-Yan process. Note that the E866 data is not included in the global analysis.
Our predicted strange quark distribution at Q 2 = 2.5 GeV 2 are shown in Fig 10 with the recent reanalysis data by HER-MES collaboration and other widely used parton distribution functions. The predicted strange quark distribution describes the experimental data well, and are consistent with the other PDFs. Our strange quark distribution is purely dynamically generated, since there is no strange quark component in the parameterized nonperturbative input. Compared to the up and down dynamical sea quark distributions, the dynamical strange quark distribution is suppressed in our approach. The suppression of the strange sea quark distribution is not hard to understand, because the current mass of the strange quark is  [12], MSTW08 [13] and CT10 [14] at Q 2 = 20 GeV 2 . The predictedū −d differences from the QCD analysis of only DIS data are shown. The Drell-Yan data from E866 experiment [58] are also shown for comparison. much heavier than that of the up or down quark. This kind of suppression are supported by the LQCD calculation. Fig. 11 shows the comparison of the charm quark distributions to the measurements by H1 and ZEUS Collaborations. The charm quark distributions are based on LO calculation of photon-gluon fusions. This method dealing with the charm quark distributions is also used in the global analysis by GRV95 and GRV98. Although it is a simple calculation un-der the FFNS, the calculation of photon-gluon fusion subprocesses basically reproduced the experimental measurements of the charm quark contribution to F 2 structure function.  In our approach, parton distribution functions at very low Q 2 are also given. We extend the input scale form Q 2 0 = 1 GeV 2 down to Q 2 0 = 0.1 GeV 2 . Our valence quark distribution at low Q 2 are shown in Fig. 12. The valence quark distributions are obviously high at large x. Fig. 13 shows the gluon distributions at low Q 2 . The gluon distributions are Regge-like and positive at even extremely low scale. On the issue of gluon distribution, the prominent advantage of the extended dynamical parton model is that there is no negative gluon density at any resolution scale, no matter how small the Q 2 is.

VII. IMPARTON PACKAGE
We provide a C++ package named IMParton to access the obtained PDFs in the wide kinematic range, in order to avoid the complicated QCD evolution with GLR-MQ-ZRS corrections and make the practical applications of the PDFs easier. The package is now available from us via email, the WWW[68], or download by the git command [69]. Two data sets of the global analysis results, called data set A (Fit 2 result) and data set B (Fit 3 result), are provide by the package. Data set A is from the three valence quarks nonperturbative in-put, and data set B is from the nonperturbative input of three valence quarks adding flavor-asymmetric sea quark components, as discussed in Sec. IV.
The package consists of a C++ class named IMParton which gives the interface to the PDFs. IMParton has a method IMParton::setDataSet(int setOption), which let the user choose data set A or data set B via setDataSet(1) or set-DataSet(2) respectively; The most important method of IM-Parton is IMParton::getPDF(int Iparton, double x, double Q2), which is the method called by users to get the PDF values. Iparton set as -4, -3, -2, -1, 0, 1, 2, 3, 4 corresponds to gettinḡ c,s,d,ū, gluon, u, d, s, c quark/gluon distribution functions respectively. The given PDF values come from the quadratic interpolations of the table grid data calculated by the DGLAP equations with GLR-MQ-ZRS corrections. The table grids are generated in the kinematic range of 10 −6 < x < 1 and 0.125 < Q 2 < 2.68 × 10 8 GeV 2 . The PDF values outside of the grid range are estimated using some sophisticated and effective extrapolation methods. The relative uncertainty of the interpolation is less than 1% in the kinematical range of 10 −6 < x < 0.9.

VIII. DISCUSSIONS AND SUMMARY
We composed some naive nonperturbative inputs inspired by the quark model and some other nonperturbative QCD models at very low Q 2 . By using DGLAP equations with GLR-MQ-ZRS corrections, PDFs generated from these nonperturbative inputs are consistent with various experiments. The obtained gluon distribution is purely dynamically produced, without even the valence-like gluon distribution. The dynamical parton distributions generated in this approach expect to have small bias as a result of the strict theoretical constraints of the method. A C++ package named IMParton is introduced to interface with the obtained PDFs. Two PDF data sets are provided. One is from the three valence quarks input, and the other is from three valence quarks with a few flavorasymmetric sea components. The obtained PDFs can be justified and updated with further investigations of many other hard processes, such as the Drell-Yan process, the inclusive jet production and the vector meson production.
By the global analysis, we find that the quark model on the proton structure has some interesting and good results. The three valence quarks can be viewed as the origin of the PDFs observed at high Q 2 . Our analysis also shows that the nonperturbatvie QCD effects beyond quark model are also needed to reproduce the experimental data in details. By adding the flavor-asymmetric sea components the quality of the global QCD fit improves significantly. This is a clear evidence of the other nonperturbative parton components of proton beyond the quark model [50][51][52][53][54][55][56][57]. It is interesting to know the fact that the sea quarks and gluons are mainly from the parton radiations of three valence quarks predicted by the quark model. However there are more degrees of freedom inside proton which needs the interpretations of the QCD theory in the future.
The nonlinear effects of parton-parton recombinations are important at low Q 2 and small x. Without the recombination processes, the splitting processes generate much steep and large parton densities because of the long evolution distance from extremely low scale. At low Q 2 , the strength of recombination processes are comparable to that of the parton splitting processes. Thus the recombinations slow down enormously the fast splitting of partons at very small x. The preliminary results show that the parton distribution measurements at high Q 2 are directly connected to the nonperturbative models at low scale with the applications of DGLAP equations with GLR-MQ-ZRS corrections. DGLAP equations with nonlinear terms is a simple tool to bridge the physics between the nonperturbative region and the perturbative region.
The last but not least conclusion we want to draw is that the partons still exist at extremely low Q 2 , although the definition/meaning of the parton distribution at low scale is not clear. The physics of partons at low Q 2 is affected by the parton-hadron duality, which still needs a lot of investigations on both experimental and theoretical sides. Based on this work, the valence quarks are the dominant partons at low Q 2 and go down fast at the beginning of QCD evolution (Fig.  12). The dynamical sea quark and dynamical gluon distributions at low Q 2 and small x are Regge-like, which have the flat forms over x (Fig. 13). The dynamical partons grow fast at small Q 2 in the evolution. The dynamical gluon distribution grows linearly with the increase of Q 2 instead of ln(Q 2 ) at low Q 2 1 GeV 2 (Fig. 13).