Nuclear Parton Distribution Functions

• Nuclear parton distribution functions are QCD distributions that describe the momentum distribution of quarks, antiquarks, and gluons in bound nucleons and nuclei.
• They are determined through global fits using data from deep-inelastic scattering, Drell–Yan processes, and hadronic collisions, incorporating constraints like sum rules and collinear factorization.
• They are essential for precision collider predictions, helping to disentangle initial-state nuclear effects from final-state dynamics in experiments at RHIC and the LHC.

Searching arXiv for the cited nPDF literature to ground the article in recent and historical papers. Nuclear parton distribution functions (nPDFs) are the QCD distributions that describe how quarks, antiquarks, and gluons are arranged inside bound nucleons and nuclei as functions of Bjorken $x$, the hard scale $Q^2$, and the nuclear species $(A,Z)$. In the collinear-factorization framework, they are the nuclear analogues of ordinary proton PDFs, but with nuclear modification factors that encode shadowing, antishadowing, the EMC effect, and Fermi motion. They are indispensable both for describing hard processes in nuclei and for controlling nuclear corrections in free-proton PDF extractions, especially in flavor separation and in the strange sector relevant for precision $W$ and $Z$ phenomenology at the LHC (Ethier et al., 2020, Kovarik, 2010, Klasen, 2024).

1. Formal definition and factorized description

A standard definition writes the nuclear PDF for a nucleus $(A,Z)$ as the proton–neutron average

$$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$

with bound-proton and bound-neutron PDFs related to free-nucleon PDFs through a nuclear modification factor, for example

$f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$

or, in another common convention,

$$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$

These relations separate isospin effects from genuine nuclear modifications and define the basic objects used in phenomenology (Ethier et al., 2020, Zurita, 2018, Klasen, 2024).

The perturbative framework is the same as for free hadrons. For deep-inelastic scattering (DIS), the nuclear structure function factorizes as

$$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$

with Wilson coefficients $Q^2$0 that are the same short-distance functions used in proton scattering. The scale dependence is governed by the DGLAP equations, so the formal distinction between proton PDFs and nPDFs lies in the input conditions rather than in a different evolution law (Klasen, 2024, Ethier et al., 2020).

Sum rules remain essential constraints. Global fits impose valence-number and momentum conservation, and many parameterizations are constructed so that the $Q^2$1 limit reproduces the free-proton PDF baseline. This $Q^2$2 boundary condition is explicit in several modern frameworks and is especially important for controlling light-nucleus behavior (Kovarik, 2010, Khalek et al., 2019, Helenius et al., 2021).

2. Characteristic nuclear modifications and their interpretation

The canonical $Q^2$3-space pattern consists of four regions. At small $Q^2$4, roughly $Q^2$5, nPDFs are typically suppressed relative to free nucleons; this is shadowing. In the region $Q^2$6, many fits display antishadowing. At $Q^2$7, the EMC effect gives another suppression, while at $Q^2$8 Fermi motion produces a rise. This qualitative structure is well established in fixed-target DIS ratios and is a standard organizing principle of nPDF phenomenology (Klasen, 2024, Paakkinen et al., 7 Jun 2026, Paukkunen et al., 2010).

The physical origin of these regions is not uniform. Fermi motion is comparatively well understood in terms of nucleon momentum distributions. The EMC effect is experimentally robust but, as emphasized in recent reviews, its microscopic origin is still not fully understood and may involve both partonic and hadronic mechanisms. Shadowing at small $Q^2$9 is the most theoretically differentiated regime: leading-twist Gribov–Glauber approaches relate it to diffraction and coherent multiple scattering, dipole models represent it through eikonal absorption of $(A,Z)$0 dipoles, CGC descriptions emphasize nonlinear small-$(A,Z)$1 dynamics and the saturation scale, and higher-twist recombination approaches treat it as a nucleus-enhanced power correction (Paakkinen et al., 7 Jun 2026, Klasen, 2024).

A recent light-cone model gives a different microscopic decomposition, writing the nuclear PDF as a sum of uncorrelated nucleon and correlated nucleon-pair terms,

$(A,Z)$2

