nCTEQ15HIX -- Extending nPDF Analyses into the High-$x$, Low $Q^2$ Region

We use the nCTEQ analysis framework to investigate nuclear Parton Distribution Functions (nPDFs) in the region of large x and intermediate-to-low $Q$, with special attention to recent JLab Deep Inelastic Scattering data on nuclear targets. This data lies in a region which is often excluded by $W$ and $Q$ cuts in global nPDF analyses. As we relax these cuts, we enter a new kinematic region, which requires new phenomenology. In particular, we study the impact of i) target mass corrections, ii) higher twist corrections, iii) deuteron corrections, and iv) the shape of the nuclear PDF parametrization at large-$x$ close to one. Using the above tools, we produce a new nPDF set (named nCTEQ15HIX) which yields a good description of the new JLab data in this challenging kinematic region, and displays reduced uncertainties at large $x$, in particular for up and down quark flavors.

ratio for carbon illustrating the nuclear correction factor across the various x regions. The black points indicate the data used in the original nCTEQ15 fit, and the red points with the solid squares represent the additional data from this original set which are now included due to the relaxed Q and W cuts. The red open squares are the new JLab DIS data included in this analysis, and the blue points are those JLab DIS data which are excluded by the current kinematic cuts.

I. INTRODUCTION
With the EIC and LHeC/FCC on the horizon, science is now entering a new era of precision in the investigation of hadronic structure enabled by a flood of data from JLab, RHIC and the LHC. Describing one of the four fundamental forces of nature, Quantum ChromoDynamics (QCD) -the theory of the strong interaction -remains deeply complex and enigmatic, although the Parton Distribution Function (PDF) framework has proven remarkably successful in describing processes with hadronic initial states .
While the study of proton PDFs has grown exceedingly precise, the need to extend this precision to the nuclear sector, involving fits with explicit nuclear degrees of freedom, has become more urgent in recent years in order to enhance the accuracy of experimental analyses involving nuclear targets. Progress in studying QCD dynamics within nuclei has been demonstrated across a number of recent nuclear PDF (nPDF) analyses [1][2][3][4][5][6][7][8][9][10][11]. A significant challenge in the determination of nPDFs has been the acquisition of empirical data from a sufficiently wide variety of experiments as to provide complementary constraints, and, e.g., specify the A dependence of the resulting nPDFs. For this reason, there is a continual need for new data sets to broaden global analyses. In the present work, we build upon the recent nCTEQ15 analysis by including recent JLab data covering an expanded kinematic range. As we shall demonstrate, this data has the potential to furnish an improved understanding of hadronic and nuclear structure and interactions, and, in turn, new insights into QCD.

A. JLab Kinematic Reach
The recent facility upgrades of the Continuous Electron Beam Accelerator Facility at the Thomas Jefferson National Accelerator Facility (JLab) have enabled the measurement of high precision electronnucleus scattering events in an extended kinematic regime. In particular, the JLab experiments provide a wealth of data in the relatively unexplored kinematic region of large Bjorken x and intermediate to low photon virtuality Q 2 . This mostly unexplored kinematic region is often referred to as the "transition" region from resonance dominated production to deep-inelastic scattering, and is of considerable interest to both the charged lepton and neutrino scattering communities.
In Fig. 1 we display a collection of DIS and DY data from selected experiments along with recent data from the JLab DIS experiments. We have also indicated typical cuts in Q 2 and W 2 = Q 2 (1−x)/x+M 2 N , with M N the nucleon mass, which are often implemented in many global analyses. The objective of these {W, Q} cuts is to remove those data which might have significant nonperturbative or higher twist contributions; unfortunately, these cuts exclude a large subset of the data, as indicated in Table I.
In this investigation we will relax the kinematic cuts to study whether we can describe the broader subset of the JLab data with nuclear PDFs using the nCTEQ global analysis framework [1]. As we expand this analysis into new kinematic regions, we will analyze whether we might be sensitive to new effects including i) target mass corrections, ii) higher twist corrections, iii) deuteron structure function modifications, and iv) large-x nuclear corrections to allow the nuclear Bjorken variable x A > 1. We will address each of these issues in turn, and then use the optimal combination to extract a set of nPDFs which extends into this new kinematic region.
We note here a recent study [9] where PDF reweighting methods were used to investigate compatibility of the current nuclear PDFs with the CLAS data.

B. The EMC Effect
Extending the nPDFs into this new kinematic region can also provide new information on x-dependent nuclear effects observed in many DIS experiments. Starting in the early 1980's, the European Muon Collaboration (EMC) [27] found that in the DIS kinematic region the per-nucleon structure functions F 2 for iron and deuterium were not only different, but that this difference also changed as a function of x.
This evidence for nuclear effects in DIS charged-lepton nucleus scattering can be summarized as in Fig. 2, which displays the F C 2 /F D 2 per-nucleon structure function ratio as a function of x; this behavior is consistent with measurements by both the SLAC eA [28][29][30] and the BCDMS µA [31,32] experiments. The observed behavior of the ratio R=F A 2 /F D 2 , can be divided into four x regions: • the shadowing region: R ≤ 1 for x 0.1, • the anti-shadowing region: R ≥ 1 for 0.1 x 0.3, • the EMC-effect region: R ≤ 1 for 0.3 x 0.7, • and the Fermi motion region: for x 0.7.
For a review of the data and theoretical interpretations see Ref. [33].
The current study concentrates on the higher x region so the effects of shadowing and anti-shadowing are not of direct impact for this specific study. However, the modifications at medium-to-higher-x in the so-called "EMC-effect region" will be especially important for these new JLab data sets. Since its discovery, there is still no universally accepted explanation for the EMC effect.
Recent studies have observed a correlation between the magnitude of the EMC effect and the relative amount of short-range correlated (SRC) nucleon pairs in different nuclei, suggesting that the EMC effect is driven by the modification of nucleons in SRC pairs [34][35][36][37][38]. Forthcoming experiments at JLab will explore the relationship between the offshellness of nucleons and the modification of the structure function [39,40].
Additionally, it has been demonstrated [41] that the EMC effect persists to lower values of W extending into the resonance region. This is an intriguing result that the nuclear structure functions in the resonance region exhibit similar behavior as in the DIS region. In Sec. V we will explore implications of removing the W cut on the EMC effect.

C. Outline
In Section 2 we present an overview of the various corrections we will apply to better describe the large x region. In Section 3 we introduce the new JLab data used in this study and summarize the various fits we will investigate. In Section 4 we show the results of the fits including the comparison between data and theory, and present the impact on the PDFs. In Section 5 we explore the implications of extending the nPDFs into the low W kinematic region. In Section 6 we provide a comparison of the new nCTEQ15HIX nPDFs with other results from the literature. In Section 7 we summarize our conclusions. Additionally, in Appendix A we tabulate the data sets used in the fit and provide the impact of the kinematic cuts.

