Global study of nuclear modifications on parton distribution functions

A global analysis of nuclear medium modifications of parton distributions is presented using deeply inelastic scattering data of various nuclear targets. Two obtained data sets are provided for quark and gluon nuclear modification factors, referred as nIMParton16. One is from the global fit only to the experimental data of isospin-scalar nuclei (Set A), and the other is from the fit to all the measured nuclear data (Set B). The scale-dependence is described by DGLAP equations with nonlinear corrections in this work. The Fermi motion and off-shell effect, nucleon swelling, and parton-parton recombination are taken into account together for modeling the complicated $x$-dependence of nuclear modification. The nuclear gluon shadowing in this paper is dynamically generated by the QCD evolution of parton splitting and recombination processes with zero gluon density at the input scale. Sophisticated nuclear dependence of nuclear medium effects is studied with only two free parameters. With the obtained free parameters from the global analysis, the nuclear modifications of parton distribution functions of unmeasured nuclei can be predicted in our model. Nuclear modification of deuteron is also predicted and shown with recent measurement at JLab.


I. INTRODUCTION
Parton distribution functions (PDFs) of nuclei are different from PDFs of free nucleon, which clearly indicates that the quark or gluon freedom inside the bound nucleon is influenced by the nuclear medium environment. The nuclear medium effect on quark or gluon distribution attracts a lot of interests on both the experimental and theoretical sides, particularly since the discovery of the EMC effect in the valence-dominant region [1]. The nuclear medium modifications on PDFs are usually depicted by the nuclear modification factor defined as the per-nucleon structure function ratio R = F A 2 /F D 2 , since deuteron is approximately viewed as a system of free proton and free neutron. The nuclear modifications are commonly categorized as the shadowing effect, anti-shadowing effect, the EMC effect, and Fermi motion smearing according to the different shapes of the ratio R in different x regions [2][3][4][5]. The shadowing effect refers to R < 1 in x 0.1 region; The antishadowing effect refers to a small enhancement of R > 1 in range of 0.1 x 0.3; The EMC effect refers to the slope of R in the valence-dominant region 0.3 x 0.7; The Fermi motion refers to the rising of R in the range of x 0.7. The xdependence and nuclear dependence of nuclear modification * email: rwang@impcas.ac.cn † email: xchen@impcas.ac.cn are so complicated that there is no commonly accepted theory which explains the nuclear corrections in whole x range for all nuclei so far. It is commonly acknowledged that different mechanisms are responsible for the different nuclear modifications in different regions of x.
Nuclear medium modification is a fundamental issue for high energy nuclear physics, because the accurate nuclear PDFs (nPDFs) of various nuclei are indispensable for simulations/calculations of relativistic nucleus-nucleus reactions on modern ion colliders -RHIC at BNL and LHC at CERN [6][7][8].
On the one hand, nPDFs provide the pre-collision condition of colliding nuclei, which plays an important role in quantifying the microscopic processes of the evolution of QCD matter during relativistic heavy-ion collision. It is also vital in the study of vector boson (like J/Ψ and Υ), jets, Z-boson, or other hadron production in nucleus-nucleus collisions. On the other hand, the precise and detailed nuclear corrections allow us to combine data across different nuclear targets and provide maximum information for the proton PDFs, such as the difference between up and down sea quark distributions. Massive targets are usually used to get sufficient statistics for neutrinonucleus deeply inelastic scattering (DIS) measurements, as the neutrino cross section is very small. The neutrino-nucleus DIS data is important for the separation of different flavor components of PDFs. Therefore, precise nuclear modifications are required if we want to include the neutrino DIS data into the global analysis of proton PDFs.
There are a lot of global analyses of nPDFs worked out so far, such as EKS [9], EPS [10,11], nCTEQ [12][13][14], nDS [15], DSSZ [16], HKN [17,18], KP [19][20][21][22], and KT [23], showing some remarkable progresses. All of these analyses were performed with the application of DGLAP equations, and the collinear factorization also works good for the nuclear PDFs case. The main difference among the widely used nPDFs is the technique of parameterizing the nuclear modifications at the initial scale of Q 2 0 > 1 GeV 2 . For EKS and EPS [9][10][11], nuclear modification factors of valence quark (R V ), sea quark (R S ) and gluon distributions (R G ) are parameterized at the input scale, under the assumptions of R uv = R dv = R V and Rū = Rd = Rs = R S . For the early analysis of EKS98, R G = R F 2 is assumed in order to get a stable evolution of nPDFs [9]. For nCTEQ analyses, the process is straight forward with the nuclear initial parton distributions of u v , d v , g,ū +d, s ands parameterized instead of the nuclear modification factors [12][13][14]. The technique of nCTEQ analysis is in close analogy to the proton global analysis. For nDS global analysis, the nuclear modification is parameterized by using a convolution method [15] with three different convolution weights for valence quark, sea quark and gluon distributions respectively. The nuclear modification parameterization of DSSZ is similar to that of EKS except that the modification factors R V , R S and R G are not independent [16]. For HKN analyses, the modification factors of u v , d v ,q and g are used. It is shown in works of HKN that the nuclear gluon distributions cannot be fixed by the present data [17,18]. KP analyses are very successful for description of the nuclear structure functions in the entire kinematical region of x and Q 2 , by developing a model which takes into account of nuclear shadowing, Fermi motion and binding, nuclear pion excess and off-shell correction [19,20]. The semimicroscopic model by KP also gives nPDFs information using impulse approximation with corrections of coherent nuclear interactions and of nuclear meson-exchange currents [21,22], which has a successful application in proton-lead collisions [24]. KT analysis of nPDFs [23] is the global analysis performed at next-to-nextto-leading order for the first time.
Gluon distribution, nuclear dependence and spatial (impact parameter) dependence of nuclear effect are three main problems needs to be address for further improving the global analysis of nPDFs [17,18,25]. Due to the challenge of not enough experimental data, limited kinematic coverage of the data points and many free parameters (> 15), gluon distribution can not be fixed and it has large uncertainty in the whole x region [17,18]. Valence quark distribution at small x and anti-quark distribution at large x in nuclei are also not clear. Accurate nuclear gluon distribution can be obtained if reli-able constraints of PDFs are applied. The dynamical parton model has a good way to constrain the gluon distribution. The gluon distribution is zero at extremely low input scale [26,27], where gluon distribution at high Q 2 is dynamically generated through the parton splitting and parton-parton recombination processes during QCD evolution. Our previous works [26,27] show that the dynamical parton model agrees well with the experimental data. The strength of parton-parton recombination or the correlation lengthR for two-parton distribution is already fixed in the global analysis of proton PDFs [26]. The only question left is the recombination enhancement due to the recombination of partons of small x between two different nucleons in a nucleus. On the process of two-parton recombination, the enhancement is proportional to the size of nucleus (A 1/3 ). To constrain nuclear gluon distribution under the dynamical parton model is the main purpose of this work.
PDFs are of nonperturbative origin. Our previous work [26] has obtained a nonperturbative input with only three valence quarks at Q 2 0 = 0.0671 GeV 2 . This provides us a good opportunity to apply the nonperturbative nuclear effects (nucleon swelling, binding effect and Fermi motion) on the nonperturbative input. Therefore, two questions are encountered for us to get nPDFs. One is how the nuclear effects modify the initial nPDFs, and the other is how the nuclear medium effects depend on the nuclear targets. This analysis is based on some models with clear mechanisms, which gives a way to know how the nuclear effects modify the nonperturbative input. In this work, the x-dependence of nuclear modification factor in whole x range is attributed to the composited influence of parton-parton recombination corrections, nucleon swelling, binding effect and Fermi motion. The last decade has seen a lot of progresses on the nuclear dependence of the EMC effect [28,29]. Local nuclear density is thought to be the most relevant quantity to explain the magnitude of the EMC effect. A novel nuclear dependence of the EMC effect [29] is used in this paper.
This article is organized as follows. The charged leptonnucleus DIS data we used are listed in Sec. II. The convolution method of Fermi smearing and the binding effect are shown in Sec. III. The nucleon swelling effect on the standard deviation of valence quark distributions is discussed in Sec. IV. DGLAP equations with parton-parton recombination corrections are shown in Sec. V, which are used to describe the Q 2 -dependence of nPDFs. In Sec. VI, procedure of the global analysis is presented. In Sec. VII, we show the global fit results, comparisons of the obtained nPDFs with experimental data and with other widely used nPDFs. In Sec. VIII, a C++ program is introduced for obtaining the nuclear modification factors of many nuclei predicted by this work. At last, some discussions and a brief summary are given in Sec. IX.