with $(A,Z)$3 a nucleon-pair PDF and with the correlated contribution used to account for the EMC effect and for nuclear corrections to the Paschos–Wolfenstein ratio. This is a model construction rather than a global-fit standard, but it illustrates a current line of work that connects nuclear modifications to nucleon–nucleon correlations (Yang et al., 2024).

The observed modifications are flavor dependent and evolve with $(A,Z)$4. Valence distributions are generally the best constrained; sea quarks are less precise; gluons remain the least constrained component over much of $(A,Z)$5, especially at small $(A,Z)$6. DGLAP evolution tends to reduce the magnitude of the nuclear effects at high $(A,Z)$7, pushing ratios closer to unity as the scale increases (0709.3038, Kovarik et al., 2013, Paukkunen et al., 2010).

3. Parameterizations, global-fit strategies, and uncertainty estimation

Global nPDF analyses follow the same general logic as proton PDF fits: choose a parameterized input distribution at some $(A,Z)$8, impose sum rules and isospin constraints, evolve in $(A,Z)$9, compute observables, and minimize a global $W$0. The main distinction is the explicit dependence on $W$1, which can be introduced either through multiplicative nuclear modification factors or through $W$2-dependent PDF parameters (Kovarik, 2010, Ethier et al., 2020, Klasen, 2024).

One common class of fits writes

$W$3

with a weight function such as

$W$4

possibly with further $W$5-dependence assigned to the coefficients. This strategy underlies HKN-type and KA-type analyses and provides a direct flavor-by-flavor modification relative to a free-nucleon baseline (0709.3038, Khanpour et al., 2016, Tehrani, 2017).

A second class parameterizes the bound-proton PDFs directly, with coefficients promoted to $W$6-dependent functions such as

$W$7

This approach is characteristic of CTEQ-derived and TUJU-style analyses and has the practical advantage that the proton limit is built into the same functional family used for nuclei (Kovarik, 2010, Kovarik et al., 2013, Walt et al., 2019, Helenius et al., 2021).

A third strategy reduces parametrization bias by using neural networks and Monte Carlo replicas. The first NNPDF-style nPDF fit, nNNPDF1.0, parameterized the singlet, octet, and gluon combinations with a single $W$8 network in $W$9, $Z$0, and $Z$1, imposed the $Z$2 boundary through a penalty term tied to NNPDF3.1, and validated the procedure with closure tests (Khalek et al., 2019).

Perturbative accuracy ranges from LO to NNLO. Early global fits were frequently LO or NLO; later analyses moved to NNLO for DIS- and DY-dominated datasets; and contemporary fits increasingly combine fixed-target and collider observables with general-mass variable-flavour-number schemes such as FONLL, while some older NNLO fits used ZM-VFNS. Open-source frameworks including xFitter and APFEL have been extended to nuclear applications, and fast-grid methods are now used for NNLO electroweak-boson observables (Tehrani, 2017, Walt et al., 2019, Khanpour et al., 2020, Helenius et al., 2021).

Uncertainty propagation is dominated by the Hessian method, although Monte Carlo is also used. The adopted tolerances differ substantially across groups. EPS09 defined error sets by $Z$3 at 90% confidence level, KA15 used $Z$4, KSASG20 used $Z$5 for 68% confidence bands, and TUJU19 used $Z$6 for the nuclear fit. These choices, together with the flexibility of the input parameterization and the dataset composition, are a major source of the spread among published uncertainty bands (Paukkunen et al., 2010, Tehrani, 2017, Khanpour et al., 2020, Walt et al., 2019).

4. Experimental constraints and representative global analyses

Historically, the backbone of nPDF extraction was fixed-target charged-lepton DIS and Drell–Yan. HKN07 performed LO and NLO global fits to $Z$7 and $Z$8 ratios, finding that valence quarks are well determined, antiquarks are determined mainly at $Z$9, and gluon modifications cannot be fixed with the then-current scaling-violation data (0709.3038). KA15 later extended this fixed-target program to NNLO with $(A,Z)$0 points, explicit $(A,Z)$1-dependent coefficients, and Hessian uncertainties, obtaining $(A,Z)$2 and providing a complete $(A,Z)$3 set (Khanpour et al., 2016).