II. STRUCTURE FUNCTION MODIFICATIONS
The goal of this study is to extend our global nPDF fit into new kinematic regions at large x and low Q. Many global analyses impose stringent cuts in both Q and W . For example, the nCTEQ15 fit was performed with Q>2 GeV and W >3.5 GeV. The motivation for these cuts is that the Q cut largely eliminates a variety of nonperturbative and/or powersuppressed corrections with the potential to complicate the extraction of leading-twist PDFs. As Q decreases, α s (Q 2 ) becomes large, and our perturbation expansion breaks down. Correspondingly, the W cut removes events in both the low Q 2 and large x region where contributions from non-factorizable higher-twist terms become large.
Imposing stringent cuts allows us to avoid kinematic regions that may be difficult to compute. However, the trade-off is that this can significantly reduce the data set excluding kinematic regions important for both the charged lepton and neutrino communities. In Table I we display the total number of data points which satisfy various kinematic cuts. For example, with the loosest set of cuts {Q>1.14 GeV, W >1.3 GeV}, the number of data points is 1839; this is in contrast to the very conservative cuts used in the nCTEQ15 fit, {Q>2 GeV, W >3.5 GeV}, which reduces the number to 708-less than half.
Clearly, it is advantageous to reduce the {Q, W } cuts as much as possible, but this region is challenging to compute. For example, if we simply take the nCTEQ15 fit and extend this into the low {Q, W } region without accommodating any additional phenomenological effects, the quality of the resulting theory to data comparison is no longer acceptable.
As we mentioned in the introduction, there are a variety of possible effects that may enter our analysis in this region, and we will examine them systematically to study the impact of each one to obtain the best description of the data.
In the remainder of this section, we will review the different effects, and outline how they are included in the global fit.

A. Target Mass Corrections
We begin by examining the operator product expansion (OPE) to determine what possible corrections we might encounter as we extend the kinematic region to high-x and low Q. The expressions for the structure functions can be derived in the OPE, and a detailed review is provided in Ref. [42]. Using this formalism, the leading, twist-2, structure function expressions will have corrections in powers of (M/Q) where M is the mass of the proton: 1 Here, ξ=2x/(1+r) is the Nachtmann variable with r= 1+4x 2 M 2 /Q 2 , and F 2 (x, Q) is the structure function in the limit where the proton mass M is set to zero. The additional terms, indicated by the dots, are given by convolution integrals for which explicit expressions are given in Eqs. (22)-(24) of Ref. [42]. The magnitude of these corrections can be sizeable, especially for lower scales, Q ∼ M , as we illustrate in Fig. 3 where we show the x dependence of the 2 (ξ, Q). As can be seen, for F 2 , the leading TMC gets modified by the additional TMC by up to a ∼25% at intermediate to large x values.
The leading TMC effects due to the Nachtmann variable ξ and the pre-factor x 2 /(ξ 2 r 3 ) are accounted for in all our calculations. However, we have not included the effect of the convolution terms. Since we are fitting only nuclear ratios such as F A 2 /F D 2 and σ A /σ D , the nuclear A dependence of these additive convolution terms is very minor. 2 One important question which arises is: do the TMCs scale as the mass of the nucleus M A /Q, or the mass of the nucleon M/Q where M A = AM ? This issue is especially important for the heavy nuclear targets such as lead (A=208) where a M A /Q scaling would dramatically impact the low energy data. The answer can be obtained by tracing the mass terms in the derivation of the master formula (e.g., Eqs. (22)(23)(24) of Ref. [42]) from the OPE. For a nucleon, the OPE expansion is in terms of M x/Q while for a full nucleus this becomes M A x A /Q. Here, x A is the Bjorken x for the nucleus defined as x A = Q 2 /(2P A ·Q), where P A is the full nucleus momentum and x A ∈ [0, 1]. Thus, we find M A x A /Q = (AM )(x/A)/Q ≡ M x/Q, (where x=A x A and x ∈ [0, A]) so even for heavy nuclei the TMC terms are suppressed by powers of M/Q where M ∼1 GeV. Consequently, the master formula of Ref. [42] holds for both nucleons and nuclei. There are certainly some subtle steps in extending the OPE and DGLAP formalism from the case of a proton to a heavy nuclear target, and a separate paper expanding the work of Ref. [42] with detailed derivations is in progress [45]. Figure 3 displays a band of curves for each fixed Q value which show the variation between D (the lowest curve) and Pb (the highest curve). The spread of 1 See Eq. (23) of Ref. [42] for details. We have implemented an alternative definition derived in the light-cone frame [43,44] which coincides with Eq. (1) to 1%. 2 The NMC F D 2 data set ID=5160 is the one case which is not a ratio. This set extends out to x 0.48, and we have checked that these TMC corrections are less than the data uncertainty. these bands can be as large as 5% for Q=1.0 GeV, and decreases quickly for larger Q values. For the current investigation, we will focus on cuts of Q>1.3 GeV and W >1.7 GeV. For example, from Fig. 1 we see that for Q=1.3 GeV, the W >1.7 GeV cut includes points out to x 0.46, and here the spread of this band (yellow) is on the order of ∼1%. This represents the maximum deviation between a deuteron (D) target and a lead (Pb) target; lighter nuclei such as carbon (C) would be less. Therefore, we see that within our kinematic cuts, the TMC correction provides a uniform shift, within about 1% or less, of all the nuclei. Since these additional TMC corrections shift both the A and D results by nearly the same factor, the F A 2 /F D 2 ratio will not be affected by the non-leading TMCs: Therefore, for our present study these leading-twist TMCs do not impact the results. However, if we were to examine the absolute structure functions (instead of ratios) or the absolute cross-sections, these TMCs must be incorporated in the high x, low Q region.

B. Higher Twist Corrections
At high x and low Q values, the sub-leading 1/Q 2 higher twist (HT) and residual power corrections can be important. To explore these effects, we will use a phenomenological x-dependent function taken from the CJ15 study [46] with the following form: where the higher-twist coefficient C HT depends on x and the nuclear A: with {h 0 , h 1 , h 2 } = {-3.3 GeV 2 , 1.9, -2.1}; h 0 carries units of GeV 2 , and the h 1,2 are dimensionless. Note that the A τ scaling was not part of the original CJ15 formulation and was included here to explore a possible nuclear dependence [47]. We have varied the τ exponent in our fits and the optimal values is indistinguishable from zero (τ ∼0 ± 0.005); hence, we find no significant A dependence for this correction. As we are primarily focused on lighter nuclei {He, Be, C} in this study, it may be useful to revisit the A τ dependence for the case of heavy nuclei.
Also note that since the h 0 term in Eq. (4) will generally be scaled by the target mass M , the separation between the higher-twist and the Target Mass Corrections (TMCs) discussed above is not unambiguous.
In Fig. 4 we display the higher twist correction as a function of x for a selection of Q values. The HT correction gives a slight reduction in the intermediate x region (x∼0.3), and then becomes large and positive at large x, depending on the Q value. For example, at x = 0.7, the correction factor is ∼1.5 for the Q = 1.3 GeV, whereas at Q = 2 GeV this is already reduced to ∼1.2. In Fig. 4 we also provide an inset plot to show the details in the intermediate x region. 3