II. EXPERIMENTAL DATA
The charged lepton-nucleus deeply inelastic scattering has been a powerful tool to study the nuclear structure and the nucleon structure for decades. The high energy lepton probe is so clean that we only include the DIS data in the analysis as a baseline for further studies. The nuclear modification data of per-nucleon structure functions or differential cross sections used in this work are taken from EMC [30][31][32][33], NMC [34][35][36][37][38], SLAC [39], BCDMS [40,41], Fermilab E665 [42,43], and JLab [44] experiments. Some kinematic cuts on the experimental data are used to make sure the data are in the deep inelastic region, which is shown in Eq. (1). (1) After the selection by these kinematic cuts, we get 56, 431, 114, 25, 15 and 60 data points from EMC, NMC, SLAC, BCDMS, Fermilab E665 and JLab respectively. The measured nuclear targets used in the global analysis are listed in Table I. The kinematic ranges of nuclear data are shown in Fig. 1 (a) and Fig. 1 (b) for isospin-scalar data and all nuclear data respectively. Since the Q 2 of nuclear data are not high, the target mass correction is performed to the DIS data. The formula of target mass correction is taken from Refs. [45,46].

III. FERMI MOTION AND OFF-SHELLNESS
Compared to free nucleon, the nucleon in heavy nuclear target is moving with average Fermi momentum and of off-shell kinematic. For nuclear experimental data, the x variable is determined in the approximation that the nucleon is at rest. Therefore the measured nuclear structure function is the convolution of the bare nucleon structure function with the momentum distribution function of a nucleon inside the nucleus. Generally the Fermi smearing and binding effect deform the nuclear structure functions mainly at large x 0.7. In principle, nuclear parton distributions are non-zero in the range of 1 < x < A, and Fermi smearing effect still exists in the region. However, nPDFs in the tail region of x > 1 are very small, which can be neglected in most cases. Therefore, we only focus on nPDFs in range of 0 < x < 1 in this analysis.
In this work, the Fermi motion and binding effect are taken into account to interpret the nuclear modification at large x. These effects were well studied by Bodek and Ritchie [47], and Frankfurt and Strikman [48]. The deformation of the initial  parton distributions is given by the convolution method [2,[47][48][49]], where (3) with k F the average Fermi momentum, m N the nucleon mass, and η A = 1 − B A /m N , in which the B A is the nuclear binding energy per nucleon. The average nucleon Fermi momentum for the nuclei is around 200 MeV/c, and the binding energy is taken from Ref. [50], which is precisely measured.
The last important thing to evaluate the Fermi motion effect is to estimate the average Fermi momenta for different nuclei. According to nuclear Fermi gas model, the Fermi momentum is related to the nuclear density by k F = (3π 2 ρ/2) 1/3 fm −1 . Since the nuclear density is almost a constant for very heavy nuclei, the Fermi momentum is basically the same for very heavy nuclei. However in terms of light nuclei, it is hard to evaluate nuclear density, and Fermi gas model fails to give the Fermi momentum with precision. In this analysis, we composed an empirical formula to describe the nucleardependence of the Fermi momentum, which is written as where Z, N and A are proton number, neutron number and mass number respectively. This empirical formula gives zero Fermi momenta for both the free proton and the free neutron. Fermi momenta of several nuclei from 6 Li to 208 Pb were measured [51][52][53] by the quasielastic electron-nucleus scattering process. The free parameters k p F , k n F , t 1 and t 2 in Eq. (4) are determined by a fit to the experimental data. The quality of the fit is good with χ 2 /N d f = 7.77/8 = 0.97, which suggests that Eq. (4) is able to estimate the Fermi momentum with good precision. k p F , k n F , t 1 and t 2 are obtained to be 365 MeV/c, 231 MeV/c, 0.479, and 0.528 respectively. The average nucleon Fermi momentum of deuteron and 3 He are estimated to be 87 MeV/c and 134 MeV/c respectively within this approach.