The inclusion of gluon-sensitive hadronic data marked a second stage. EPS08 supplemented DIS and DY with inclusive high-$(A,Z)$4 hadron production in $(A,Z)$5Au collisions at RHIC, especially forward BRAHMS data, and found markedly stronger gluon shadowing than earlier global analyses while maintaining a very good simultaneous description of the data (0802.0139). nCTEQ15 incorporated 616 charged-lepton DIS points, 92 DY points, and 32 RHIC pion-production points, obtained $(A,Z)$6, and emphasized that allowing independent $(A,Z)$7 and $(A,Z)$8 nuclear corrections affects bound-proton PDFs more strongly than the full nuclear combinations (Kusina, 2016).

A broader universality-oriented NLO fit was presented in 2012, combining charged-lepton DIS, DY, neutrino DIS, and inclusive pion production in $(A,Z)$9Au. With $$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$0 data points and $$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$1, it found that one set of universal nuclear modification factors reproduces the main features of all included datasets without significant global tension (Florian et al., 2012).

Modern analyses increasingly incorporate collider measurements. EPPS16 enlarged the parameter freedom from 15 to 20 parameters and, for the first time in that series, included LHC $$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$2 data in addition to NC DIS, CC DIS, DY, and RHIC pion production. A central conclusion of that work was that including charged-current data did not introduce significant tension, while the more flexible parameterization broadened the uncertainty bands (Zurita, 2018). TUJU21 added $$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$3 $$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$4Pb $$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$5 points on top of DIS data and found that the total $$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$6 improves from $$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$7 at NLO to $$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$8 at NNLO in the nuclear fit, demonstrating that precise electroweak-boson data can make NNLO effects phenomenologically visible (Helenius et al., 2021). KSASG20 combined charged-lepton DIS, neutrino DIS, DY, deuteron data, and new JLab measurements with TMC and HT corrections, obtaining $$f_i^{(A,Z)}(x,Q^2) = \frac{Z}{A} f_i^{p/A}(x,Q^2) + \frac{A-Z}{A} f_i^{n/A}(x,Q^2),$$9 at NLO and $f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$0 at NNLO (Khanpour et al., 2020).

Analysis	Main inputs	Characteristic result
HKN07 (0709.3038)	DIS ratios and DY ratios	Valence well determined; gluon modifications cannot be fixed
nCTEQ15 (Kusina, 2016)	616 DIS, 92 DY, 32 pion points	$f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$1; pion data constrain gluons
TUJU21 (Helenius et al., 2021)	DIS plus 74 $f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$2Pb $f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$3 points	$f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$4 improves from 0.94 to 0.84 at NNLO

5. Universality tests, process dependence, and disputed nuclear corrections

The central theoretical assumption of global nPDF work is universality: the same nuclear PDFs should describe different hard processes once the appropriate coefficient functions are used. Several analyses support this operationally. EPS09 reviewed the accumulated DIS, DY, and hadronic evidence as lending support to QCD factorization in bound nucleons; the 2012 NLO global fit reported no significant overall tension among charged-lepton DIS, neutrino DIS, DY, and $f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$5Au pion data; and EPPS16 found that the addition of charged-current data did not generate significant incompatibility within the fit (Paukkunen et al., 2010, Florian et al., 2012, Zurita, 2018).

At the same time, neutrino–nucleus DIS has long been the principal source of stress tests. A dedicated neutrino–iron analysis extracted iron PDFs directly from NuTeV differential cross sections and dimuon data without applying pre-assumed nuclear corrections, defining for each observable

$f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$6

It found that, except for very high $f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$7, the charged-current neutrino–iron correction factors differ in both shape and magnitude from model-based corrections and from charged-lepton–iron corrections, with a flatter behavior and no reproduction of the classic antishadowing enhancement around $f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$8 seen in charged-lepton data (0806.0723). A subsequent pair of NLO analyses, one based on charged-lepton DIS plus DY and another on neutrino DIS, confirmed a visible discrepancy in the intermediate-$f^{p,n/A}(x,Q^2)=R_f^A(x,Q^2)\, f^{p,n}(x,Q^2),$9 region and obtained only partial compatibility in a compromise fit where the neutrino contribution to $$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$0 was down-weighted by $$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$1 (Kovarik, 2010).