C. Deuteron Nuclear Structure
A large fraction of the nuclear DIS data incorporated in global nPDF fits, including the present one, are expressed as a F A 2 /F D 2 ratio of structure functions on a heavy nuclear target and on a deuteron target, see Tables III-VI in the Appendix. A careful treatment of the denominator is therefore called for in order to analyze the effects of the nuclear dynamics on the nuclear quark and gluons, as most recently highlighted in Ref. [53].
The deuteron is the lightest and least bound of all compound nuclei. Therefore it is typically considered as composed of a free proton and a free neutron, and the deuteron structure function is computed as the sum 3 Note for the neutrino community.
Early theoretical considerations [48,49] of higher-twist (HT) effects indicated there could be real differences in these effects between weak and electromagnetic scattering. Indeed, early experimental investigations [50][51][52] of higher twist in neutrino-nucleon scattering support this conclusion and indicate that the HT contribution does not become positive at higher-x as in Fig. 4, but remains negative in the high x, lower Q 2 region. This neutrino HT behavior needs further exploration, however it suggests consideration of this difference in applying the results of this study to neutrino scattering in the relevant hi-x, low Q 2 region.
of the structure functions of its components. Extensive studies, however, indicate that the deuteron structure is far more complex than a simple combination of proton and neutron PDFs might suggest [33,46,[54][55][56][57][58]. Yet, the deuteron differs from the heavier targets considered in this paper because its binding energy is far smaller, making this a loosely bound system of two nucleons rather than a tightly bound and noticeably denser system of many nucleons interacting with each other [59]. Thus its nPDFs may not follow the A-scaling parametrization utilized in this work.
For these reasons -and to avoid double counting with the proton PDF fits that we use to anchor the nuclear PDFs parametrization used in this paper which already include deuteron target data -we choose not to fit the deuteron nPDF. Instead, we build upon the comprehensive studies available in the literature, and calculate the deuteron structure function, using the PDFs and nuclear dynamics simultaneously extracted in the recent CJ15 global QCD analysis [46]. At large x, this nuclear dynamics calculation is based on two main components: i) a baryonic smearing function evaluated on the basis of the AV18 nuclear wave function [60,61] that accounts for Fermi motion and nuclear binding effects, and ii) a fitted parametric function quantifying the effect of the offshellness of the bound nucleon. At smaller x, shadowing effects are calculated according to Ref. [62].
In Fig. 5 we display the ratio of F D 2 to the proton's F p 2 , which we use in the present nPDF analysis. For comparison, we also display the ratio to the isoscalar combination, F N 2 = F p 2 + F n 2 , which illustrates the size and x dependence of the nuclear effects in the deuteron system, after removing the underlying isospin symmetry effects. We can see that the deuteron-to-proton ratio, F D 2 /F p 2 , dips to ∼70% at x∼0.7 before dramatically increasing at larger x. Comparing this with the ratio to the isoscalar structure function, F D 2 /F N 2 , we see much of the previous effect was in fact due to the differing charge factors weighting the quark PDFs in the isoscalar combination, but the residual dynamical nuclear effects are not negligible, particularly at large-x values.
Note that the correct treatment of the deuteron structure function will impact, in principle, a large fraction of the DIS data, see Tables III-V. The nuclear dynamics in the deuteron may reasonably be neglected in typical nPDF fits because, after the usual kinematic cut on W 3.5 GeV, most of the DIS data lies in a small-to-medium x region where the corrections are small. However, if we want to extend our fit into the large x region, Fig. 5 demonstrates that it is essential to go beyond pure isospin symmetry considerations.
Operationally, we incorporate the deuteron's nuclear dynamics by rescaling the F A 2 /F D 2 ratio data by a F D 2 /F p 2 factor calculated using the CJ15 PDFs and deuteron a) We display the CJ15 calculation (Ref. [46]) of the structure function ratio of the deuteron F D 2 to the free proton F p 2 as a function of x for selected values of Q. Note that these curves are effectively fits to the experimental data. Additionally we display the ratio of F D 2 to the isoscalar combination F N 2 = F p 2 + F n 2 to demonstrate the nuclear dynamics in the deuteron beyond isospin symmetry effects. The determination of the neutron structure function in the denominator depends on the details of the nuclear correction model fitted in Ref. [46]. The colored bands are 90% confidence level uncertainty arising from the global CJ15 fit parameters.
correction model (Fig. 5a). Thus, the transformation is: and this effectively produces a nucleus-to-proton ratio data set that is then passed to the fit program. We provide the data and associated corrections in our supplementary material. Note that the F D 2 /F p 2 CJ rescaling factor is computed using the structure function results of the CJ15 global analysis [46,59]. CJ15 includes data from a variety of processes including DIS, Drell-Yan and jets from hadronhadron collisions at the Tevatron; this combination covers a broad kinematic range which contributes to improved precision. Additionally CJ15 includes data from the BONuS experiment [63,64] which measured the neutron to deuteron F 2 structure function ratio using an innovative spectator tagging technique. The BONuS data provides unique information on the neutron structure function, and this allows substantive improvement of the up and down quark flavor separation in the CJ15 analysis.
In principle we could also use CJ15 to correct to an isoscalar quantity (e.g., F D 2 /F N 2 ) as illustrated in Fig. 5-b). The difference between this approach and Eq. (5) depend on the details of the CJ15 PDFs, including the specific F 2 neutron structure function extracted. For this reason, we have chosen to implement Eq. (5); but, we have verified that this choice does not substantively impact our conclusions for the valence distributions, which is the focus of this study.
As in the DIS case, the p+D Drell-Yan cross section is modified by the nuclear dynamics in the deuteron target compared to calculations that treat the deuteron as composed of a free proton and a free neutron. The nuclear smearing model discussed above can also be applied to this case [65]. However, we have verified that the ensuing deuteron corrections amount to less than 1% in the kinematics covered by the available DY data, and can be neglected compared to the experimental uncertainties.

III. FITS TO JLAB DIS DATA SETS
Now that we have surveyed the various phenomena which enter the fit in the large x and low Q region, we will implement these effects into our fit so we can compare the relative effects.
In this section, we will first review the experimental data that enters the fit, and then provide an overview of the exploratory fits we will use to investigate the effects discussed in the previous section.

A. The Experimental Data
Historically many experiments have included isoscalar corrections (accounting a for different number of protons and neutrons in nuclei) in their data, e.g. structure function ratio data [28,66]. Since these corrections are model dependent, for the purpose of the current analysis, we have removed all isoscalar corrections from the used data sets following the specific prescriptions used by the collaborations. We provide supplemental materials with the tabulated data used in our analysis along with these corrections. This is available in a separate text file.