IV. IN-MEDIUM NUCLEON SWELLING
The EMC effect in the valence quark dominant region is the most interesting subject about the nuclear modification for both nuclear and particle physicists. Up to date, there are so many models which can explain the principal features of the effect. One can look at the reviews [2][3][4][5] for a good overview of the progresses on this issue. The fact we know is that the valence quark distributions of bound nucleon are modified by the nuclear medium environment. Zhu and Shen [54][55][56][57] tried to calculate the nuclear medium deformation of quark distribution functions in the constituent quark model, and achieved a big success to well reproduce the experimental data with a few parameters.
The modification of valence quark distributions is speculated to be related to the in-medium nucleon swelling [49,54]. The in-medium nucleon swelling here refers to the enlargement of the confinement scale of a valence quark, which is a basic dynamic in strong interaction [58][59][60][61][62]. Based on our previous work [49,63], the deformations of initial valence quark distributions from the in-medium nucleon swelling effect can be evaluated by the Heisenberg uncertainty principle, where δ A is defined as the in-medium swelling factor, since the spatial uncertainty of valence quark scales with the radius of the nucleon. Here R is the radius of nucleon. σ(x q ) in Eq. (5) is the standard deviation of x, which is expressed as follows: The initial valence quark distributions of free proton are well determined so far [26,27]. In this analysis, the initial valence quark distributions of proton are taken from the Fit 2 result of Ref. [26]. Under the assumption that < x A q >=< x N q > and using the beta function form parametrization as Fig. 2.) The nucleon swelling factor δ A arises from the nuclear force and depends on the local nuclear environment, which is complicated to calculate. Our previous work [29] finds that the strength of the EMC effect linearly correlated with residual strong interaction energy (RSIE) per nucleon. The nuclear force surely plays an important role in the modification of confinement radius of a in-medium valence quark. Therefore, a simple assumption is that the nucleon swelling factor linearly scales with the RSIE, which is written as with α a free parameter which can be determined in the global fit to experiments of various nuclei.