The interpretation of this discrepancy is itself disputed. EPS09 argued that later work suggests the NuTeV inconsistency is an artifact of the NuTeV data rather than a failure of universality (Paukkunen et al., 2010). By contrast, more recent reviews still describe a persistent tension between charged-current and neutral-current DIS and note that only a subset of neutrino data is typically retained in modern global fits (Klasen, 2024). Taken together, these results suggest that universality remains the working framework of nPDF phenomenology, but neutrino-based nuclear corrections continue to be a sensitive methodological problem rather than a fully closed issue.

A related issue arises already for deuterium. Deuteron DIS is often used as a near-free reference in both proton and nuclear analyses, yet several studies stress that Fermi motion, weak binding, and off-shell effects are non-negligible at intermediate and large $$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$2. In the CJ12 proton-PDF analysis with explicit deuteron corrections, three nuclear models—CJ12min, CJ12mid, and CJ12max—led to

$$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$3

showing that nuclear modeling can be comparable in importance to PDF-fit uncertainty for the large-$$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$4 $$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$5-quark (Owens et al., 2012). This is directly relevant to nPDF work because many nuclear ratios are normalized to deuterium and because deuteron corrections propagate into free-proton baselines (Klasen, 2024).

6. Collider phenomenology and future directions

nPDFs are indispensable for heavy-ion and proton–ion phenomenology at RHIC and the LHC. They provide the cold-nuclear-matter baseline needed to disentangle initial-state nuclear effects from final-state medium dynamics in $$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$6Pb and PbPb collisions, and they enter precision predictions for jets, electroweak bosons, prompt photons, heavy flavor, quarkonia, and top production (Tehrani, 2017, Klasen, 2024, Paakkinen et al., 7 Jun 2026).

The LHC has transformed the field. Proton–lead measurements now constrain kinematic regions far beyond fixed-target coverage, reaching down to $$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$7 and up to $$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$8. Among the most important observables are normalized dijets from CMS, which provide strong gluon constraints; isolated photons, which access the gluon through the QCD Compton channel; electroweak bosons, which probe flavor combinations with minimal hadronization ambiguity; and open heavy flavor and quarkonia, which extend gluon sensitivity to very small $$R_i^A(x,Q^2) \equiv \frac{f_i^A(x,Q^2)}{f_i^p(x,Q^2)}.$$9 (Klasen, 2024).

Ultraperipheral collisions have become especially valuable for small-$$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$0 gluons. Recent reviews emphasize coherent $$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$1 photoproduction in PbPb UPCs as direct evidence of gluon shadowing, with

$$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$2

over roughly $$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$3 at $$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$4 (Paakkinen et al., 7 Jun 2026). This channel has become one of the cleanest non-DIS demonstrations of nuclear gluon suppression.

Future facilities are expected to change the precision frontier again. EIC pseudo-data studies based on EPPS16-type re-fits reported that inclusive EIC data reduce the uncertainty in the nuclear gluon ratio by about a factor of 4 at low $$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$5, while charm-tagged pseudo-data improve the high-$$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$6 gluon region by up to a factor of 8 (Zurita, 2018). In the Monte Carlo nNNPDF framework, projected EIC measurements were found to extend direct nPDF constraints down to $$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$7 (Khalek et al., 2019). The LHeC and AFTER@LHC are likewise identified as complementary future programs, and current reviews place special emphasis on the EIC as the decisive facility for distinguishing among competing explanations of nuclear shadowing and for sharply reducing the remaining uncertainty in small-$$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$8 quarks and gluons (Zurita, 2018, Paakkinen et al., 7 Jun 2026).

In contemporary QCD phenomenology, nPDFs therefore occupy a dual role. They are at once a practical ingredient of precision collider calculations and a compact representation of unresolved nuclear many-body dynamics at the partonic level. Their global determination now rests on a broad experimental foundation and on increasingly sophisticated statistical methods, yet the smallest-$$F_2^A(x,Q^2) = \sum_i f_i^{(A,Z)}(x,Q^2)\otimes C_{2,i}(x,Q^2),$$9 gluon, flavor separation in nuclei, deuteron and neutrino corrections, and the microscopic origin of the EMC and shadowing regions remain the principal open problems (Ethier et al., 2020, Klasen, 2024, Paakkinen et al., 7 Jun 2026).