Data Used in nCTEQ15
For our fits, we will use the full set of DIS and DY data from the nCTEQ15 analysis [1] with isoscalar corrections removed. These data are listed in Tables IV-VI in  Appendix A. Additionally, the nCTEQ15 analysis includes RHIC inclusive pion data from PHENIX [67] and STAR [68] collaborations for nuclear modification R π DAu . The pion data helps to constrain the nuclear gluon PDF, especially in the lower x regions. As our focus will be on the high-x region, we have not included the pion data in the current fits which simplifies our analysis as we don't need to consider fragmentation functions. However, as a cross check, we separately compute the χ 2 values for the pion data to ensure our results are compatible with these data sets. These comparisons are reported in Sec. IV C 4 and show very good agreement with these data.

The JLab Data
In this analysis, we include new DIS data from the Jefferson Lab CLAS [38] and Hall C [69] experiments taken during the 6 GeV electron beam operation [70].
These data sets, in addition to spanning a wide range of A, provide high-precision constraints for nuclear PDFs at high-x and low-Q 2 . Nuclei included are 3 He, 4 He, 9 Be, 12 C, 27 Al, 56 Fe, 64 Cu, 197 Au, 208 Pb, spanning x from 0.2 to 0.95. Full details of these data sets are summarized in Appendix A as well as in Refs. [38,69].

B. The Parameters
As before, our nCTEQ PDFs are parameterized as [1]: where the nuclear A dependence is encoded in the c k coefficients as: 4 where k={1, ..., 5}. This parameterization has the advantage that in the limit A→1, c k (A) reduces to proton value p k .
The a k parameters control the magnitude of the modification to the p k proton "boundary condition," and the b k parameters control the power of the nuclear-A term. Specifically, within the nCTEQ15 framework we parameterize the combinations {u v , d v , (ū+d), (d/ū), s, g}, and we impose the boundary condition 5 of s=s=κ(ū+d)/2 with κ=0.5 at A=1.
We will then perform fits using a total of 19 free parameters. This includes 5  seven parameters for the gluon This set of 19 parameters includes the original 16 parameters used in the nCTEQ15 PDF fit, and adds the 3 parameters {a uv 3 , a dv 3 , a dv 4 } so that we have all the a k parameters open for the up and down valence to fully explore the parameter space. While our focus is on the high-x region where the up-and down-quark PDFs dominate, there is some interplay between the sea quarks and the gluon PDFs; hence, this will influence the g and d +ū parameters.
For our analysis we use our new nCTEQ++ framework which has a modular structure and links to a variety of external tools including a modified version of HOPPET [72] (extended to accommodate grids of multiple nuclei), APPLgrid [73], and MCFM [74]. Additionally, we have used FEWZ [75] for benchmarking our WZ calculations, and xFitter [76] for various cross checks.

C. The Fits
Having specified the required inputs, we now present a brief overview of the fits we will investigate. These are summarized in Table II and are outlined in the text below.
nCTEQ15*: This is the original set of nuclear PDFs as computed in Ref. [1]. The new data has not been fit; we have simply computed the χ 2 values.
BASE: This BASE fit will serve as our primary reference fit as it simply extends the nCTEQ15 fit by reducing the kinematic cuts to Q>1.3 GeV and W >1.7 GeV. In particular, this fit does not include any higher-twist corrections, or deuteron modifications.
HT: This fit extends the BASE fit by adding the higher twist (HT) power corrections as described in Eq. (3).
DEUT: This fit extends the BASE fit by including the deuteron modifications as described in Sec. II C.
nCTEQ15HIX: This fit extends the BASE fit by including both the deuteron correction (DEUT) and also the higher twist (HT) corrections.
Note, that all of the above fits include the TMCs originating from the scaling variable ξ and from the prefactor, as explained in the TMC discussion, Sec. II A. Table II also lists the {Q, W } cuts used in each fit. The original nCTEQ15 fit used cuts of Q>2.0 GeV, and W >3.5 GeV. For the new fits, we will relax these cuts to Q>1.3 GeV, and W >1.7 GeV to expand the kinematic range covered.
The W >1.7 GeV cut is chosen to avoid resonance contributions. However, there are indications that it may be possible to further reduce this cut to extend coverage of the Shallow Inelastic Scattering (SIS) region [77], and we examine this briefly in Sec. V. Now that we have outlined the fits, we present the results in the following section.

IV. RESULTS OF THE FITS
A. The χ 2 and χ 2 /N dof Values The χ 2 /N dof values provide a succinct measure of the quality of the fit, and these are summarized in Table II. Additionally, Fig. 6 displays the χ 2 /N dof for each individual experiment for the separate fits, and the experiment IDs are listed in Tables III, IV, V, and VI. nCTEQ15*: Starting with nCTEQ15, we first present the values for the original fit which used cuts of Q>2 GeV and W >3.5 GeV and obtained an excellent χ 2 /N dof of 0.81. However, if we take this nPDF and compute the χ 2 /N dof including the new data (and new less stringent cuts)-without fitting-we obtain a very large χ 2 /N dof of 1.70. Inspecting Fig. 6 we see that this large value reflects the fact that these new data, which are beyond the kinematic bounds of the original fit, are not well described. 6 These two entries in Table II provide a useful benchmark so we can see how our new fits improve in comparison.
HT: The HT fit applies the higher-twist correction of Sec. II B, and yields moderate improvement compared to the BASE fit (1489 vs. 1525), as well as associated improvements in the individual experiments of Fig. 6.
DEUT: The DEUT fit applies the deuteron correction of Sec. II C, and yields significant improvement compared to BASE (1331 vs. 1525). Additionally, in Fig. 6 we observe that the fit to the JLab data sets has improved and all the individual χ 2 /N dof values are now below 2.0, in contrast to both the BASE and HT fits.
nCTEQ15HIX: Finally, the nCTEQ15HIX fit applies both the deuteron correction and the higher twist power correction. This combination yields the best fit (1297 vs. 1525), and reduces the χ 2 /N dof to 0.84. Furthermore, examining the χ 2 /N dof values of Fig. 6 we see each experiment is well fit, with the exception of set 5108 which is a known outlier. 7 Comparing to the BASE fit, we observe that some of the original nCTEQ15 DIS data sets show a reduced χ 2 /N dof indicating the HT and DEUT corrections also improve the fit to these data.
Summary: In summary, we find that the HT modifications provide some improvements (∼3% of χ 2 /N dof ), the DEUT modifications provide larger improvements (∼10% of χ 2 /N dof ), and the combination of the HT and DEUT modifications yield the best fit with an improvement of ∼15% of χ 2 /N dof .