V. SCALE-DEPENDENCE AND GLUON RECOMBINATION EFFECTS
Generally, the Q 2 -dependence of nuclear modification factor is very weak, which is difficult to detect. However the nuclear shadowing at small x and high Q 2 shows clear Q 2 -dependence for heavy nuclei [38]. For protons, the Q 2 -dependence of PDFs is precisely described by DGLAP equations [64][65][66]. It was found that the Q 2 -dependence of nuclear PDFs also obeys the leading twist DGLAP-evolution by Eskola, Kolhinen and Ruuskanen for the first time [67], which well interprets the experimental observations. Since then, all the global analyses on nuclear modifications by different collaborations are performed with DGLAP equations, yet with different techniques constructing the nuclear initial modifications or PDFs modified by various nuclear effects including the nuclear shadowing effect. These works imply that the collinear factorization also works for universal nPDFs.
Experimental measurements of nuclear shadowing at high Q 2 indicate that the nuclear shadowing is of partonic origin, which is different from the real photon absorption reaction. Mueller and Qiu [68] suggest that the shadowing in deeply inelastic scattering off nuclei is attributed to the gluon recombination process. The parton recombination process is the first non-trivial higher-twist correction. The nuclear shadowing implies that the nucleons inside the nucleus are not completely independent. The longitudinal size of a sea quark or gluon at small x can be larger than the size of a nucleon. A parton from one nucleon which has large spatial uncertainty could leak into a neighbor and fuse with one of partons in the neighbor nucleon. Therefore the parton fusion process is enhanced in nuclear target. The modifications to structure functions arising from the parton fusion reduce the structure functions at small x. The scenario of parton recombination explains the x-dependence of nuclear shadowing naturally. The smaller x of a parton, the bigger spatial uncertainty the parton has, which consequence in the stronger shadowing. Works of Qiu et al. [69][70][71][72][73] proved that the gluon recombination process is the main source of nuclear shadowing, which is strong. Our previous work [49] showed that only parton recombination process gives a good description of the nuclear experimental data. In this work, gluon fusion process is taken to characterize the nuclear shadowing.
DGLAP equations with GLR-MQ-ZRS [68,[74][75][76][77] corrections are used in this work. The gluon recombination process in GLR-MQ-ZRS corrections slows down the gluon splitting process, resulting in lower dynamical gluon distribution in nuclei. Therefore the dynamical sea quarks from gluon splitting are reduced, which leads to the observed nuclear shadowing. The DGLAP equations with gluon fusion corrections we used are 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, andR is the correlation length of the two interacting partons. Note that the integral terms as 1/2 x in above equations should be removed when x is larger than 1/2. Σ in Eq. (10) is defined as The splitting functions of the linear terms are given by DGLAP equations, and the recombination functions of the nonlinear terms are taken from Refs. [26,49,76]. The value ofR = 3.98 GeV −1 is taken as the Fit 2 result in Ref. [26]. A e f f in Eq. (9) and (10) is the gluon recombination enhancement coefficient for gluon recombination involving two nucleons in nuclear target.
Since A e f f linearly scales with the nuclear size, a parameterized formula is composed as in this analysis. (A 1/3 − 1) is the number of the shadowed nucleons. β is a free parameter, which can be determined in the global fit to nuclear data.

VI. QCD ANALYSIS
The initial valence distribution functions of free proton, initial scale Q 0 and two-gluon correlation lengthR are taken as the Fit 2 result of our previous work [26]. In this analysis, the input parton distributions at Q 0 of both free nucleon and nuclei are parameterized as the beta function form as Ax B (1− x) C . To get the initial valence quark distributions of nuclei, it is simply just to perform some modifications by Fermi motion smearing and off-shell effect, and in-medium nucleon swelling effect as discussed in Sec. III and IV. Since the nonperturbative input of free nucleon is already determined, the nuclear nonperturbative input is also determined if the average Fermi momentum, binding energy and swelling parameter δ A are all known. The initial valence quark distributions of nuclei are evolved to high Q 2 by DGLAP equations with gluon fusion process discussed in Sec. V.
The running coupling constant α s and the quark masses are the same as that in our previous work [26]. The fixed flavor number scheme is used to deal with heavy quarks in this work. The suppression of dynamical strange quark is implemented to model the flavor-dependence. Since the nuclear modification of heavy quark is not clear in theory and their contributions are trivial at not very high Q 2 , only up, down and strange quarks are used to calculate the nuclear modification factor of structure function. In this analysis, the isospin-scalar corrected per-nucleon structure functions [39] are calculated for all nuclear targets in order to compare with the experimental data. In theory, the isospin-scalar corrected per-nucleon nuclear structure function is simply just (F p in A 2 This analysis is based on the Leading Order (LO) calculations of theory. DGLAP equations with QCD corrections of parton splitting and parton-parton recombination are at LO, since parton-parton recombination evolution kernels at Nextto-Leading Order are not available so far. The running coupling α s in DGLAP equations, and the Wilson coefficient functions are also the LO results, for the consistency of the analysis.
In this work, only two free parameters α and β are used to describe the nuclear-dependence and the x-dependence of nuclear modification. The free parameters α and β are determined by the least squares method. The χ 2 function is defined as where N e is the number of data points in experiment e, D i is a data in an experiment, T i is the predicted value from QCD evolution, and σ i is the total uncertainty combing both statistical and systematic errors. Two separated global fits are performed to extract the nuclear medium corrections. One is the global fit to the data of isospin-scalar nuclei (Fit A), and the other is the global fit to all nuclear data (Fit B). In this paper, the nuclear corrections from Fit A and Fit B are called Set A and Set B respectively. Since proton and neutron have obviously different structure functions, the isospin-scalar corrections are applied for non-isospin-scalar nuclei to remove the proton-neutron non-equality [39]. The isospin-scalar corrections strongly depend on the precise structure function information of free neutron which have big uncertainties. The neutron structure functions at different Q 2 are rarely known until now with a development of the spectator tagging technique [78]. The aim of Fit A is to eliminate the uncertainty of isospin-scalar corrections, though the number of data points is cut down. The quality of the fits is good, with both χ 2 /N d f around 1.3, which is shown in Table II. The χ 2 for each experimental data set are also calculated and shown in Table I. Our model by modifying the nonperturbative input with different nuclear effects agrees well with the experimental data. The parameters α and β are obtained and shown in Table II, which are used to describe the nuclear dependence of nucleon swelling and shadowing respectively.

VII. RESULTS
The global fit results compared with the experimental measurements are shown in Figs. 3, 4 and 5. The obtained data sets are consistent with the experimental data from light nuclei (Fig. 3) target to heavy nuclei target (Fig. 4). The most precise nuclear shadowing data were measured by NMC Collaboration using two solid target sets [37]. Fig. 5 shows the result of Fit A and Fit B compared to the precise NMC data of per-nucleon structure function ratios of one nucleus to another The global fit results of structure function ratios are compared to experimental data from EMC [30][31][32], NMC [35,36], SLAC [39], E665 [43], BCDMS [40], and JLab [44] for light nuclei targets. The Q 2 of both Set A and Set B data are 5 GeV 2 in this figure.
nucleus. It is also clearly demonstrated that the nuclear dependence and x-dependence are well reproduced in the global analysis in these figures. The simple formulas expressed as Eqs. (7) and (11) are good approximations to model the nuclear dependence of nuclear modifications. The x-dependence of nuclear shadowing is well reproduced by the parton-parton recombination corrections. Fig. 6 shows the Q 2 -dependence of the nuclear modification. The predictions by the DGLAP equations with nonlinear corrections are consistent with the NMC data of 118 Sn. The scale-dependence of structure function ratio is very weak, because the parton evolution at high Q 2 is almost the same for two different nuclei. In our approach, the only difference of evolutions for different nuclei comes from the higher-twist corrections to DGLAP equations. So far the Q 2 -dependence of nuclear shadowing around x = 0.01 is in agreement with the DGLAP equations with parton-parton recombination corrections, which suggest that the DGLAP equations are uni- The global fit results of structure function ratios are compared to experimental data from EMC [31][32][33], NMC [35], SLAC [39], E665 [42,43], and BCDMS [40,41] for heavy nuclei targets. The Q 2 of both Set A and Set B data are 5 GeV 2 in this figure.
versal for hadron systems. More experimental data at small x with wide kinematic range are needed to further check the DGLAP evolution of nPDFs. Fig. 7 shows the predicted flavor-dependence of nuclear modifications on parton distributions of this work. In this analysis, no strong flavor-dependence of nuclear modifications is shown. Nevertheless, the nuclear modification of gluon shows some difference. The shadowing of gluon is a little weaker than that of sea quarks, and its anti-shadowing is quite small. In small x region (x < 10 −3 ), the strengths of shadowing effect of both quark and anti-quark are the same. However, the shadowing magnitude of quark distribution around x = 0.01 is weaker than that of anti-quark distribution, which is due to the weak shadowing effect of valence quark distribution contributing to the shadowing of quark distribution. For the EMC effect, we find that the magnitude  [34,35,37,38]. The Q 2 of both Set A and Set B data are 5 GeV 2 in this figure. of the EMC effect of gluon is smaller than that of valence quarks, and the magnitudes of the EMC effect of sea quarks are smaller than those of gluon. In our approach, the nucleon swelling effect is applied on the nonperturbative input which merely consists of valence quarks. Therefore the gluons generated by the radiation of valence quarks have weaker EMC effect, and the sea quarks from the gluon splitting have even weaker EMC effect. In the semimicroscopic model by KP, the sea quarks also have very small EMC effect [22]. The valence quarks also have the strongest anti-shadowing effect, and, the anti-shadowing effect of sea quarks is stronger than that of gluon. The dynamical sea quarks are generated from the gluon splitting in the QCD evolution. Therefore the modification factors ofū,d, ands are the same.   from different groups in terms of the nuclear modification factors. For the structure function ratio, our result is very close to the prediction of EPS09, except the shadowing effect. The structure function shadowing of this analysis at small x is a little stronger than that of EPS09 analysis, and it is close to the result of nCTEQ15 around x 0.01. For the gluon shadowing effect, this analysis gives similar prediction as that of  [11], nCTEQ15 [14], and HKN07 [18]. The resolution scales for these nuclear modifications are all at Q 2 = 5 GeV 2 .
EPS09 and nCTEQ15 around x ∼ 10 −4 . However, for antishadowing effect of gluon distribution, the prediction of this work is special, which has very weak anti-shadowing. The predictions of HKN07, EPS09, and nCTEQ15 all have very big anti-shadowing for gluon distribution. This is an interesting window to test the nuclear dynamical gluon model of this analysis. For the shadowing effect of sea quarks, our predictions are stronger than those of HKN07 and EPS09, but they are weaker than predictions by nCTEQ15. In terms of the EMC effect for valence quark and sea quark distributions, the predictions of this analysis are close to that of EPS09. The EMC effects of 208 Pb predicted by HKN07 and nCTEQ15 are stronger than our result. More nuclear DIS data, especially in small x region (x 10 −3 ) are needed to distinguish the predictions by different nPDFs analyses.
With the obtained parameters α and β, the nuclear correction of deuteron is predicted under the framework of Eq. (7) and (11). Although deuteron is a very loosely bound nucleus, its structure function is slightly different from that of free proton and free neutron. The state-of-art measurement of the structure function ratio R d EMC = F d 2 /(F n 2 + F p 2 ) is performed at JLab [79] in the kinematic region of W > 1.4 GeV and Q 2 > 1 GeV 2 . Fig. 9 shows the comparison of the predicted nuclear correction of deuteron with the JLab data. In the EMC effect region, the prediction of the global analysis is in good agreement with experimental data. The predicted Fermi motion smearing effect is weaker than the measured data. Fermi momentum of 87 MeV for deuteron used in the calculation maybe is small. The realistic deuteron wave functions are needed, or we need a more complicated formula for Fermi smearing.

VIII. NIMPARTON PACKAGE
We provide a C++ library named nIMParton to access the obtained nuclear modification factors of various nuclei for practical applications, so as to avoid the complicated nuclear effect calculations and the slow QCD evolution with partonparton recombination corrections. The C++ package is now available from us via email, the WWW[80], or downloading by the git command [81]. Two data sets of the global analysis result, called data set A and data set B, are provided by the package, which is discussed at length in Sec. VII. The given nuclear modification factor by nIMParton is the result from a interpolation of the grid table calculated by the models.
The package consists of a C++ class named nIMParton which gives the interface to the nPDFs. The constructor function nIMParton(Z, A) is used to choose a nuclear target. nIM-Parton has two important methods getRToN(Iparton, X, Q2) and getRToD(Iparton, X, Q2), which are used to get the parton distribution ratios of a nucleus to free nucleon and deuteron respectively, and suggested to be called in users' programs. Iparton set as -3, -2, -1, 0, 1, 2, and 3 corresponds to getting parton distribution ratios Rs, Rd, Rū, R g , R u , R d , and R s respectively. Another important method of nIMParton is set-DataSet. setDataSet(1) corresponds to use the data set A, and setDataSet(2) corresponds to use the data set B. For isospinscalar nuclei, the data set A is recommended. Nuclear modifications on heavy quarks are not determined in this work. We suggest using R c = R b = R g if they are needed, since heavy quarks are mainly produced by gluons. This is a feasible solution, as this analysis shows relatively weak flavor-dependence of nuclear modifications.
The nIMParton package gives only the nuclear modification factors which are calculated as the parton distribution ratios of a bound nucleon to a free nucleon. Under the assumption that the nuclear modification factors of parton distributions are the same for both the bound proton and the bound neutron, the nPDFs can be constructed by the following formula, where i, A and Z are the flavor index, the mass number and proton number of a nucleus respectively. The proton PDFs f p i are precisely determined by many collaborations through decades of development, and the neutron PDFs can easily be given by the isospin symmetry. For modern PDFs of free proton, we can look at the work [26] and the references therein.

IX. DISCUSSIONS AND SUMMARY
Two data sets of nuclear modification factors are given by a global analysis of nuclear DIS data worldwide. Data set A is from the global fit to the data of only isospin-scalar nuclei, and data set B is from the global fit to all nuclear data. Both data sets are in excellent agreement with the measured structure function ratios. The small difference between set A and set B is on the shadowing in the small x region (x < 0.01). The predicted nuclear correction of deuteron is consistent with the state-of-art measurement at JLab (not included in the global fit). Comparisons with other nPDFs are shown, such as HKN07, EPS09, and nCTEQ15. Nuclear modification factors of parton distributions from different nPDFs analyses show many differences, especially at small x. More nuclear data at small x are needed to improve the precision of nuclear corrections at small x. Our prediction of the strength of quark shadowing is weaker than the prediction by nCTEQ15, and stronger than predictions by EPS09 and HKN07. On the shadowing of gluon distribution at small x, the prediction of this work is close to that of EPS09 and nCTEQ15. The interesting characteristic of this analysis is the gluon anti-shadowing. The gluon ratio of a heavy nucleus to deuteron is almost one in the anti-shadowing region around x = 0.1. KT analysis also gives no anti-shadowing of gluon distribution, but it also gives no noticeable gluon shadowing [23]. The obtained data sets are expected to be the good options as the nPDFs input for d-A or A-A collisions. A C++ package is provided to interface with the obtained nuclear modification factors. The nPDFs are constructed by Eq. (13).
There are a lot of free parameters to describe the xdependence and nuclear dependence of nuclear modification for the common global analyses of nPDFs. Same as the KP analysis [19,20], the global analysis are based on a combination of different nuclear models in this paper. By this method, the number of free parameters to describe the nuclear modifications is reduced to only two. Based on the dynamical parton model and models of nuclear effects, the method eliminates the arbitrary of parametrization and of the solution of the optimal fit in high dimensional space of parameters. The nuclear modification factors are more reliable, since they have more theoretical constraints.
The nuclear dependence for nuclear modifications in this work is very predictive, which will be tested by measurements of unmeasured nuclei at future EIC (Electron-Ion Col-lider) project. The Q 2 -dependence is weak. However the observed weak Q 2 -dependence at small x and low Q 2 support the description of the DGLAP equations with parton-parton recombination corrections. The flavor-dependence of nuclear modification factors in this analysis is not big. The gluon shadowing is just a little weaker than that of quarks. The anti-shadowing effect of valence quarks is stronger than that of sea quarks, and the anti-shadowing effect of sea quarks is stronger than that of gluons. The magnitude of the EMC effect of valence quarks (or quarks at large x) is stronger than that of gluon, and the magnitude of the EMC effect of gluon is stronger than that of sea quarks. This flavor-dependence is sensitive to test the nuclear models used in this analysis.
The nuclear modification of gluon distribution can be fixed in the global analysis under the dynamical parton model using combined nuclear models. The obtained nuclear gluon distribution is expected to be predictive and more reliable. The nuclear gluon distribution is vital for the vector boson, jet or other hadron productions at high energy. One window to test the nuclear modification of gluon distribution in this work is the J/Ψ production in d-A or A-A collision. We will discuss it elsewhere.