Dynamical parton distributions from DGLAP equations with nonlinear corrections

Determination of proton parton distribution functions is presented under the dynamical parton model assumption by applying DGLAP equations with GLR-MQ-ZRS corrections. We provide two data sets, referred to 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 (Q2 > 2 GeV2). Furthermore, our analysis shows that the input with flavor-asymmetric sea components better reproduces the structure functions at high Q2. Generally, the parton distribution functions obtained, especially the gluon distribution function, are good options for inputs to simulations of high energy scattering processes. The analysis is performed under the fixed-flavor number scheme for nf = 3, 4, 5. Both data sets start from very low scales, around 0.07 GeV2, where the nonperturbative input is directly connected to the simple picture of the quark model. These results may shed some lights on the origin of the parton distributions observed at high Q2.


Introduction
Hadrons are complex systems consisting of quarks and gluons, which has made precise understanding of hadron structure a long and difficult path. Thanks to the collinear factorization theorem [1][2][3] in quantum chromodynamics (QCD), the calculation of high energy hadron collisions has become much more 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 or match the PDFs of the proton. PDFs in wide kinematic ranges of Q 2 and x are an important tool to give some theoretical predictions of high energy hadron collisions and simulations of expected interesting physics in modern colliders and JLab high lu-minosity experiments.
Determination of the proton PDFs has attracted a lot of interest on both the theoretical and experimental sides. To date, the most reliable and precise PDF data has come from the global QCD analysis of experimental data. There has been a lot of effort made and progress achieved on this issue [11][12][13][14][15][16][17][18]. In the global 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. The least squares regression method is usually chosen for this procedure. Finally, PDFs in a wide kinematic range are given with the obtained optimized nonperturbative input. Although a lot of progress has been made, the gluon distribution at small x is still poorly estimated, and has large uncertainties [19,20]. Even worse, the gluon distributions from different collaborations exhibit large differences. The gluon distribution needs to be more quantitative, with a number of physics issues relating to its behavior [21][22][23].
PDFs at low resolution scale are always confusing since they are in the nonperturbative QCD region. However, they are related to the nucleon structure information measured at high resolution scale. Therefore the nonperturbative input gives some valuable information about 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 do 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 sea quarks and valence-like gluons, 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 worth pointing out that the valence-like input and PDFs generated from it are positive. In some analysis of MSTW [27][28][29], negative gluon density distributions are allowed for the nonperturbative input in order to fit the small-x behavior observed at high scale. The MMHT2014 PDFs [18] supersede the MSTW2008 PDFs, making some changes to the parameterisation to stop negative distribution in some regions of x, and including LHC and 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 was developed and extended to even lower scales around Q 2 0 ∼ 0.06 GeV 2 in our previous works [30,31]. The naive nonperturbative input [30] with only three valence quarks was realized, which is the simplest input for the nucleon. In 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 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.
The DGLAP equations [32][33][34], based on the parton model and perturbative QCD theory, successfully and quantitatively interpret the Q 2 -dependence of PDFs. They have been so successful that most PDFs up until now have been extracted by using the DGLAP equations. The common way of improving the accuracy of the determined PDFs is to apply the higher order calculations of the DGLAP equations. However, there are many QCDbased evolution equations and corrections to the DGLAP equations [35][36][37][38][39] being worked out. It is worth applying 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 are expected to be more reliable at small x. The second purpose is to connect the quark model picture of the 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 the GLR-MQ-ZRS corrections by determining the value of parton correlation length R.
The organization of this paper is as follows. Section 2 lists the experimental data we used in the analysis. Section 3 discusses the QCD evolution equations, which are the most important tool to evaluate the PDFs. The nonperturbative input inspired by the quark model and other nonperturbative effects are discussed in Section 4. The other details of the QCD analysis are explained in Section 5. Section 6 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 7 introduces the IMParton package, which gives the interface of the obtained PDFs. Finally, a simple discussion and summary is given in Section 8.

Experimental data
The deep inelastic scattering (DIS) of charged leptons on nucleons has been a powerful tool to study nucleon structure for a long time. The quark structure of matter has been 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 the SLAC [41], BCDMS [42], NMC [43], E665 [44] and HERA (H1 and ZEUS) [17,45] collaborations.
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. (1) For neutral-current DIS, the contribution of the Z-boson exchange cannot be neglected at high Q 2 . Therefore we compose another kinematic cut to reduce the influence of the Z-boson exchange contribution, as shown in Eq.
With these kinematic requirements, we get 469, 353, 258, 53 and 763 data points from the SLAC, BCDMS, NMC, E665 and HERA experiments respectively. SLAC was the first to perform 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 implements some improved treatments of radiative correction and the value of R = σ L /σ T . The minus four-momentum transfer squared Q 2 of the 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 relatively 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. Precise measurements of the structure function F 2 were made by the experiments at CERN, Fermilab, and HERA at DESY. Both BCDMS and NMC data were collected from muon-proton DIS with the CERN SPS muon beam but with radically different detectors. The BCDMS data were taken at beam energies of 100, 120, 200 and 280 GeV, and the NMC data were taken at beam energies 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 a combined analysis of the 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 a run with 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 kinematics of all the data covers 3 orders of magnitude in both x and Q 2 . Since the SLAC and the NMC data are distributed at 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.