B. Data and Theory Comparison
Having evaluated the χ 2 /N dof for the separate experiments, we will examine some sample comparisons of our predictions with the data.
In Figure 7 we display the nucleus-to-deuteron F A 2 /F D 2 ratio for the selected data sets sorted by nuclei.   8. We display the F A 2 /F D 2 · F D 2 /F p 2 CJ ratio for selected data sets sorted by nuclei. We also overlay the theoretical prediction of nCTEQ15HIX in blue. The theory predictions have been calculated at averaged Q values where data sets overlap.
In Figure 8 we display the nucleus-to-proton (F A 2 /F D 2 )· (F D 2 /F p 2 ) CJ ratio, again sorted by nuclei. Here, we have multiplied by the ratio (F D 2 /F p 2 ) CJ taken from the CJ15 study [46], shown in Fig. 5, to approximately convert the results of the previous figure to ∼(F A 2 /F p 2 ). Note that the introduction of the x-dependent multiplicative (F D 2 /F p 2 ) CJ factor visually suppresses the A-dependent change in slope seen in Figure 7. However, a check of the values of (F A 2 /F p 2 ) at x∼0.3 and x∼0.7 for 4 He and 197 Au confirms that the A-dependent change in slope has been maintained.
In Fig. 8 we also display the corresponding theoretical calculations (blue line) obtained with the nCTEQ15HIX PDFs. We can see that they provide a very good description of the fitted data. Iron PDF Ratios to nCTEQ15 (Q = 2 GeV) FIG. 11. Iron ( 56 Fe) PDFs ratio compared to nCTEQ15 at Q = 2 GeV. We show the uncertainty bands for nCTEQ15 (blue) and nCTEQ15HIX (yellow) computed with the Hessian method. Lead PDF Ratios to nCTEQ15 (Q = 2 GeV) FIG. 12. Lead ( 208 Pb) PDFs ratios compared to nCTEQ15 at Q = 2 GeV. We show the uncertainty bands for nCTEQ15 (blue) and nCTEQ15HIX (yellow) computed with the Hessian method.

C. The PDFs
Finally, we examine the impact of the various corrections on the underlying PDFs. Although the up and down flavors dominate in the high-x region, these PDFs will also influence the gluon and sea-quark distributions indirectly. Thus, we must examine all the flavors to obtain a complete picture of the impact of the JLab data.
We will present a sample of three representative nuclei: 56 Fe, 208 Pb}. As the bulk of the new data is at lower A values we will focus our attention on the 12 C PDFs, but we will want to examine the heavier nuclei to infer the various trends in A as well as to verify that our extrapolations to Pb are reasonable. For these plots, we have also displayed the uncertainty bands computed with the Hessian method for both nCTEQ15HIX and the original nCTEQ15 PDF sets. Note, that when plotting nuclear PDFs, we always display the full nuclear distribution, f A , constructed out of the bound proton (neutron) distributions, f p(n)/A , as: with f n/A constructed using the isospin symmetry.
In Figs. 9 and 10 we present the results for 12 C. Figure 9 shows the magnitude of the PDF x f (x, Q) at Q = 2 GeV on a log-scale. We observe that the separate fits are generally similar across the full x kinematic range and that there is no abnormal behavior in either the low or high x region. Additionally, this figure underscores the extent to which the PDFs are decreasing as we move to larger x, especially for the gluon and sea-quarks {g,ū,d, s=s}; obtaining large event statistics in the highx region is challenging.

The Small x Region
In Fig. 10 we display the ratio of the fits compared to the nCTEQ15 nPDFs, we use here a log-linear scale to better illustrate the full x range. Our first observation is that all of our new fits to the JLab data yield very comparable results in the small x region. This is a reasonable result as the different corrections we are applying in the fits (e.g., HT, DEUT) all primarily impact the large x region. We also observe that all these fits are uniformly below the nCTEQ15 PDFs by ∼6%.
Thus, it appears the fits are effectively hardening the high-x PDFs of select parton flavors (especially those of the u-and d-quark), while somewhat suppressing PDFs in the small-x region. This behavior suggests increased shadowing and merits further investigation.

The Large x Region
From Fig. 10 it is clear that the fits are hardening the nuclear PDFs at larger x to better describe the new data sets. In particular, we see the up, down and gluon PDFs all are larger than nCTEQ15 while the sea quarks {ū,d, s=s} are smaller. For example, if we look in the region of x∼0.7 we see theū andd PDFs are reduced by as much as ∼10% while the gluon is increased by ∼5% and the up and down PDFs are increased by ∼10%.
Comparing the individual fits, we see a very clear pattern in theū andd PDFs. As we move from the nCTEQ15 fit, to the BASE, HT, DEUT, and finally the nCTEQHIX fit we see theū andd PDFs are uniformly decreasing in the large x region. Recall that this ordering of the fits also shows a monotonically decreasing value of χ 2 /N dof , c.f., Table II.
We also observe that the BASE and the HT fits are quite similar, as the addition of the HT corrections, that are substantially smaller for structure function ratios than for either the numerator or the denominator, see Similarly, the DEUT and the nCTEQ15HIX are also quite similar as these fits only differ by the inclusion of the HT corrections (for the nCTEQ15HIX fit). This difference is barely discernible in the figures as the curves are generally within a percent, but nCTEQ15HIX generates a 3% improvement in χ 2 /N dof compared to DEUT. These two fits also produce the largest downward shift of theū,d and strange PDFs while yielding a smaller increase of the up and down PDFs. Recall that DEUT and nCTEQ15HIX have the lowest χ 2 /N dof values.

Large A Results
Having investigated the carbon results in detail, we now want to explore the larger A values to see how the above effects scale to heavier nuclei. In Figures 11 and 12 we display the PDF ratios (compared to nCTEQ15) for iron and lead PDFs on a log-linear scale.
The small x behavior is quite similar to the carbon PDF result. At x=10 −3 all the PDFs are shifted downward by ∼6% to ∼8%, and the shift of the individual fits is quite uniform.
In the large x region, the behavior is also similar to the carbon PDFs, but the size of the shifts is slightly larger. For example, at x∼0.7 theū andd PDFs for carbon were shifted downward by ∼10% compared to nCTEQ15, while the iron and lead are shifted by up to ∼30%. We again see that the gluon, up and down PDFs are increased at large x. Again comparing at x∼0.7, while the up and down PDFs of carbon are increased by ∼10%, the iron and lead are shifted by up to ∼20%.
Clearly the new data has significant impact on the PDFs in this kinematic region. We also observe that the impact of the individual corrections {HT, DEUT} are rather uniform, and we obtain the best fit using both the HT and DEUT corrections which give us the nCTEQ15HIX fit. Therefore in the following, we will now focus on this fit.

