Nuclear parton distribution functions with uncertainties in a general mass variable flavor number scheme

In this article we obtain a new set of nuclear parton distribution functions (nuclear PDFs) at next-to-leading order and next-to-next-to-leading order accuracy in perturbative QCD. The common nuclear deep-inelastic scattering (DIS) data analyzed in our study are complemented by the available charged-current neutrino DIS data with nuclear targets and data from Drell-Yan cross-section measurements for several nuclear targets. In addition, the most recent DIS data from the Jefferson Lab CLAS and Hall C experiments are also added to our data sample. For these specific datasets, we consider the impact of target mass corrections and higher twist effects which are expected to be important in the region of large $x$ and intermediate-to-low $Q^2$. Our analysis is based on a publicly available open-source tool {\tt APFEL}, which has been modified to be applicable for our analysis of nuclear PDFs. Heavy quark contributions to nuclear DIS are considered within the framework of the {\tt FONLL} general-mass variable-flavor-number scheme. The most recent {\tt CT18} PDFs are used as baseline proton PDFs. The uncertainties of nuclear PDFs are determined using the standard Hessian approach. The results of our global QCD analysis are compared with existing nuclear PDF sets and with the fitted cross-sections, for which our set of nuclear PDFs provides a very good description.

A precise determination of nuclear PDFs is crucial for studies of the strong interaction in high-energy scattering processes in heavy-ion collisions, such as proton-lead (p-Pb) and lead-lead (Pb-Pb) collisions at the CERN-LHC. Furthermore, nuclear PDFs are important for highenergy neutrino interactions with heavy nuclear targets, which are sensitive to the separation of up-and downtype quarks, and hence, could provide important information for the decomposition of quark flavors in a QCD analysis [4].
Several collaborations have recently presented new determinations of nuclear PDFs using the available experimental data, improved theoretical assumptions and advanced methodological settings. For the most recent determination of nuclear PDFs, we refer the reader to the analyses by the nNNPDF Collaboration [1,2], KA15 [3], EPPS16 [5], TUJU19 [4], RKPZ [7], AT12 [10], KP14 [11], DSSZ [12], HKN07 [17] and nCTEQ15 [18]. Some of the mentioned nuclear PDF determinations are based on nuclear DIS data only. By using these data alone with a rather limited kinematic coverage, significant simplifying assumptions for the nuclear PDF parameterizations need to be taken into account. Hence, the constraints on the extracted quark and gluon nuclear PDFs are rather limited in these analyses.
Nuclear PDFs at next-to-next-to-leading order (NNLO) accuracy in pQCD have been studied for the first time by KA15 [3] in the zero-mass variableflavor-number scheme (ZM-VFNS). The more recent work by nNNPDF1.0 [1] is also performed at NNLO applying the NNPDF methodology and the resulting nuclear PDFs are determined by a Neural Network (NN) in the general-mass variable-flavor-number scheme (GM-VFNS). The most recent study by TUJU19 [4] is performed at NNLO accuracy as well, but based on the open-source xFitter package [24] within the nCTEQ framework.
It is important to perform this study considering a different framework and to define the bound nucleon PDFs relative to a free nucleon baseline using the most recent proton PDF determination. The work presented in our paper focuses on the determination of new nuclear PDF sets, which we refer to as KSASG20, at NLO and NNLO accuracy in pQCD. All available and up-to-date neutral current charged-lepton nuclear DIS, charged-current neutrino DIS experimental data, and the Drell-Yan (DY) cross-section ratios for several nuclear targets are used. The former two datasets are sensitive to the flavor composition.
Our analysis also incorporates the most recent DIS data for several nuclei at high x from the Jefferson Lab CLAS and Hall C experiments. The JLab experiments provide a wealth of nuclear DIS data in the kinematic region of large Bjorken x and intermediate to low photon virtuality Q 2 . Hence, we particularly consider the impact of target mass corrections (TMCs) and higher twist effects (HT) which are expected to be important in the kinematic range of the JLab data.
The work presented in this paper is based on the publicly available open-source APFEL package [30], which has been modified in order to accommodate the data from nuclear collisions, i.e., neutral-current charged lepton and charged-current neutrino DIS on nuclear targets. For the heavy-quark contributions, we use the FONLL-B and FONLL-C implementations of the GM-VFNS at NLO and NNLO, respectively. The standard 'Hessian' approach is used to estimate the uncertainties of nuclear PDFs for quarks and gluons due to experimental errors. The resulting uncertainties are examined in view of the sparse kinematic coverage of the available data.
For the free proton baseline, we use the most recent PDF analysis by CT18 [31], which is based mainly on the most recent measurements from the LHC and a variety of available world collider data. The CT18 PDFs are consistent with our assumptions and the kinematical cuts made for our nuclear PDF analysis. The nuclear PDFs presented in our study are available via the standard LHAPDF library in order to provide an open-source tool for phenomenological applications.
We should mention here that a large amount of new and precise data from the LHC in proton-lead and leadlead collisions became recently available [25][26][27][28][29]. These high precision data, especially the data on W and Z boson production in proton-lead collisions obtained by the ATLAS and CMS collaborations at center-of-mass energies of 5.02 TeV and 8. 16 TeV, could provide further constraints on nuclear PDFs, especially for the case of the nuclear gluon PDF. Their impact on nuclear PDFs has been extensively studied in Refs. [2,5]. In addition, an analysis of the impact of available experimental data for proton-lead collisions from Run I at the LHC on nuclear modifications of PDFs is reported in Ref. [32] where the Bayesian reweighting technique [33][34][35] was used. In Ref. [36] the impact of the single inclusive D 0 meson production data from LHCb [37] in proton-lead collisions on nuclear PDFs is quantified by the Hessian reweighting method.
It will be very interesting to repeat the analysis described here and determine a new set of nuclear PDFs by adding the data from proton-lead collisions at the LHC [25][26][27][28][29]. In terms of future work, we plan to revisit this study and consider also the hadron collider data from the LHC.
The remainder of the paper is organized as follows: In Sec. II, we briefly review the general theoretical formalism for a global QCD analysis of nuclear PDFs and our assumptions for the input parameterization. This section also describes how we include target mass corrections, higher twist effects, and heavy flavor contributions in the nuclear PDF analysis. The charged lepton-nucleus DIS, very recent nuclear DIS data from JLab experiments, the neutrino(antineutrino)-nucleus DIS data, and the Drell-Yan cross-section measurements analyzed in our study are listed and discussed in Sec. III. Then, in Sec. IV, the procedure of χ 2 minimization and the estimation of nuclear PDF uncertainties are presented. In Sec. V, we show and discuss in detail the global fit results and compare with other nuclear PDFs available in the literature. This section also includes our discussions of the fit quality and the data-theory comparison. Finally, the discussion and a summary of the main results are given in Sec. VI.