Nonlinear corrections to DGLAP equations
The DGLAP equations [32][33][34] are an important and widely used tool to describe the Q 2 dependence of quark 053103-3 and gluon densities. The equations are derived from perturbative QCD theory using the quark-parton model instead of the rigorous renormalization group equations, and offer an illuminating interpretation of the scaling violation and the picture of parton evolution with 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 have been 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 was initiated by Gribov, Levin and Ryskin (GLR) [35], and followed by Mueller and Qiu (MQ) [36], and 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 becomes so large that the quanta of partons overlap spatially. One simple criterion to estimate this saturation region is xf 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, time-ordered perturbative theory (TOPT) is used instead of the AGK cutting rules [37]. The corrections to the DGLAP equations are calculated in the leading logarithmic (Q 2 ) approximation, and extended to the whole x region, which satisfies the momentum conservation rule [38].
In this analysis, DGLAP equations with GLR-MQ-ZRS corrections are used to evaluate the PDFs of the proton. The GLR-MQ-ZRS corrections are very important to slow down the parton splitting at low Q 2 < 1 GeV 2 . To date, ZRS have derived all the recombination functions for gluon-gluon, quark-gluon and quarkquark 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 the 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 1/2 x in the above equations should removed when x is The splitting functions of the linear terms are given by the DGLAP equations, and the recombination functions of the nonlinear terms are written as [38] P gg→g (x, y) = 9 64

Quark model and nonperturbative inputs
The quark model has achieved remarkable success in explaining hadron spectroscopy and some dynamical behavior of high energy reactions with hadrons involved. The quark model uncovers the internal symmetry of hadrons. Moreover, it implies that 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 053103-4 three valence quarks. There are some model calculations of the initial valence quark distributions at some low Q 2 0 from the MIT bag model [4,5], the Nambu-Jona-Lasinio model [6] and maximum entropy [7] estimation.
Inspired by the quark model, an ideal assumption is that the proton consists of only three colored quarks at some low scale Q 0 . This assumption results in the naive non-perturbative input -the three valence quark input. At the input scale, the sea quark and gluon distributions are all zero. This idea was 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 a 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 dominant valence quarks and small quantities of valence-like components. Partons produced by QCD evolution are called dynamical partons. The input scale for the valence-like input is around 0.3 GeV 2 [11,12] and the evolution of the valence-like input is performed with the DGLAP equations. In our works, the dynamical PDF model is developed and extended to even lower 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, the three valence quark input is called input A, and the one with flavorasymmetric 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 the valence quark distribution is the time-honored canonical parametrization f (x) = Ax B (1 − x) C , which is found to depict the valence distribution at large x well. 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 the proton. Hence, we have the momentum sum rule for valence quarks in the naive input, With the above constraints, there are only three free parameters left for the parametrizations of the naive input. The naive input (Eq. (9)) is the simplest nonperturbative input for the proton, and 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 These 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 flavorasymmetric sea components accurately, the following constraint from the E866 experiment [58], shown in Eq. (14), is taken in this analysis.
Therefore, there are only 7 free parameters left for the parametrization of input B. For better discussion of the quantity of flavor-asymmetric sea, we define δ as 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 the coupling constant. In this analysis, the initial scale Q 0 is viewed as a free parameter which can be determined by experimental data.

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 have already been 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) are not considered as massless partons within the nucleon. The number of active flavors n f in the DGLAP evolution and the corresponding Wilson coefficients are fixed at n f = 3 (only u, d and s light quarks). The heavy quark flavors are produced entirely perturbatively from the initial light quarks and gluons. The FFNS predictions agree very well with the DIS data [12,24]. In this analysis, only the charm quark distribution is given, since the bottom and top distributions are trivial. The charm quark distribution comes mainly from the gluon distribution through the photongluon 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 nucleon structure studies [63]. As discussed in Section 3, the flavor-asymmetric sea componentsū AS andd AS result in thed−ū difference naturally. As found in experiments [64,65] and predicted by 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]. u 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 squares 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 point in a experiment, T i is the predicted value from QCD evolution, and σ i is the total uncertainty combining both statistical and systematic errors.