PDF Uncertainties
Finally, we now compare the uncertainties of the nCTEQ15 and nCTEQ15HIX PDF sets as displayed in Figs. 9-12. These results should be interpreted carefully as we must be attentive to parameterization bias and other imposed constraints. For example, our parameterization ties the strange PDF to the average (ū +d) distribution, so this is not a fully independent distribution. 8 Examining Fig. 9 showing the 12 C PDFs, both nCTEQ15 and nCTEQ15HIX have comparably sized uncertainty bands, with the exception of the gluon. The larger uncertainty band for the nCTEQ15HIX gluon may be more representative of the true uncertainties as it can be difficult to fully explore the loosely constrained gluon in the small x region where minimization methods often encounter troublesome flat parameter directions. In this sense, the larger gluon uncertainty band of, for example, the EPPS16 nPDF set [4] may be more characteristic (c.f., Fig. 18).
Note that the nCTEQ15 fit includes data on single inclusive pion production from the STAR and PHENIX experiments, which we do not include in the nCTEQ15HIX fit. However, we have verified these experiments are compatible with our nCTEQ15HIX fit, and they yield a χ 2 /N dof ∼ 0.5. The pion data can help constrain the nuclear gluon PDF, especially in the midto lower-x region, and to some extent, the slight increase of the gluon uncertainty of nCTEQ15HIX can thus be attributed to this.
Turning to the PDF ratio plots of Figs. 10-12, again these should be interpreted carefully, especially in the large x region where we have a ratio of small numbers. We observe that both the up and down PDFs appear to have reduced uncertainty bands in the larger x region. This makes sense as the bulk of our new data is constraining these flavors in this kinematic region. The gluon uncertainty is increased compared to nCTEQ15, presumably for the same reasons discussed above. Theū andd uncertainties are roughly comparable between the nCTEQ15 and nCTEQ15HIX, as we expect only indirect constraints on these flavors.
While the uncertainties of the deuteron structure functions displayed in Fig. 5 are not used in the computation of the error bands, their impact can be assessed by comparing the fits without (BASE, HT) and with (DEUT, nCTEQ15HIX) the deuteron corrections. 8 Note, our recently released nCTEQ15WZ nPDF analysis included LHC W/Z production data which allowed to provide additional freedom (and constraints) for the strange PDF [71].
These differences lie within the uncertainty bands in all cases, except for the up and down nPDFs at intermediate x values (x∼0.5); hence, our uncertainties in this region may be slightly underestimated. Finally, the strange PDF also shows a reduced uncertainty band. As mentioned in Sec. III B, this distribution is tied to the (ū +d) sea quarks and we have not introduced any free strange parameters. For a more representative uncertainty estimate of the strange PDF, compare with the error bands of nCTEQ15WZ (c.f., Ref. [71], Fig. 12).

D. Parameter Scans
We have performed parameter scans for the 19 parameters of the nCTEQ15HIX fit, and selected results are displayed in Figs. 13 and 14. As expected, the JLab data exerts strong constraints on the parameters of the up and down quarks, while the impact on the sea quarks and gluon is reduced. We also note that for the sea quarks (not shown) and gluon, the value of the parameters at the minimum χ 2 is quite similar for both the nCTEQ15 set and the new JLab set. However, for the up and down parameters, where the JLab data has more constraining power, we find in some cases (e.g., a uv 2 ) the fit minimum is a compromise between the two data sets. This is a common feature in global PDF fits with large data sets. To gain additional insight on which data sets are driving this result we now examine Figs. 13 and 14 in detail.
In Fig. 13 we display the profile for the a uv 2 parameter. The first plot, Fig.13-a), shows the total χ 2 shift, as well as the individual shifts due to the original nCTEQ15 data set and the new JLab data set. Clearly the new JLab data help constrain this parameter compared to the nCTEQ15 data alone. We also observe the original nCTEQ15 data prefers slightly larger values for the a uv 2 parameter as compared to the JLab data; this type of behavior is common for global PDF fits. In Fig. 13-b) we display the χ 2 shift for the most significant data sets constraining the a uv 2 parameter. We note the data sets with the largest impact are the 56 Fe/D data which is pulling a uv 2 to higher values, and the 12 C/D data which is pulling a uv 2 to lower values. In Fig. 14 we display the profile for the a g 1 parameter. Fig.14-a) again shows the total χ 2 shift, as well as the individual shifts due to the original nCTEQ15 data set and the new JLab data set. We observe the new JLab data provides a weaker constraint on the gluon as compared to the valence quarks. We also note the preferred minima of both the original nCTEQ15 data and the JLab data coincide for the gluon parameter. Turning to Fig. 14-b) we see the χ 2 shift for the most significant data sets constraining the a g 1 parameter. In contrast to the above case of a uv 2 , we note the data sets with the largest impact are 119 Sn/ 12 C and 208 Pb/D; in this case, the 56 Fe/D data has less impact on this parameter. 14. PDF parameter profiles of the the gluon parameter a g 1 . Fig. a) We display the χ 2 shift due to the original nCTEQ15 data set (blue dots), the new JLab data set (red dashes), and the total (black solid). Fig. b) We display the χ 2 shift for the most significant projectile-target combinations constraining the a g 1 parameter.

V. THE LOW-W KINEMATIC REGION
The W >1.7 GeV cut on the data fitted in the nCTEQ15HIX analysis discussed above was imposed to restrict the study to the DIS regime, where we can reasonably assume QCD factorization to work, and avoid resonance excitations in electron-scattering events.
However, it has been shown that the characteristic x dependence of nuclear structure-function ratios displayed, for example, in Fig. 7 persists into the resonance region at low W values. Specifically, Fig. 3 of Ref. [41] overlays data from the resonance region (1.1 < W < 1.7 GeV, Q∼2 GeV) with W > 1.7 GeV DIS data from the SLAC E139 [28], SLAC E87 [30], and BCDMS [32] experiments, and finds that "the agreement of the resonance region data with the DIS measurement [...] is quite striking." This is likely a manifestation of the quark-hadron duality phenomenon, abetted in the case of nuclear scattering by Fermi motion of the bound nucleons, which results in a further, effective averaging of the nucleon structure functions over multiple resonances; together, these dynamics may permit a description of nuclear structure functions in terms of partonic degrees-of-freedom, even in kinematic regions where resonance excitation is the dominant effect. It is therefore interesting to remove the W cut, and explore to what extent the resonance region data can be described in terms of nPDFs. A. nPDFs at large x An important difference between nuclear PDFs and proton PDFs is that, in principle, the momentum fraction x of a parton with respect to the average nucleon momentum in a nucleus can range from [0, A] whereas that of a proton is confined to lie in the range [0, 1]. Phenomenologically, at the photon energies under consideration, nuclear partons belong to bound nucleons and share a relative fraction x N ∈ [0, 1] of their momentum. As the nucleons interact among themselves, they exchange momentum; thus, their partons can receive a boost in the photon's direction, causing their momentum fraction x to shift and exceed 1. Since the nucleon-nucleon interaction is on average soft, the shifts are moderate and the nuclear parton distributions are still predominantly within x ∼ [0, 1] with only a small tail exceeding 1.
This effect is currently not captured in the nPDF parametrization adopted in this paper, see Sect. III B, which is confined to x ∈ [0, 1] and needs to be suitably extended. In fact, a full description of nPDFs in the x ∈ [0, A] range would also entail a generalization of DGLAP evolution equations and careful consideration of the sum rules [78]. 9 In this paper we limit ourselves to an exploratory study and demonstrate, as a proof of principle, that a simple generalization of the adopted parametrization can lead to a reasonable description of electron-nucleon collisions in terms of nuclear PDFs even in the resonance region. We will leave more detailed studies for future work.