II. THEORETICAL FORMALISM AND INPUT DISTRIBUTIONS
This section presents the theoretical framework used in our analysis. First, we present the parameterization of the KSASG20 parton distributions of the nucleus. Then we discuss our method to include the heavy flavor contributions in the nuclear DIS processes. A.
The parton distributions of the nucleus In this section, we present our strategy to parameterize the KSASG20 nuclear PDFs at the input scale. Similarly to our previous analysis [3] and as in other QCD analyses available in the literature [15,17], we will work within the conventional approach which defines the nuclear PDFs, xf , for a bound nucleon in a nucleus with the atomic mass number A with respect to those for a free nucleon, xf N i (x, Q 2 0 ), through a multiplicative nuclear modification factor, W i (x, A, Z): where i is an index to distinguish the distribution functions for the valence quarks u v and d v , the sea-quarksd, u, the strange quark s ands, and the gluon g.
Using the PDFs for a bound nucleon inside a nucleus presented in Eq. (1), one can obtain the PDFs for a general nucleus (A, Z) by averaging over the number of protons and neutrons inside the nuclei. It is given by, The bound neutron PDFs, f , by assuming isospin symmetry in Eq. (1).
For the nuclear modification functions, W i (x, A, Z), we follow the QCD analyses described in Refs. [3,10,13,17,[38][39][40] and assume the following cubic-type modification function, The advantage of the cubic form with the additional term d i in contrast to a quadratic-type function, i.e. with d i = 0, is that the nuclear modification becomes flexible enough to accommodate both shadowing and antishadowing in the valence quark distributions. For a detailed investigation of these functions, we refer the reader to Refs. [3,16,17]. The explicit A-dependent pre-factor in W i in Eq. (3) is constructed in such a way that for the proton (A = 1), one recovers the underlying free proton PDFs. The parameter α is considered to be fixed at α = 1 3 . The two terms of the pre-factor 1 − A −α describe nuclear volume and surface contributions [16,41]. For the valence quark distributions, β v is fixed at 0.81, and for the sea quark and gluon densities, we fix βq = β g = 1 [16].
In general, all the coefficients a i , b i , c i and d i could carry an A-and a Z-dependence [3]. Since the experimental data do not provide sufficient information to perform a stable fit for such a general ansatz, we therefore consider the d parameter for the gluon density to be fixed at zero, i.e., d g = 0. In order to give further and enough flexibility to our modification factor, we assume the following A-dependence of the parameters a i for sea-quarks, b i and c i for all parton flavors, and d i for valence and sea quarks, The nuclear and neutrino DIS and Drell-Yan experiments analyzed in this study are performed on a variety of nuclear targets, and hence, allow us to constrain such an A dependence of our fit parameters. The coefficients a i (A) in Eq. (4) depend on the atomic number A, but not all of them are free parameters that can be fitted. Among them only aq needs to be determined from the QCD fit. There are three constraints for a uv , a dv and a g due to the sum rules for the nuclear charge Z, the baryon number A, and from momentum conservation [3,17]. The nuclear charge is given by the baryon number is expressed as and finally the momentum sum rule readŝ We emphasize that with these prescriptions the parameters a i for the valence up-and down-quark distributions are different, i.e. a uv = a dv . Hence, the nuclear correction factors for the up-and down-valence distributions are not exactly the same, but expected to be similar in shape since the other parameters in the nuclear modification functions W uv and W dv are assumed to be the same. The remaining parameters in Eq. (3) are obtained by a global χ 2 analysis.
In the KSASG20 nuclear PDF analysis, we use for the free proton PDFs the most recent PDF set of CT18 [31] at the input scale Q 2 0 = 1.69 GeV 2 , i.e., Considering the discussion above and following Eq. (2), the KSASG20 nuclear PDFs for all parton species can be explicitly written as [3,17]: As one can see from our parameterizations, we have assumed flavor dependent sea-quark densities, i.e., f , and the small differences between them come from the underlying free proton PDFs and the different number of protons and neutrons in different nuclei. The nuclear DIS data which we include in our analysis are not sensitive enough to constrain the seaquark flavor decomposition, but the neutrino DIS data and the Drell-Yan cross-section measurements are sufficiently sensitive to the separation of up-and down-type quarks [4,5,12,18]. Since we use CT18 as baseline proton PDFs, the strange quark distributions in nuclei are assumed to be flavor symmetric, f . We show that the parametrization presented above is sufficiently flexible to allow for a good fit quality to the available datasets.
In our analysis, we define the momentum fraction x with respect to the bound nucleon, x = Q 2 /[2(q · p N )], where q and p N are the photon and nucleon momenta, respectively. With this convention, x is allowed to vary in the interval 0 < x < A. However, recent studies available in the literature have shown that the nuclear structure functions fall off rapidly for x > 1 [42,43] and, hence, can be neglected [12,18]. Therefore, to facilitate comparisons to other analyses, we follow the same path as other nuclear QCD analyses available in the literature and neglect the x > 1 region.

B. Target Mass Corrections
In this section, we describe in detail how we include target mass corrections in our analysis. TMCs are formally sub-leading 1/Q 2 corrections to the leading twist structure functions, where Q 2 is the squared four-momentum transfer to the hadron. TMC effects are most pronounced in the region of large-x and moderateto-small values of Q 2 , which coincides with the region where nuclear PDFs are not very well determined from fits to the data. Since we include the neutral-current charged lepton DIS data from Jefferson Lab Hall C [44] and CLAS [45] experiments, a reliable extraction of nuclear PDFs from these high-x and low-Q 2 JLab data is expected to require an accurate consideration of TMCs in our QCD analysis.
The leading contributions to the TMC have been computed in Refs. [46,47] and have been used in several studies in the literature [6,[48][49][50][51][52][53]. The target-mass corrected structure function, F TMC 2 (x, Q), is given by where η = M 2 p /Q 2 and M p is the mass of the proton. In the above equation, ξ = 2x/(1 + τ ) refers to the Nachtmann variable [54] with τ = 1 + 4ηx 2 , which shows that the TMCs vanish in the limit M 2 p /Q 2 → 0. The superscript in F 0 2 (ξ, Q) indicates the limit, when the proton mass M p is set to zero, and the second term in Eq. (10) gives the convolution term. It is found that the magnitude of these corrections could be sizeable for lower values of the photon virtuality Q ∼ M p , and hence, one needs to take TMCs into account in the high-x and low-Q region in a QCD fit. In our analysis of the JLab data we include the leading target mass corrections due to using the Nachtmann variable and the convolution term proportional to η, but omit terms of order O(η 2 ) which are also known [46].

C. Higher Twist Corrections
The inclusion of higher twist corrections is particularly important at high x and low Q 2 . To include these effects, we use the common phenomenological x-dependent function from the study by CJ15 [52] and other studies available in the literature [49,[55][56][57]. It is given by where F (LT,A) 2 indicates the leading twist structure function. By assuming a simple A 1/3 scaling of HT effects for light nuclei, the function C HT (x, A) can be written as [58], In fact, we will need this only for carbon.

D. Heavy flavor contributions
We note that the correct treatment of heavy quark mass contributions is important for global PDF analyses. To account for the mass dependence in the KSASG20 nuclear PDF analysis for the charm and bottom PDFs, we treat heavy quarks within the GM-VFNS. We refer the reader to Refs. [59,60] for a detailed overview. For the KSASG20 analysis we use the heavy-quark schemes implemented in the public APFEL package [30]. At NLO we choose the scheme FONLL-B which implements the NLO massive scheme calculation with NLO PDFs, while at NNLO we believe that the scheme FONLL-C is the better choice since it combines the NLO massive scheme calculation with NNLO PDFs [61,62]. We refer the reader to Ref. [48] for more details of these schemes.
The CT18 proton PDFs [31], which we consider as a baseline in our study, is based on the Aivazis-Collins-Olness-Tung (SACOT-χ) heavy-quark scheme [63,64], which is widely used in the CTEQ family of PDF fits. However, in our theoretical calculations and determination of nuclear PDFs, we use the FONLL GM-VFN scheme. In order to examine the potential mismatch between these two schemes at NLO, we have calculated the inclusive DIS cross sections for a proton target using the SACOT-χ and FONLL-B mass schemes. The results showed that the differences between these two schemes on the calculated cross sections are rather small, especially for x > 0.01, i.e. the range which is covered by a large amount of nuclear data. In addition, the bulk of the data sets that we used in our study corresponds to cross-section ratios rather than to absolute cross-sections resulting in even smaller sensitivity to the choice of the scheme. Moreover, the heavy-quark production does not play a significant role in the inclusive cross sections at x > 0.2 [6].
In summary, the choice of the heavy flavor and mass scheme is not particularly critical for our analysis and we expect that it will have a negligible effect on the calculated cross sections in the kinematic regions covered by the nuclear DIS data. We refer the reader to the findings highlighted in Ref. [6] for detailed discussions.
In order to remain consistent with the CT18 [31] baseline proton PDF analysis, for both the NLO and the NNLO fits, the heavy-quark masses are fixed at m c = 1.30 GeV and m b = 4.75 GeV. The strong coupling constant is set equal to α s (M Z ) = 0.118 [65] for both the NLO and NNLO fits.

III. NUCLEAR DIS DATASETS
In this section, we discuss in detail the datasets which we have used in the KSASG20 nuclear PDF analysis.
First, we start by presenting the neutral-current charged lepton-nucleus ( ± A) DIS data. These include the bulk of the datasets in our analysis which help to extract well-constrained valence and sea quark distributions; however, they provide only a limited sensitivity to the gluon distribution and to the distinction of different quark flavors. Then we discuss the charged-current neutrino(antineutrino) DIS experimental data with nuclear targets. They depend on different combinations of quark flavors, compared with the neutral current case. The combination of neutral-and charged current data is, however, expected to be sufficiently sensitive to the flavor composition of non-isoscalar nuclei [4]. The neutralcurrent charged lepton DIS data from Jefferson Lab Hall C and the most recent CLAS data measured by JLab will also be discussed. Finally, we present the Drell-Yan dilepton pair production cross-section measurement by the Fermilab experiments E772 and E886, which are important to separate different quark flavors. This section also includes a detailed investigation of our data selection and kinematical cuts that need to be made when performing the KSASG20 global QCD analysis. A.
Neutral-current charged lepton nucleus DIS The neutral-current charged ± A DIS process is a powerful tool to study the nuclear structure and to extract nuclear PDFs. Hence, we consider the nuclear DIS data in the analysis as a baseline. The DIS of charged-leptons off nuclear targets, which initiated all studies of nuclear PDFs, provide the best constraints on nuclear modifications of the quark densities. These data are usually presented as a ratio of structure functions for two different nuclei and span the range from 0.005 to 0.95 in momentum fraction x with a maximum photon virtuality of Q 2 max = 123 GeV 2 . The nuclear DIS data at lower momentum fractions, namely x < 0.01, are sensitive to the nuclear modifications of sea quarks, Wq. The data at medium-to-large x mainly probe the valence quark densities. A separation between quarks and antiquarks is not possible with these data alone. Other available data such as Drell-Yan dilepton production and (anti)neutrino collisions off nuclear targets should be used to provide a discrimination between valence and sea quarks.
All available modern inclusive DIS measurements of neutral-current structure functions on nuclear targets are considered in our KSASG20 analysis. In particular, we use the nuclear DIS data from the NMC, EMC, and BCDMS experiments at CERN [66][67][68][69][70][71][72][73], measurements from SLAC [74][75][76][77][78], HERMES measurements at HERA [79], as well as the data from E665 experiment at the Fermilab [80,81]. The measurements of the nuclear structure functions in such experiments are typically presented as ratios of two different nuclei In Tables I and II, the measured nuclear targets used in our KSASG20 QCD analysis are listed. For each dataset, we indicate the nuclei A 1 and A 2 , which are used to construct the above structure function ratios. In addition, the experiments, the corresponding number of data points after cuts, and the published references are shown as well. In order to judge the quality of the fits, the values of χ 2 extracted from our NLO and NNLO analyses are also presented. As can be seen from Tables I and II, the number of available data points varies for different nuclei. A very large number of data points is available for the deuteron (Table I) and for heavier nuclei, such as carbon (Table II). For other nuclei, such as e.g. lithium, only a few data points are available (Table II).
In order to remain consistent with the CT18 [31] baseline proton PDF analysis, we consider a kinematical cut on the momentum transfer Q 2 : Our choice for the cut on Q 2 is the same as that of the EPS09 [15] and EPPS16 [5] nuclear PDF analyses. We do not impose any cut on the invariant square of the final state mass, W 2 , again in agreement with what was done in Refs. [5,15]. After imposing the kinematical cut on Q 2 as presented in Eq. (14), we end up with a total of N data = 1075 data points. As one can see from Table I, a large amount of these points correspond to ratios of heavy nuclei with respect to deuterium. In this work, we treat the deuteron as a nucleus. Hence, in addition to the nuclear DIS data discussed above, our analysis also includes the deuteron structure function F D 2 measurements from NMC [82], BCDMS [83,84], HERMES [85], and finally the data for the deuteronproton ratio F D 2 /F p 2 from NMC [86]. The deuteron data help to extract information on the flavor asymmetric antiquark distributions,d =ū. Therefore, these data are essential for a successful QCD fit and for extracting information on the modification factors for the deuteron. For these specific datasets, a Q 2 ≥ 2 GeV 2 cut on the momentum transfer is considered. In Table III, we list the measured deuteron structure function F D 2 and deuteronproton ratio F D 2 /F p 2 used in the KSASG20 QCD analysis. After the cuts are applied, 529 data points remain, i.e. 373 for F D 2 and 156 for F D 2 /F p 2 . Compared with other analyses, we use a large amount of deuteron data in our fit. As one can see from Table III, the inclusion of higher-order QCD corrections significantly improves the description of NMC-96 [82] and HERMES [85] data. A slight improvement is achieved for the case of BCDMS [84] at NNLO, however, we obtained rather poor χ 2 values both at NLO and NNLO. The χ 2 values for the NMC-96 [86] and the BCDMS data from Ref. [84] deserve a separate comment. Despite of relatively good χ 2 values for these specific datasets, the inclusion of higher-order QCD correction does not improve the agreement of data with theory. The variation for the χ 2 values that can be seen in the table is possibly due to the tension among some of the datasets included in our QCD analysis, and the choice of input parameterization. It is important to note that two out of five deuteron-target data sets in Table III, namely, the BCDMS F D 2 data [84] and the NMC-96 F D 2 /F p 2 data [86], have already been used in the determination of the CT18 PDFs [31] that we used as a baseline. This might result in certain double fitting of nuclear PDFs. We have checked that our KSASG20 PDFs and the CT18 PDFs lead to very similar predictions for the deuteron structure function F D 2 (x, Q 2 ) in the kinematic range of the available data which shows the internal consistency of our approach. However, a more detailed study would be required to prove that the resulting uncertainties of our KSASG20 nuclear PDFs are not underestimated by this double fitting issue.
The kinematic coverage of the world data for nuclear DIS used in the present global QCD analysis is shown in Fig. 1 for some selected experiments. The applied kinematic cut on nuclear DIS data is illustrated by the dashed line in the plot. The data points lying below the line are excluded in the present QCD analysis.

B. Nuclear DIS data from JLab experiments
JLab experiments have recently provided a wealth of nuclear DIS data in the region of large Bjorken x and intermediate to low values of photon virtuality Q 2 , which span a wide range of A. These new JLab datasets that we consider in our analysis are for C, Al, Fe and Pb and they could provide high-precision constraints for the nuclear PDFs in the region of low Q 2 and high-x (see Fig. 1). The neutral-current charged lepton DIS experimental datasets for F A 2 /F D 2 from Jefferson Lab Hall C [44] and CLAS [45] measured by JLab during the 6 GeV electron beam operation [87] are listed in Table IV. In total, the number of data points from JLab is 199, among which 96 data points come from the CLAS collaboration.
In our NLO and NNLO QCD analyses of the JLab DIS data, we take into account the effect of TMCs as explained above in Sec. II B, and HT corrections as presented in Sec. II C. The values of χ 2 for the individual data for different nuclei extracted from our NLO and NNLO analyses are presented in Table. IV. As one can see, the χ 2 values presented in this table reflect the fact that the new Hall C and CLAS data measured by JLab including the TMC and HT corrections are described well by our QCD fits. This agrees with the results of Ref. [6], where a PDF reweighting method is used to in-     Figure 1: The kinematic coverage of the world data for nuclear DIS data used in the present global QCD analysis. The kinematic cut, which we apply in our fit, is illustrated by the dashed line in the plot.  vestigate the compatibility of the available nuclear PDFs studies, EPPS16, TuJu19 and nCTEQ15, with the recently published CLAS data. Note, however, that global QCD fits, which explicitly include a particular dataset, e.g. the JLab data considered in our analysis, deliver more direct constraints on nuclear PDFs than the statistical reweighting of this data.

Nucleus Experiment Number of data points
C.

Charged-current (anti)neutrino nucleus DIS
In addition to neutral-current DIS of charged leptons, we also include the data from neutrino-nucleus chargedcurrent DIS experiments. Including the neutrino-nucleus DIS data in the analysis improves the nuclear PDF determination because it has additional sensitivity to the flavor composition of the PDFs due to the different couplings to down-and up-type quarks [4].
In the KSASG20 nuclear PDF analysis, the cross section measurements of the CDHSW ν andν Fe experiment [88] and the CHORUS ν andν Pb experiment [89] have been included.
In the single-W exchange approximation, the cross section for neutrino collisions with nuclear targets is given by where x and y are the standard kinematical variables for a DIS process, m = ν andν, and A denotes the nucleus. The coupling factor N for (anti)neutrino-nucleus collisions can be written as The cross section in Eq. (15) and for antineutrino scattering by Full expressions including higher-order corrections can be found in Ref. [48]. As can be seen from Eqs. (17) and (18), F mA 2 is proportional to a particular non-singlet combination of quark distributions. Hence, it is sensitive to both valence and sea quark densities. In addition, F mA 3 provides additional sensitivity to the flavor composition since it depends on a different linear combination of quark and antiquark PDFs. By combining the nuclear and neutrino DIS data, one can arrive at a considerably improved valence and sea quark separation in the entire region of x, where the data overlap.
In addition to the CDHSW and CHORUS neutrino DIS data, there are more neutrino scattering datasets available in the literature, namely the cross sections for an iron target measured by the NuTeV collaboration [90] and the data from the CCFR collaboration [91]. For the CCFR measurements, the quantities Q 2 and x were not publicly available for the cross sections. In addition, only the averaged structure functions F 2 and xF 3 for neutrino and anti-neutrino scattering on iron nuclei are available, which have less sensitivity to the flavor composition [4]. Hence, we do not consider the CCFR data in our analysis.
Several studies in the literature have found some unresolved tension between the NuTeV measurements and other lepton-nucleus data [92,93]. A similar tension was also reported in Refs. [14,94] when taking into account the neutrino DIS data from the CHORUS and CCFRR measurements. Detailed studies presented in Refs. [34,95] have shown that the tension with other data was specifically due to the inclusion of data from the NuTeV experiment. Due to this tension, we have excluded the NuTeV neutrino DIS data from our QCD analysis.
The νA DIS data used in our analysis are presented in Table V. The number of data points, the respective reference and the specific nuclear target are listed as well. These datasets are subject to the same standard cut as presented in Eq. (14). In total, we include 2458 data points for neutrino-nucleus collisions from CHORUS and CDHSW. The χ 2 values obtained in our NLO and NNLO fits are shown as well. Figure 2 shows the kinematic coverage of these data. Our kinematic cut is illustrated by the dashed line in the plot. The data points lying below the line are excluded from the present QCD analysis.

D. Drell-Yan cross-section data
The analysis presented in our paper also incorporates Drell-Yan dilepton pair production. The expression for the differential cross-section for the DY process is given in the literature [17,[96][97][98][99][100][101]. The proton-nucleus differential Drell-Yan cross-section can be written as a sum of two parts: the contributions of qq annihilation and the quark-gluon Compton process, Typically, the Drell-Yan cross-sections are measured in terms of the Feynman variable x F = x 1 − x 2 , where the variables x 1 and x 2 refer to the momentum fractions of the involved partons. The subprocess crosssections dσ qq(0) and dσ qq(1) for qq annihilation, and dσ qg for qg andqg scattering, combined with PDFs are given   by [17,101] dσ A qq As one can see from Eqs. (20) and (21), the Drell-Yan cross sections depend on the charged-weighted sum of all quark-antiquark flavors. This is again a different combination of quark flavors compared to the total inclusive DIS cross section. Therefore the combination of DIS and DY data may in principle help to separately constrain up and down quark distributions. We should stress here that Eqs. (20) and (21) provide the DY cross section at NLO accuracy. For the calculation of the DY cross section at NNLO, we use a NNLO/NLO K-factor [3]. In order to obtain this K-factor, we have computed the NNLO DY cross sections using the DYNNLO package [102,103] and then corrected the NLO cross sections applying this K-factor during the fit.
The datasets for the Drell-Yan cross section ratios σ A DY /σ A DY analyzed in the KSASG20 global QCD fit are listed in Table VI. The specific nuclear targets, the number of data points, the extracted χ 2 values for the NLO and NNLO QCD analysis, and the related references are listed in this table, as well. In total, we include 92 DY data points.

IV. χ 2 MINIMIZATION AND UNCERTAINTY ESTIMATION
The optimal values for the nuclear PDF parameters defined in Eq. (9) are extracted from the nuclear DIS and and Drell-Yan data using the global function χ 2 global ({ξ i }) given by where m labels the experiment. w m allows us, in principle, to include datasets with different weight factors. However, we always use the default value w m = 1 for all experimental datasets [104][105][106][107][108][109]. Each experiment contributes with χ 2 m ({ξ i }) to the global χ 2 . These terms epend on the fit parameters ({ξ i }), which are identified with the parameters of the bound proton PDFs at the initial scale. χ 2 m ({ξ i }) is calculated as Here j runs over data points, m indicates a given individual dataset, and N data m corresponds to the total number of data points in this set. In the above equation, O data j is the value of the measured data point for a given observable, and ∆ data j is the experimental error calculated from the statistical and systematic errors added in quadrature. The theoretical predictions for each data point j are represented by T theory j ({ξ i }), which has to be calculated at the same experimental kinematic point x and Q 2 using the DGLAP-evolved nuclear PDFs with given parameters ({ξ i }). We use the CERN subroutine MINUIT [106] to determine the independent fit parameters of nuclear PDFs f (A,Z) i (x, Q 2 ; A, Z) by minimizing χ 2 global ({ξ i }). In Eq. (23), ∆N m describes the overall normalization uncertainty for each charged lepton DIS experiment. We include the normalization N m of different experiments as a free parameter along with other independent fit parameters ({ξ i }). First, we determine their values in a global pre-fit, then we fix them on their best-fitted values when we determine the uncertainties of the nuclear PDF parameters.
The quality of the QCD fit can be estimated from the resulting χ 2 /N data , where N data indicates the number of data points.
After describing the method to obtain the central value of the KSASG20 nuclear PDFs by minimizing the χ 2 global ({ξ i }) function, we are in a position to present our method to estimate the uncertainties of our nuclear PDFs. There are three established methods, namely the Hessian method [107,110], the Monte Carlo (MC) method [1,111] and Lagrange multiplier method [105], which can be used for the error analysis. The analysis of the uncertainties in KSASG20 is done using the standard 'Hessian' approach [5,60,107,110], which we will briefly describe in the following. For the uncertainty estimate, we follow the notation adopted in Refs. [107,110] and refer the reader to these publications for a detailed discussion of the Hessian formalism.
The KSASG20 nuclear PDF uncertainties are estimated by using the Hessian matrix, H, defined by  24) is the tolerance for the required confidence interval. It is calculated so that the confidence level (CL) becomes the one-σ error range, i.e. 68% CL, for a given number of independent fit parameters. In an ideal case, with the standard 'parameter-fitting' criterion for one free parameter, one would choose the tolerance criterion T 2 = ∆χ 2 global = 1 for a 68%, i.e. one-sigma CL, or T 2 = 2.71 for a 90% CL [60]. In the KSASG20 nuclear PDF analysis, the tolerance for χ 2 global is based on the method presented in Refs. [5,17,60]. In our study with 18 free fit parameters it becomes ∆χ 2 global = 20 at the 68% CL.
We note that other groups, e.g., nCTEQ15, EPPS16 and TUJU19, base their results on different values of the tolerance criterion: ∆χ 2 = 35 for nCTEQ15, ∆χ 2 = 52 for EPPS16, and ∆χ 2 = 50 for TUJU19.  The Drell-Yan cross section ratios σ A DY /σ A DY measured at FNAL used in the KSASG20 nuclear PDFs global fit. The specific nuclear target, the number of analyzed data points, and the related reference are listed, as well.

V. RESULTS AND DISCUSSIONS
In the following section we present the main results and findings of our QCD analysis. We first discuss the main features of the KSASG20 nuclear PDF parameters. Then, we assess the stability of our NLO and NNLO results with respect to the perturbative order. We present a detailed comparison with the recent NLO and NNLO  nuclear PDF analyses available in the literature. Finally, the section is concluded with a discussion of the quality of our fit results by comparing the resulting structure function ratios with the nuclear DIS experimental data; we compare our theoretical predictions with the neutrino DIS and Drell-Yan data as well.

A. Best fit parameters
In this work, we analyze nuclear PDFs using the CT18 proton PDF set as a baseline [31]. The nuclear modification factors in Eq. (3) are extracted from QCD fits to the nuclear and neutrino(antineutrino) DIS, and Drell-Yan data. Our best fit parameters obtained in the KSASG20 NLO and NNLO fits at the initial scale Q 2 0 = 1.69 GeV 2 are presented in Table VII along with their errors. Values marked with an asterisk (*) in this table are fixed at the given particular value since the analyzed nuclear and neutrino DIS and Drell-Yan data could not constrain these parameters well enough. The fixed values of β g = 1 and β q = 1 for gluon and sea quarks, as well as the value β v = 0.81 for valence quark densities is motivated by the HKN07 analysis [17]. Freeing these parameters can easily lead to unphysical fit results, and hence, we have decided to keep them fixed at this stage. As we mentioned before, the heavy-quark masses are fixed at m c = 1.30 GeV and m b = 4.75 GeV to be consistent with the CT18 proton PDFs. The strong coupling constant is taken as α s (M Z ) = 0.118 [65].
As discussed in Sec. II A, the nucleus-dependent parameters a i (A, Z) for the sea quark densities need to be determined from the fit to data. The parameters a i (A, Z) for the valence quark and gluon densities depend on the mass number A and atomic number Z in general and are extracted from the three constraints in Eqs. (5), (6) and (7). The numerical values for these parameters are presented in Tables VIII and IX at NLO and NNLO accuracy, respectively.
Regarding the best fit parameters and their errors listed in Table VII, some comments are in order. The obtained parameters for the nuclear valence-quark distributions reflect the fact that the nuclear DIS data analyzed in this study constrain these distributions well enough. In addition, neutrino DIS and Drell-Yan data also play an important role in obtaining a consistent behaviour for the up and down valence quarks. As can be seen from Table VII, some fit parameters of our sea-quark and gluon densities come with larger errors, especially dq and cg . In addition, we fixed d g to zero. The nuclear DIS and Drell-Yan data only loosely constrain the gluon nuclear modifications because they cover a too limited range in Q 2 .
To further constrain the nuclear sea-quark and gluon densities and reduce their uncertainties, other observables will have to be taken into account. Once more data are included, for example the data from the LHC and a future eA collider such as LHeC or FCC-he, it should be possible to relax some of the assumptions mentioned above.

B.
KSASG20 nuclear PDFs and their uncertainties In the following, we discuss the KSASG20 nuclear PDFs including the nuclear modification functions and present a detailed comparison between our NLO and NNLO analyses.
In Fig. 3, we show representations of different types of nuclear modification functions for some selected nuclei, deuterium (D), beryllium (Be), iron (Fe) and gold (Au) at the scale Q 2 = 2 GeV 2 . The nuclear modification functions are shown for the valence-quark W uv and W dv , sea-quark Wq and gluon W g at NLO (top row) and NNLO (bottom row). We repeat here that, in this work, we have treated deuterium as a nucleus in the fitting procedure. Hence, as one can see from Fig. 3, small deviations from the baseline CT18 proton PDFs are found for deuterium. The deviations from the CT18 PDFs become larger with increasing atomic mass, and significant effects are found for heavier nuclei, such as gold. We should notice here that our results in the small-x region, i.e. x < 10 −2 , are not directly constrained by the nuclear and neutrino DIS and Drell-Yan data, but determined by extrapolation based on our parametrization.
As can be seen from Fig. 3, the typical nuclear modification effects, such as anti-shadowing, shadowing and EMC suppression, are visible at the Q 2 = 2 GeV 2 for the up and down valence quark modifications, W uv and W dv . For the NLO analysis, the nuclear modification for the gluon shows a rapid rise with increasing x, x > 0.1. This trend repeats itself for the NNLO analysis, a behaviour which is similar to what one can observe in the analyses by HKM01 [16], HKN07 [17] and KA15 [3]. This behaviour may be an artefact due to the used parametrization in these analyses. Figure 3 also shows that for the case of sea-quarks one can observe the typical nuclear modifications. However, the magnitude of these effects slightly differs at different perturbative orders. Our nuclear modifications for quark and gluon densities are flat in the small-x region. This behaviour is similar to the analysis of HKN07 [17].
We continue with the discussion of the nuclear modification factors and their uncertainties for the case of lead (Pb) as an example of a large nucleus. Lead is particularly relevant for the present and future heavy-ion program at the LHC for p-Pb and Pb-Pb collisions. The nuclear modification factors for lead along with their uncertainties at 68% CL with ∆χ 2 = 20 are shown in Fig. 4. The results are shown at the scale Q 2 = 2 GeV 2 at NLO and NNLO accuracy. The uncertainty bands at very small and large x are not directly constrained by data. They are affected by the restricted flexibility of the considered parametrization and the limitations of the fitting framework. Due to the limited sensitivity of the available nuclear data to the sea-quark and gluon densities,  one has to limit the number of shape parameters. One can expect that additional new data will provide better constraints and one can envisage to consider more fit parameters allowing a larger flexibility of the modification functions.
In the following we discuss the effect of TMC and HT corrections on the fit quality and the shape of the extracted nuclear PDFs. At NLO, we observe that the inclusion of TMCs slightly increases χ 2 /N dof from 1.02 to 1.06. This can be traced back to the JLab Hall C and CLAS data. Including the HT corrections of Sec. II C yields a small improvement: χ 2 /N dof is reduced from 1.06 to 1.05 and the TMC and TMC+HT fits are quite similar. The same behaviour is seen in the case of the NNLO fit. We conclude that the inclusion of TM and HT corrections is not crucial for a good description of the presently available high-x data [6].
The effect of the inclusion of TMC and HT on the extracted nuclear PDFs is apparent from Fig. 5, where we compare our fit without TMC and HT (KSASG20 without TMC+HT) with our final fit, which includes both of these corrections (KSASG20). The first observation emerging from the comparisons presented in this plot is that while the TMC and HT affect the large x region, they do not significantly change the results at small x. As one can see, the KSASG20 valence quark PDF is slightly larger than (KSASG20 without TMC+HT) and this effect is mostly localized at large values of x as expected. Also the gluon and the charm-quark PDFs are modified by TMC+HT effects at large, while somewhat suppressed at intermediate values of x.
The resulting nuclear PDFs are presented in Fig. 6 for iron (left) and lead (right) at Q 2 = 10 GeV 2 to show the effects of the Q 2 evolution. As we mentioned in Sec. II A, the up and down nuclear PDFs have been assumed to be flavor dependent, i.e., xd = xū. For the strange quark x Figure 5: The ratio of KSASG20 nuclear PDFs at NLO which include TMC and HT effects compared to the corresponding results without such corrections for lead at Q 2 = 10 GeV 2 . We also show the uncertainty bands computed with the Hessian method.
distributions in the nuclei, we assume as usual xs = xs. The perturbatively generated heavy quark densities, xc and xb, are obtained through DGLAP evolution. As can be seen in the figure, all the gluon and sea-quark densities come with relatively large error bands at small x, reflecting the fact that there are not enough data constraints below x ∼ 0.01. We find only very small differences for the fittedū andd nuclear PDFs; the corresponding error bands shown in Fig. 6 are difficult to distinguish. The small differences betweenū andd nuclear PDFs mainly come from the underlying free proton PDFs and the different number of protons and neutrons in different nuclei. From our definition, theū andd PDFs are equal for isoscalar nuclei such as calcium and carbon.
In the following, we compare our NLO and NNLO results. In Fig. 8, the NLO and NNLO nuclear PDFs are compared at Q 2 = 10 GeV 2 for the valence quark xq v , gluon xg, strange quark xs, sea quarks xū and xd, and finally the perturbatively generated charm quark xc density. It is worth noticing here that the magnitude and the shape of the nuclear PDFs for a given flavor at some arbitrary scale Q 2 depends on the chosen set of reference PDFs for the free proton. Due to the limitations of the applied fitting framework and the limited sensitivity of the analyzed nuclear and neutrino(antineutrino) DIS, and Drell-Yan data to the gluon and sea-quark nuclear PDFs, the provided uncertainty bands are rather large, especially for small values of x.
A few remarks concerning the comparison between our NLO and NNLO analyses are in order. For both lead and gold nuclei, the valence quark, xu v and xd v , and strange quark xs PDF densities at NLO and NNLO accuracy are very similar in size. The sea quark densities, xū and xd, are slightly different at NLO and NNLO accuracy in the region of small x, x < 0.01. A significant difference can be found for the gluon xg and the perturbatively generated charm quark xc density at NLO and NNLO accuracy. As can be seen from Fig. 8, the NNLO gluon and charm quark PDFs are smaller than at NLO at small values of x.
To quantify the magnitude of NNLO corrections, we also present ratios of nuclear PDFs obtained in the NNLO fit over those of the NLO fit. This comparison is shown in Fig. 9 for lead. The results are displayed at the scale Q 2 = 10 GeV 2 and we include the one-σ uncertainty bands for ∆χ 2 = 20. As can be seen, the uncertainty for the nuclear gluon density slightly decreases when going from NLO to NNLO accuracy due to the improved overall fit quality when higher-order QCD calculations are taken into account. However, the differences between NLO and NNLO nuclear PDFs are rather small for all other parton species. These findings are consistent with the perturbative convergence of the global χ 2 discussed in Sec. IV and listed in Tables I, II, III, V, IV and VI. Concerning the fit quality of the total nuclear and neutrino DIS, and the Drell-Yan datasets, the most noticeable feature is a small improvement upon inclusion of higher-order QCD corrections. The inclusion of NNLO QCD corrections affects the nuclear PDFs uncertainty and improves the description of the data as well.

C.
Comparison with other nuclear PDF sets In this section, we present a comparison with the most recent nuclear PDF determinations available in the literature, namely nCTEQ15 [18], EPPS16 [5] and TUJU19 [4]. Since the nCTEQ15 and EPPS16 analyses were performed only at NLO accuracy, we limit the comparison to this perturbative order. All the comparisons presented in this section have been generated by using the standard LHAPDF6 library [117] and the published grids.
Each of these nuclear PDF analyses is based on a set of assumptions, for example, the form of the input parameterization at the initial scale, the choice of the proton baseline PDFs, the included datasets and the kinematical cuts applied to the data, the perturbative order, and the scheme for the heavy quark contributions.
In order to compare our results with other nuclear PDF sets, we begin with the detailed comparisons of nuclear modifications in lead. In Fig. 10, we compare the KSASG20 nuclear modification factors in lead at NLO accuracy to those of nCTEQ15 [18], TUJU19 [4] and EPPS16 the 68% uncertainty estimation with ∆χ 2 = 20 obtained using the Hessian method. One should keep in mind that the nCTEQ15 [18] analysis is based on the tolerance criterion ∆χ 2 = 35, while in the analysis by TUJU19 [4] the condition ∆χ 2 = 50 is used and EPPS16 [5] have performed their error calculation for ∆χ 2 = 52. As one can see in Fig. 10, in the case of valence up-quarks, we find a behavior similar to nCTEQ15 over the small and intermediate values of x, but a stronger large-x suppression.
In the case of down-valence and gluon modifications, differences both in shape and uncertainty bands can be seen between KSASG20 and nCTEQ15. The obtained valence modifications for KSASG20 are very similar both in shape and error bands. For the case of strange-quark modification, the result is compatible with nCTEQ15 within the estimated uncertainties at medium values of x. A comparison of KSASG20 with TUJU19 is also presented in Fig. 10. The up-valence and sea-quark modifications for the two fits can be considered compatible since the error bands always overlap over the whole range of x. For the sea quark, the TUJU19 uncertainties appear clearly larger than those of KSASG20. As can be seen, the central values for the gluon and down-valence are rather different. For both distributions, the TUJU19 uncertainties are clearly larger than those of KSASG20, except in the small-x region for the gluon modification. Our smaller uncertainties may be due to the fact that our parametrization is less flexible, especially in the case of the gluon density. In comparison with EPPS16, we find several differences and similarities. As one can see from Fig. 10, the uncertainty for the sea-quark density for EPPS16 is much larger than the one of KSASG20, and also larger than the results of other groups. For the gluon density we find compatible results at intermediate and small values of x, but KSASG20 and EPPS16 are different both in shape and central value at high-x. For the case of valence-quark nuclear modifications, one can again see that the central values are compatible only for x < 0.1. The valencequark uncertainty bands for KSASG20 are tighter than for EPPS16.
In spite of some differences for the nuclear modification factors, when calculating the total nuclear PDFs (see the next section), we find good agreement with the other analyses available in the literature.
For the full nuclear PDFs, we begin with a detailed comparison with the most recent nuclear PDF determination by TUJU19. Regarding the experimental DIS data to determine the nuclear PDFs, both KSASG20 and TUJU19 are based on the same datasets with different kinematical cuts. However, in addition to the deuteron structure function F D 2 from NMC [82], we also enrich our analysis with the F D 2 data from BCDMS [83,84] and HER-MES [85], and data for the ratio F D 2 /F p 2 from NMC [86]. Our KSASG20 analysis also incorporates the data from Drell-Yan cross-section ratios for several nuclear targets [112,113] and the most recent DIS data from the Jefferson Lab CLAS and Hall C experiments [44,45]. For the JLab data, target mass corrections and higher twist effects are taken into account in KSASG20. The TUJU19 nuclear PDF sets are based on a CTEQ proton baseline fitted within the same framework. TUJU19 also assumed flavor symmetric sea quark densities, i.e.ū =d = s =s.
The uncertainties for both analyses are obtained using the Hessian method, and TUJU19 calculated the uncertainty for ∆χ 2 = 50.
The comparison with TUJU19 is presented in Fig. 11 at Q 2 = 100 GeV 2 for lead at NNLO accuracy. Concerning the shapes of these nuclear PDFs, a number of interesting differences between the two sets can be seen from the comparison presented in this figure. Small disagreements are found for the valence and sea-quark densities, however, the two sets still agree at the one-σ level. A moderate difference is observed for the strange quark density below x < 0.1. A more pronounced difference in shape is observed for the gluon, charm and bottom quark PDFs, for which the KSASG20 distributions are more suppressed at medium to small values of x. The differences in shape among these three densities are more marked in the case of the gluon density and bottom quark PDFs. One should remember that neither of these analyses includes data which are directly sensitive to the gluon distribu-tion. The origin of the differences between KSASG20 and TUJU19, at medium to low x for the gluon and sea-quark densities, and for medium x in the case of valence quark PDFs, is likely to be due to the input parameterization and a larger number of data points for the deuteron included in the KSASG20 analysis. Another origin of these differences is due to the inclusion of Drell-Yan data and the most recent nuclear DIS data from JLab in our analysis.
In Figs. 12 and 13, our nuclear PDFs at the scale Q 2 = 100 GeV 2 are presented for lead at NLO accuracy. The most recent NLO nuclear PDF determinations available in the literature, namely from nCTEQ15 [18], EPPS16 [5] and TUJU19 [4] are also shown for comparison. We should mention here that the nCTEQ15 analysis is based on the tolerance criterion ∆χ 2 = 35, while EPPS16 presented their results for ∆χ 2 = 52. The uncertainty bands for the KSASG20 nuclear PDFs are obtained using the Hessian method with ∆χ 2 = 20, and are related to the parameters in nuclear modification factors.
In the EPPS16 analysis, the bound nucleon PDFs are defined relative to the free nucleon baseline CT14 PDFs [118], as for the case of our previous study KA15 [3]  where we considered the same framework and used JR PDFs [119]. The EPPS16 analysis was the first study which used data from the LHC for Z and W ± boson [26][27][28] and dijet production [25] in proton-lead collisions. These collider data provide further constraints for the gluon nuclear modifications and for the flavor separation.
In the nCTEQ15 analysis, the nuclear PDFs are parameterized by a polynomial functional form, in which all the A dependence is encoded in the coefficients of the parameterization. nCTEQ15 assumed s =s and the strange quark density is related toū +d by an additional A-dependent factor, s =s = (κ(A)/2)(ū +d). The TUJU19 analysis considered flavor symmetry for the sea quark densities, u =d = s =s, while in EPPS16 and KSASG20 only s =s was used as a constraint.
After presenting the main properties of these recent nuclear PDFs, we compare their results at NLO accu-racy. As can be seen from Fig. 12, for the xu v and xd v densities, the KSASG20 results are in agreement in size with TUJU19, nCTEQ15 and EPPS16, and well within their uncertainties, despite differences in the dataset and parametrization, etc. For the valence quark densities we find that both KSASG20 xu v and xd v slightly tend to stay below all other results at medium values of x ∼ 0.1. Moderate differences for the gluon density and significant differences for the sea-quark densities are observed. For both cases, EPPS16 exhibits relatively wider error bands compared with other analyses. The gluon distributions from nCTEQ15 and EPPS16 fall below the one of KSASG20.
In Fig. 13 we show a comparison for the sea-quark distributions (upper row) and the charm and the bottom distributions (lower row). The latter two are perturbatively generated. Moderate differences between KSASG20, nCTEQ15 TUJU19, and EPPS16 can be seen. We find The nuclear modification factors KSASG20 in lead at NLO accuracy compared to the nuclear PDF sets nCTEQ15 [18], TUJU19 [4] and EPPS16 [5] shown at the scale Q 2 = 100 GeV 2 . The comparison is presented per parton flavor i.
The bands for the KSASG20 show the 68% uncertainty estimation with ∆χ 2 = 20 obtained using the Hessian method.
that the agreement between the results of KSASG20 and EPPS16 for the central values of xū, xd and xc for the whole range of x is slightly better than with the other two groups. The resulting uncertainties for KSASG20 are somewhat smaller. We should stress again that all results presented here are based on the choice ∆χ 2 = 20 for the tolerance. Choosing the larger tolerance value ∆χ 2 = 50, as preferred by other groups, the error bands of our nuclear PDFs would increase by a factor of ∼ 1.5.

D. Fit quality and comparison of data and theory
In Tables I, II, III, V, IV, VI and VII presented above in section III we have shown the values of χ 2 per data point for each individual dataset, both for the NLO and the NNLO fits.
In total, we find for the KSASG20 fit at NLO, χ 2 = 4582 with 4353 data points. With 18 free parameters this leads to a χ 2 /d.o.f = 1.05 which indicates a relatively good fit. At NNLO, the KSASG20 fit leads to χ 2 /d.o.f = 1.04. This is only a moderate improvement upon inclusion of the higher-order QCD corrections.
While most datasets for nuclear and neutrino DIS experiments and Drell-Yan data satisfy the goodness of fit criterion, there are some experiments which stand out as having a poor fit. We also notice that for some datasets, χ 2 is poor even for the NNLO fit. In addition, for some individual datasets χ 2 increases as higher-order QCD corrections are included. A similar observation was made in previous analyses by TUJU19 [4] and nNNPDF1.0 [1].
As one can see from Table V, with the exception of νPb data from CHORUS, for all other (anti)neutrino-nucleus DIS data the inclusion of higher-order QCD corrections leads to a better fit quality.
In order to illustrate our discussion, we present selected comparisons of the datasets used in this study to the corresponding NNLO theoretical predictions obtained using the KSASG20 NNLO nuclear PDFs. In the following plots, we combine experimental data for each nucleus for a wider range of Q 2 -values (roughly between 1.8 and 67 GeV 2 ) where scaling violations are observed to cancel in the considered ratio. Our NNLO results have been calculated at Q 2 = 5 GeV 2 . In Fig. 14  parison is displayed for the nuclear DIS data for the ratio F A 2 (x, Q 2 )/F C 2 (x, Q 2 ) as a function of x. In Fig. 15 we also compare our NNLO results for F A 2 (x, Q 2 )/F D 2 (x, Q 2 ) and F A 2 (x, Q 2 )/F Li 2 (x, Q 2 ) as a function of x with some selected nuclear DIS data. The data shown in these plots are measured by the NMC and E139 Collaborations. The bands show the 68% uncertainty estimates with ∆χ 2 = 20. The comparisons presented in these plots demonstrate that the agreement between our NNLO theoretical predictions and the nuclear DIS experimental measurements varies between different data for different nuclei. Apart from a few data points in the small-x region, the agreement with most of the data published by the NMC and E139 Collaborations is excellent.
A detailed comparison with the CHORUS data on neutrino(antineutrino) lead collisions [89] and the CDHSW data on neutrino(antineutrino) iron collisions [88] are shown in Figs. 16 and 17, respectively. The results are shown as a function of Q 2 for some selected bins of x and y. The incident beam of neutrino(antineutrino) energies are not high enough to reach small values of x. Here we consider the range from x = 0.125 to x = 0.65, corresponding to the range between y = 0.3 and y = 0.5. Very good agreement is achieved for the neutrino(antineutrino)-nucleus data presented in these plots for the whole region of x and Q 2 .
The data from JLab experiment [45] included in this study provide important additional constraints on the shape and the uncertainty bands of nuclear PDFs in the high-x and low-Q 2 regime [6]. In order to show the agreement between our results using the extracted nuclear PDFs, in Fig. 18, we present the data/theory comparison for carbon (C), aluminum (Al), iron (Fe), and lead (Pb) nuclei. The bands in this plot show the 68% uncertainty estimates with ∆χ 2 = 20. We have again combined experimental data for different Q 2 -values (between 1.8 and 46 GeV 2 ) while our theory results are calculated at Q 2 = 3 GeV 2 . As one can see, our NLO and NNLO theory results nicely describe the recent JLab CLAS [45] data.
We now turn to the Drell-Yan cross section ratios σ A DY /σ A DY . In Fig. 19, we display the cross section ratio measured by the Fermilab experiment E866 for some selected nuclear targets [112]. Our NLO and NNLO results along with their uncertainties are shown as well. The Drell-Yan cross sections are presented in four bins of the invariant mass and are shown as a function of the proton momentum fraction x 1 . As one can see, except for some isolated points with relatively large errors, the KSASG20 NLO theory results describe the data well. In addition, the results show that the uncertainty bands at NNLO are slightly smaller than at NLO.

VI. SUMMARY AND CONCLUSIONS
In summary, in this work, we have introduced a new set of nuclear PDFs at NLO and NNLO accuracy in pQCD. KSASG20 nuclear PDFs are obtained from most up-to-date experimental data, including neutral-current nuclear DIS with several nuclear targets, charged-current neutrino DIS, and Drell-Yan cross-section measurements. The combination of these data are sensitive to the flavor decomposition of nuclei. Our analysis also incorporates the most recent DIS data from the Jefferson Lab CLAS and Hall C experiments. For these specific data sets, we take into account target mass corrections and higher twist effects which are mainly important in the region of large x and intermediate-to-low Q 2 . Heavy quark mass effects are included in the FONLL general-mass variableflavor-number scheme for charm and bottom quarks.
Our determination of nuclear PDFs includes error estimates obtained within the Hessian method with the tolerance criterion of ∆χ 2 = 20. The effects arising from the inclusion of higher-order QCD corrections are investigated. We found only small differences between the NLO and NNLO QCD fits, both for the shape and size of the uncertainty bands. The inclusion of higher-order QCD corrections slightly improves the description of the nuclear data analyzed in this study, but there is no strong indication that NNLO corrections are required by the data. The largest differences appear for the gluon nuclear PDF, where our NNLO fit leads to a slight decrease in the uncertainties.
TUJU19 [4] which is based on the CTEQ framework, we parameterize the nuclear PDFs considering nuclear correction factors. To this end, a modern set of parton distribution functions for free protons, namely CT18, is considered as reference. Our results are consistent, within uncertainties, with the previous determinations of nuclear PDFs available in the literature, in particular for nCTEQ15 [18], EPPS16 [5] and TUJU19, which are based on a different selection of data sets and assumptions. However, we found a number of differences which only occur in regions without any constraints from data. These differences can be attributed to different assumptions such as the input parametrization of the nuclear modification factors.
The nuclear PDFs analysis presented in this article represents the first step of a broader study. A number of improvements are foreseen for the near future. While the data used in the present study allowed us to determine the nuclear quark and anti-quark densities, the nuclear gluon PDF is only loosely constrained. This is actually the most important limitation of the KSASG20 nuclear PDFs analysis. To resolve this limitation, we plan to include additional datasets, especially from present and future measurements of proton-lead and lead-lead collisions at the CERN-LHC, which are expected to provide direct information on the nuclear gluon modifications and more stringent flavor-dependence constraints. These include, for example, the data from PHENIX and STAR for inclusive pion production in deuteron-gold (d-Au) collisions [114,115], neutral pion production data from PHENIX [114], and the charged and neutral pion data from the STAR experiment [115,116]. Also data from proton-lead (p-Pb) collisions from the ATLAS and CMS Collaborations at the LHC, which have been included in the analyses of nuclear PDFs performed in Refs. [2,5], are expected to improve the constraints of the nuclear gluon density at large momentum fractions. With more precise data, it might be necessary in the future to consider more flexible parametrizations.
We can also expect that the electron-ion collider (EIC) [120], the Large Hadron-Electron Collider (LHeC) [121] or a Future Circular Collider (FCC) [122,  123] could provide precise data for nuclear PDF analyses. It would be very interesting to repeat the analysis described here using additional observables from hadron colliders such as the LHC. The NLO and NNLO nuclear PDF sets presented in this work are available in the LHAPDF format [117] for all relevant nuclei from A = 2 to A = 208 and can be obtained from the authors upon request. x N M C -9 5 x F 2 ( C a ) / F 2 ( L i ) x F 2 ( C ) / F 2 ( L i ) Figure 15: Comparison of our NNLO results calculated at Q 2 = 5 GeV 2 for the ratios F A 2 (x, Q 2 )/F D 2 (x, Q 2 ) and F A 2 (x, Q 2 )/F Li 2 (x, Q 2 ) as a function of x with some selected nuclear DIS data from E139 and NMC-95.  Figure 17: Same as Fig. 16, but for the neutrino(antineutrino) iron DIS data from the CDHSW measurements [88].  Figure 19: Comparison of the KSASG20 NLO and NNLO results along with their uncertainties for the Drell-Yan cross section ratios σ A DY /σ A DY with data for some selected nuclear targets from the Fermilab experiment E866 [112].