Results
Two separate fits are performed for input A, which consists of only three valence quarks. One is a fit to the whole x range (Fit 1) and the other is a fit to the data excluding the region 2 × 10 −3 < x < 0.15 (Fit 2). The results of the fits are listed in Table 1. 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 1. The obtained R values are smaller than the proton radius, which is consistent with 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 fit is bad if we use the DGLAP equations without partonparton recombination corrections, because the parton splitting process only generates very steep and high parton distributions at small x [30]. Parton-parton recombination corrections cannot 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 a similar shape to that measured in experiments, as shown in Fig. 2(a). However, it fails to depict the experimental data in detail 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 the three valence quark input miss a peak-like component in the transition region from valence-domain to sea-domain. The three valence quark input needs to be modified and developed. This discrepancy is expected to be removed by the intrinsic light quarks, 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 quark input from the analysis. This is why we performed Fit 2 to input A. Fit 2 is in excellent agreement with the experimental data at both large x and small x, as 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 over the whole x range (Fit 3). The quality of the fit is improved greatly compared to input A, as shown in Table  1. The additional flavor-asymmetric 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 flavorasymmetric 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 found to be 0.1. The obtained parameters Q 0 and R are shown in Table 2, 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 Section 4. 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 larger 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 the 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. Figures 7 and 8 show the valence quark, sea quark and gluon distributions at high Q 2 compared to other widely used parton distribution functions. The valence quark distributions exhibit some differences between our result and other recent global analyses. This discrepancy suggests that we need a more complicated parametrization for valence quark distributions beyond the simple beta function form. Sea quark distributions are consistent with each other. Our gluon distribution 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 produced purely dynamically in the QCD evolution. We argue that this gluon distribution is more reliable since no arbitrary parametrization of input gluon distribution is involved. 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 distributions at Q 2 = 2.5 GeV 2 are shown in Fig. 10 with the recent reanalysis data by the HERMES collaboration and other widely used parton distribution functions. The predicted strange quark distributions describe the experimental data well, and are consistent with the other PDFs. Our strange quark distributions are 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 much heavier than that of the up or down quark. This kind of suppression is supported by LQCD calculations. Figure 11 shows comparisons of the charm quark distributions with the measurements by the H1 and ZEUS Collaborations. The charm quark distributions are based on LO calculation of photon-gluon fusion. This method of dealing with the charm quark distributions is also used in the global analysis by GRV95 and GRV98. Although it is a simple calculation under the FFNS, the calculation of photon-gluon fusion subprocesses basically reproduces the experimental measurements of the charm quark contribution to the F 2 structure function.  In our approach, parton distribution functions at very low Q 2 are also given. We extend the input scale from Q 2 0 = 1 GeV 2 down to Q 2 0 = 0.1 GeV 2 . Our valence quark distributions at low Q 2 are shown in Fig. 12. The va-lence 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, 053103-11 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.

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 quark nonperturbative input, and data set B is from the nonperturbative input of three valence quarks adding flavor-asymmetric sea quark components, as discussed in Section 4.
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 lets the user choose data set A or data set B via set-DataSet(1) or setDataSet(2) respectively. The most important method of IMParton 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 gettingc, 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 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.

Discussions and summary
We have composed some naive nonperturbative inputs inspired by the quark model and some other nonperturbative QCD models at very low Q 2 . By using the 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 are expected 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 quark input, and the other is from three valence quarks with a few flavor-asymmetric sea components. The obtained PDFs can be justified and updated with further investigations of many other hard processes, such as the Drell-Yan process, inclusive jet production and vector meson production.
By a global analysis, we find that the quark model of 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 nonperturbative QCD effects beyond the quark model are also needed to reproduce the experimen-tal data in detail. By adding the flavor-asymmetric sea components, the quality of the global QCD fit improves significantly. This is clear evidence of the other nonperturbative parton components of the proton beyond the quark model [50][51][52][53][54][55][56][57]. It is interesting to know that the sea quarks and gluons are mainly from the parton radiations of the three valence quarks predicted by the quark model. However, there are more degrees of freedom inside the proton, which need the interpretations of 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 very steep and large parton densities because of the long evolution distance from extremely low scale. At low Q 2 , the strengths of the recombination processes are comparable to those 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 the DGLAP equations with GLR-MQ-ZRS corrections. The DGLAP equations with nonlinear terms are a simple tool to bridge the physics between the nonperturbative region and the perturbative region.
Last but not least, 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 investigation 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, and 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).