Convolution kernel
The parton momentum shift caused by nucleonnucleon interactions can be effectively described by a convolution of the nucleon-level PDF f with a smearing function S A : where S A represents the probability that the parton f belongs to a nucleon of momentum fraction y compared to the average nucleon momentum. Such representation is, in fact, grounded on the so-called Weak Binding Approximation to the calculation of nuclear structure functions, where S A can be obtained starting from the nucleon wave function in the corresponding nucleus [79][80][81][82][83][84]. In the Weak Binding Approximation, the smearing function is in general scale-dependent at finite Q values. Therefore f A does not satisfy the DGLAP evolution equations except at Q → ∞ where S A becomes scale independent, or at small x where the smearing effect is marginal. In the context of nPDF fits, however, we rather regard the convolution of Eq. (9) as a generalization of the nPDF parametrization discussed in Sect. III B, and would only apply it at the initial scale Q 0 .
At x values not exceeding 1, the nucleon-nucleon interaction is soft and the smearing function can be approximated by a normalized Gaussian, where the effect of nucleon-nucleon interactions is described by the Gaussian width ∆ A . Furthermore, the nucleon binding causes a (small) shift of the Gaussian peak towards smaller x values, and this is controlled by the δ A term. Note that the larger the nucleus, the more the nucleons that interact with each other (hence a larger ∆ A ) and the larger the binding energy (hence a larger δ A ).
The net result of this convolution is to deplete the partons in the intermediate x region, and enriching the large x 1 region. Hence the nuclear ratio F A 2 /F B 2 with A > B will also display a dip below unity at intermediate x values (the EMC region) and a relatively steep rise above unity in the x → 1 Fermi motion region, as seen experimentally. With a suitable choice or fit of the ∆ A and δ A parameters one can therefore expect to quantitatively describe the observed nuclear effects in the Fermi motion region 10 .
As an example, in Fig. 15 we show the ratio of convoluted to unmodified PDF for a toy model f (x) = x −1.5 (1 − x) 3 , which reproduces the main features of the up and down quark PDFs. The parameters have been chosen to approximate the behavior of the Carbon structure function data.

x-rescaling
While the smearing approach has a close connection to the underlying nuclear dynamics and guarantees that the PDF sum rules are satisfied, the convolution is computationally expensive. Hence, in this exploratory study we will try a simple remapping of the parametrization considered so far, by an A-dependent shift of the x variable.
We want this remapping to approximate the rise of F A 2 /F p 2 in the large x region without distorting the phenomenology in the intermediate to small x region. To ensure these conditions are satisfied, the remapping should (i) only impact the very large x region (x∼1), and (ii) ensure the various momentum sum rules are satisfied within uncertainties.
The specific remapping we shall use is 10 In principle, by generalizing the Gaussian ansatz in Eq. (10) for the smearing function and including, for example, power law tails [85], the convolution formalism of Eq. (9) can also be used to extend the analysis into the x>1 region. We leave this for future work.
where f A are the PDFs determined in the nCTEQ15HIX fit, and As the PDFs are typically decreasing functions of the parton's fractional momentum, the negative shift ensures that the transformed function is larger than the unmodified one, and non-vanishing as x → 1. The ε parameter controls the overall size of the rescaling effect, and the x κ factor ensures that only the large-x region is modified. The log 10 A term ensures an increasing modification across the full range of nuclear A values from proton (A = 1) to lead (A = 208).
In Fig. 16, we show the impact of the rescaling on the ratio f A /f p for A = 208 and a selection of ε and κ values. For simplicity, we have used here the same f (x) = (1 − x) 3 /x 1.5 toy PDF considered in the nuclear smearing discussion. Clearly, the remapping does not preserve the PDF normalization within the x ∈ [0, 1] range, and therefore breaks the PDF sum rules. However, in all cases the breaking amounts to a small increase of the momentum sum by less than 1%, which is well within the uncertainty of the nPDFs.

B. nPDFs in the x→1 Limit
Overall, the remapping accomplishes our stated goals, and can be used to investigate how well it can describe the data at the highest x, that was excluded from our nCTEQ15HIX fit due the W > 1.7 GeV cut.
To simplify the discussion, we fix the small-x suppression parameter to k = 10 and will study the effect of varying the parameter. One can appreciate that this choice is reasonable by looking at Fig. 15, where the A = 12 remapped f A /f p , plotted in red, is compared with the blue convolution model discussed earlier. By design the remapping only deforms the PDF at large x, and does not reproduce by itself the intermediate-x dip displayed by the convolution model. This however is small and can be easily taken care of in a fit by the standard nPDF parameters discussed in Sect. III B. On the contrary, the large-x steep rise is pretty well approximated by the remapped PDF. We have furthermore verified that k = 10 is also a good choice for other nuclear targets, given the assumed logarithmic A dependence of the rescaled x A momentum fraction.
We can now investigate the impact of rescaling the nCTEQ15HIX nPDFs and compare these to the highx data. In Fig. 17 we display a selection of data sets that have data points in the x → 1 region, and compare them to calculations performed with shifted nPDFs for fixed κ = 10 and a representative choice of parameters in the [0, 0.1] range. Note that in the 56 Fe and 64 Cu panels, the empty symbols represent data points at W < 1.7 GeV, that were excluded from the nPDF fit. We also display the computed χ 2 values in the plot legend. These are computed both for the points that satisfy the W > 1.7 GeV cut (labeled χ 2 cut ), and for all the points (labeled χ 2 nocut ). Note that for these χ 2 calculations we have simply applied the x rescaling when computing F 2 (x, Q) using the nCTEQ15HIX nPDFs; this is not a separate fit.
If we focus on the top two panels of Fig. 17, 12 C and 27 Al, we see all these data points satisfy the W > 1.7 GeV cut. Turning our attention to the bottom two panels of Fig. 17, 56 Fe and 64 Cu, some data points do not satisfy the W > 1.7 GeV cut, and these points give rise to the differing χ 2 cut and χ 2 nocut values. Overall, we observe there is a small but significant variation in χ 2 cut , and a dramatic variation in χ 2 nocut . In both cases, the fit shows preference for ≈ 0.03. With this parameter, we can describe very well the turnover of the structure function ratio at x ∼ 0.7, and the subsequent large-x rise over the whole A range. Overall, we obtain a very reasonable description of the large-x data in the resonance region.

C. Extending nPDF fits to low W
We have found that the x rescaling provides a simple means of mimicking the behavior of the nuclear PDF in the x → 1 region. While we have used it here postfacto, it can also be effectively utilized to generalize the nPDF parametrization discussed in Sect III B with only 2 additional free parameters. As the x rescaling already increases the F A 2 /F D 2 ratio in the Fermi motion region by about the right amount compared to the baseline nCTEQ15HIX nPDFs, the remaining parameters will more easily adjust to the data, and we can expect an improvement in the fit quality in the high x, low W region.
Furthermore, this analysis suggests that it may be possible to expand the kinematic reach of the fit in the {x, Q}-space to include data in the low W < 1.7 GeV resonance region. This is an interesting possibility, and we reserve the details for a future study.

VI. COMPARISONS WITH OTHER nPDFS
In this section we compare our new nCTEQ15HIX nPDFs with other results from the literature.
Examining the ratios of Fig. 19 we note some features similar to the comparison between nCTEQ15 and nCTEQ15HIX of Fig. 10. Specifically, for the up, down, and gluon PDFs, nCTEQ15HIX is generally above the other PDFs in the high x region. Correspondingly, for the anti-up and anti-down PDFs, nCTEQ15HIX is Carbon PDF Ratios to nCTEQ15HIX (Q = 2 GeV) FIG. 19. We compare our nCTEQ15HIX results to other nPDF sets from the literature including EPPS16 [4], nNNPDF2.0 [6], and TUJU19 [5]. We plot the nPDF ratio for carbon 12 C at Q = 2 GeV compared to nCTEQ15HIX on a log-linear scale.
generally below the other PDFs in the high x region. This is similar to the pattern observed in Fig. 10, although these difference are well within the uncertainty bands. It will be of interest to follow these differences as the uncertainty bands are reduced in future analyses. We do not display a comparison with our recent nPDFs nCTEQ15WZ [71] where the nCTEQ15 analysis was supplemented by the available W/Z boson production data from the LHC. This is because here we concentrate on the high-x region of quark distributions, and the W/Z data were mostly relevant for the strange and gluon distributions at lower x values. Hence, they had minimal impact on the high-x behavior of the nPDFs. Therefore, the high x behavior of the nCTEQ15WZ PDFs was quite similar to the original nCTEQ15 nPDFs, and this was already presented in Figs. 9-12.

VII. CONCLUSION
PDFs encode the dynamics of the strong interaction and provide a crucial link between experimental measurements and theoretical models.
Thus, the updated nCTEQ15HIX nPDFs presented here will contribute to increased precision of experimental analyses, which can thereby yield further insights into the QCD. The novel aspects of this study consisted of several elements.
The newly-included JLab data provide important additional constraints on the nuclear PDFs in the high-x and low-Q regime. We have explored this region with the nCTEQ15HIX set by relaxing our kinematic cuts to Q>1.3 GeV and W >1.7 GeV.
In the high-x region, there are various theoretical corrections which must be taken into account, especially for the Q ∼ few GeV region now included with our relaxed kinematic cuts. Among these, the Higher-Twist (HT) modifications discussed in Sec. II B reduce χ 2 /N dof by ∼3%, while the Deuteron (DEUT) modifications of Sec. II C reduce χ 2 /N dof by ∼10%; ultimately, the combination of both (which was used in nCTEQ15HIX) reduces χ 2 /N dof by ∼15%.
We also investigated the low-W kinematic region including the impact of extending the nuclear PDFs into the region of very large x.
We used a xrescaling technique that mimics the parton momentum shift to higher x induced by the nucleon's Fermi motion.
We find that the x-rescaling provides a simple (albeit approximate) means to extend the nPDF parametrization to larger x values than currently included in our main fit, and to address data falling within the resonance region. Calculations based on a simple two-parameter x rescaling were indeed able to improve the χ 2 /N dof for data at large-x not only inside our current cuts, but -crucially and more markedlyalso for those data that are currently outside and do not contribute to the nPDF fits. This preliminary analysis suggests it may be possible to further reduce our W cut, extend the fits into the resonance region, and further improve the precision of the nPDFs at large momentum fraction x.
The structure function ratio (F A 2 /F D 2 )·(F D 2 /F p 2 ) CJ (Fig. 8) displays the characteristic EMC shape, and the qualitative x dependence of the nuclear ratio in the EMC Region was found to be consistent for the range of A considered in this analysis.
Our parameter scans (Figs. [13][14] show that while the new data provide substantive constraints on the up-valence and down-valence distributions, there are nonetheless some tensions between individual data sets. This behavior is common for global PDF fits, but some of the details, such as the tension between the Fe and C data sets (Fig. 13) may warrant additional investigation.
A general feature of the new nCTEQ15HIX fit is a relative hardening of the nuclear PDFs involving shifts of select distributions toward higher x, an effect realized in the enhancement of the up-valence, down-valence, and gluon distributions, and corresponding depletion of the sea-quark PDFs, {ū,d,s}.
This pattern is apparent when comparing nCTEQ15HIX to nCTEQ15 (Figs. 9-12). Additionally, we find similar behavior when comparing nCTEQ15HIX to other nPDFs from the literature including EPPS16, nNNPDF2.0, and TUJU19 (Figs. 19). While these variations are within the nPDF uncertainty bands, it will be interesting to see to what extent this behavior persists with improved precision.
As PDF global analyses represent a primary computational tool available to describe hadronic interactions and structure in the context of QCD, the nCTEQ15HIX nPDFs will enable a new level of precision in the exploration of nuclear dynamics and particle phenomenology. In a broader sense, these tools can validate features of the standard model to nextgeneration precision, aid in the search for discrepancies which may signal undiscovered phenomena, and thereby yield deeper insights into the QCD theory and the structure of hadronic matter.   Table III we indicate the total number of data in the sets, and the remaining number after the nCTEQ15HIX cuts. In Tables IV, V and VI we indicate the total number of data in the sets, and the remaining number after the nCTEQ15 kinematic cuts (first) as well as the nCTEQ15HIX cuts (second).
NMC F D 2 : The NMC F D 2 data set ID=5160 of Table IV is the one case which is not a ratio of structure functions. It serves as a cross-check of how well we describe absolute structure functions in the deuteron for which we apply a dedicated correction, instead of fitting as is done for heavier nuclei. Recall from Sec. III B that our proton parameters p k are fixed, and for the special case of the deuteron (A = 2) we construct this using isospin symmetry. We then apply the deuteron structure function modifications described in Sec. II C for the DEUT and nCTEQ15HIX fits, and find this improves χ 2 /N dof for F D 2 by ∼1%.
Data Files: We have provided a full list of the data used in this fit in the text file "fitted_data.txt" which is included in the arXiv version of this document. 11 11   IV. The DIS F A 2 /F D 2 data sets used in this fit. The table details the specific nuclear targets, the references, and the total number of data points in the set. The last column shows first the number of data points remaining with the nCTEQ15 cuts of Q>2 GeV and W >3.5 GeV, and second the number of data points remaining with the nCTEQ15HIX cuts of Q>1.3 GeV and W >1.7 GeV. The NMC data set ID=5160 is F D 2 ; this is the one case which is not a ratio.