74 2024

# Nuclear PDFs After the First Decade of LHC Data

Michael Klasen1 and Hannu Paukkunen2,3 1Institute for Theoretical Physics,
University of Münster, Münster, Germany; email<EMAIL_ADDRESS>2Department of Physics, University of Jyväskylä, Jyväskylä, Finland; email:
<EMAIL_ADDRESS>3Helsinki Institute of Physics, Helsinki, Finland

###### Abstract

We present a review of the conceptual basis, present knowledge and recent
progress in the field of global analysis of nuclear parton distribution
functions (PDFs). After introducing the theoretical foundations and
methodological approaches for the extraction of nuclear PDFs from experimental
data, we discuss how different measurements in fixed-target and collider
experiments provide increasingly precise constraints on various aspects of
nuclear PDFs, including shadowing, antishadowing, the EMC effect, Fermi
motion, flavor separation, deuteron binding, target-mass and other higher-
twist effects. Particular emphasis is given to measurements carried out in
proton-lead collisions at the Large Hadron Collider, which have revolutionized
the global analysis during the past decade. These measurements include
electroweak-boson, jet, light-hadron, and heavy-flavor observables. Finally,
we outline the expected impact of the future Electron Ion Collider and discuss
the role and interplay of nuclear PDFs with other branches of nuclear,
particle and astroparticle physics.

###### doi:

10.1146/annurev-nucl-102122-022747

###### keywords:

quantum chromodynamics, nuclear structure, parton distribution functions,
collider physics, future developments

††journal: Ann. Rev. Nucl. Part. Sci.

###### Contents

  1. 1 Introduction
  2. 2 Theoretical foundations
    1. 2.1 Factorization in the OPE and QCD-improved parton model
    2. 2.2 Target mass corrections
    3. 2.3 Higher twist corrections
  3. 3 Global analyses of nuclear PDFs
    1. 3.1 The nCTEQ framework
    2. 3.2 The EPPS framework
    3. 3.3 The nNNPDF framework
    4. 3.4 NNLO and model-dependent approaches
    5. 3.5 Overview of global nuclear PDF analyses
  4. 4 Fixed-target data
    1. 4.1 Early DIS data and constraints on quarks
    2. 4.2 DIS at high $x$ and nuclear effects in the deuteron
    3. 4.3 Drell-Yan process and constraints on antiquarks
    4. 4.4 Neutrino DIS data and flavor separation
  5. 5 Collider data
    1. 5.1 Electroweak bosons
    2. 5.2 Photons
    3. 5.3 Light hadrons
    4. 5.4 Jets
    5. 5.5 Heavy quarks and quarkonia
    6. 5.6 Exclusive and inclusive observables in ultraperipheral collisions
  6. 6 Other developments
    1. 6.1 Electron Ion Collider
    2. 6.2 Lattice QCD
    3. 6.3 Relations to other phenomena
  7. 7 Conclusion

##  1 Introduction

The nuclear structure at high energies is an important current research topic,
relevant not only for our understanding of the fundamental quark and gluon
dynamics in protons and neutrons bound in nuclei, but also for elucidating the
formation, properties and evolution of a deconfined state of hadronic matter
that existed in the early Universe, the so-called quark-gluon plasma (QGP).
Nuclear parton distribution functions (PDFs) encode cold binding effects in
nuclei [2, 3], determine the pre-equilibrium phase and initial-state phase
transition to the QGP [4], and are correlated with the precise extraction of
important QGP properties such as its temperature from photons [5] and the
final-state phase transition during chemical freeze-out [6]. While the
evolution of the PDFs with the scale $Q^{2}$, at which they are probed, can be
computed in perturbative Quantum Chromodynamics (QCD) [7, 8, 9, 10], their
dependence on the longitudinal parton momentum fraction inside the hadron is
non-perturbative and must be fitted to experimental data. Traditionally, deep-
inelastic scattering (DIS) of charged leptons or neutrinos and Drell-Yan (DY)
dilepton production with fixed targets have provided the bulk of the data.
However, over the last decade collider data from the Large Hadron Collider
(LHC) at CERN and also from the Relativistic Heavy Ion Collider (RHIC) at BNL
have led to significant improvements in our knowledge of collinear,
unpolarized nuclear PDFs, which are reviewed here.

In the naive parton model, the double-differential charged-lepton DIS cross
section per nucleon

$\frac{d^{2}\sigma^{lA}}{dxdy}=\frac{2\pi\alpha^{2}}{Q^{4}}s\left[1+(1-y)^{2}\right]F_{2}^{A}(x)$
(1)

is directly related to the nuclear structure function
$F_{2}^{A}(x)=x\sum_{q}e_{q}^{2}f_{q}^{A}(x)$ and thus the quark PDFs
$f_{q}^{A}(x)$ in the nucleus $A$, while the gluon PDF and $Q^{2}$-dependence
enter only in the QCD-improved parton model (see Sec. 2). Here, $\alpha$ is
the electromagnetic fine-structure constant, $e_{q}$ the fractional charge of
quark $q$, $Q^{2}$ the virtuality of the exchanged photon, $s=Q^{2}/(xy)$ the
hadronic center-of-mass energy, $y$ the lepton inelasticity, and $x$ the
Bjorken variable. The nuclear structure functions differ from free nucleon
structure functions not only due to an admixture of protons and neutrons, but
also in various other ways depending on the region in $x$. With increasing
$x$, a depletion or shadowing is observed in $F_{2}^{A}(x)$ at $x\lesssim
0.05$, followed by an enhancement or antishadowing at $0.05\lesssim x\lesssim
0.3$, again a depletion – the famous EMC effect – at $0.3\lesssim x\lesssim
0.7$, and ultimately an enhancement due to Fermi motion at $0.7\lesssim x$.
The ratios of large-$A$ isoscalar structure functions to the structure
function of the deuteron (D) – a loosely bound isoscalar state of a proton
($p$) and a neutron ($n$) – can be parameterized and fitted, albeit with
considerable uncertainties, to SLAC [11, 12] and NMC [13] data in the form
[14]

$R(x)=1.10-0.36x-0.28e^{-21.9x}+2.77x^{14.4},$ (2)

i.e. without a dependence on the nuclear mass number $A$ or the hard scale
$Q^{2}$. However, a logarithmic decrease in $A$ has been observed in the EMC
region in the 1984 SLAC data and a logarithmic increase in $Q^{2}$ in the
shadowing region in the more precise NMC data [15]. To first approximation,
this is then also how quark PDFs in a bound nucleon $N$ differ from the free-
nucleon PDFs.

The focus of this review is on model-independent global fits of nuclear PDFs
and the progress due to LHC data over the last decade. To obtain a first
impression of the nuclear dynamics, it is nonetheless interesting to discuss
the main historical experimental measurements and theoretical interpretations
in the different $x$ regions.

Shadowing had been known to be present in real and virtual photon scattering
on nuclei since the 1970s [16, 17]. At the hadron level, it can be interpreted
by assuming that the photon fluctuates from its pointlike state into a
superposition of vector mesons ($\rho$, $\omega$, $\phi$), which then interact
strongly with the nucleons on the surface of the target nucleus (vector meson
dominance). Such nucleons absorb most of the incoming “hadron” flux and
therefore cast a shadow onto the inner ones [18]. At the parton level and in
the nuclear rest frame, the photon can be seen to split into a quark-antiquark
dipole with lifetime $\tau$, that scatters coherently from multiple partons in
the nucleus, if $\tau\geq R_{A}\sim A^{1/3}$ fm or $x\leq 1/(2M_{N}R_{A})\sim
0.1A^{-1/3}$, resulting again in a reduced nuclear cross section [19, 20, 21].

The relative motion of nucleons inside the nucleus was considered first in the
1970s for the deuteron [22] and in the 1980s also for heavier nuclei near the
Fermi surface [23, 24, 25]. The structure function of a nucleon in a nucleus
can then be expressed as a convolution

$\displaystyle
F_{2}^{A}(x_{N})=\int_{x_{N}}^{A}dyf_{A}(y)F_{2}^{N}\left(\frac{x_{N}}{y}\right),$
$\displaystyle\,{\rm where}\,$ $\displaystyle
f_{A}(y)\sim\frac{1}{\sqrt{2\pi\Delta_{A}}}\exp{\left\\{\frac{-[y-(1-\delta_{A})]^{2}}{2\Delta_{A}^{2}}\right\\}},$
(3)

of the nucleon structure function $F_{2}^{N}(x_{N})$ with the nucleon momentum
distribution $f_{A}(y)$. Its peak is shifted away from unity due to soft
nuclear interactions by an amount $\delta_{A}\sim 0.04$, which corresponds
roughly to the ratio of nucleon separation energy over its mass. The width
$\Delta_{A}$ is determined by the Fermi momentum $k_{F}\sim 10$ MeV divided by
the nucleon mass $M_{N}\sim 1$ GeV and thus small. Eq. 3 can therefore be
approximated by a simple rescaling
$F_{2}^{A}(x_{N})=F_{2}^{N}(x_{N}/(1-\delta_{A}))$. The net result of this
rescaling is to deplete the partons in the intermediate $x_{N}$ region,
implying $F_{2}^{A}/F_{2}^{B}<1$ for $A>B$, and to enrich the large $x_{N}\sim
1$ region with $F_{2}^{A}/F_{2}^{B}>1$. The region $x_{N}>1$ can also be
modeled by modifying the Gaussian ansatz for $f_{A}(y)$ in Eq. 3 (e.g. with a
power law tail [26]) and also in deconfinement or cluster models [3], but this
region is usually neglected [27].

The discovery of a suppression of $F_{2}^{\rm Fe}/F_{2}^{\rm D}$ at
$x_{N}=0.65$ of $\sim 0.89$ in 1983 by the EMC collaboration in muon DIS [28]
and its confirmation in reanalyzed iron and aluminum SLAC data from the early
1970s [29, 30] came as a big surprise, since Fermi motion models predicted an
enhancement of $\sim 1.25$ at this value of $x_{N}$ and a suppression only for
$x_{N}<0.5$. It triggered many theoretical explanations at both the nuclear
and the partonic level, and a consensus still has to emerge [31]. Models with
nucleons as the only degrees of freedom in the nucleus must be incomplete,
since the convolution in Eq. 3 violates baryon number and momentum sum rules.
The missing momentum could be carried by pions, whose pair exchanges lead to
an intermediate-range nuclear attraction of $300-500$ MeV, that is canceled by
short-distance vector exchanges of almost equal size. The net effect is an
average binding energy of 8 MeV per nucleon as observed [32]. However, one
would then expect an enhancement of antiquarks and therefore of the DY
process, which has not been seen [33]. The failure of nucleon-only and
nucleon-pion models indicates that the nucleon structure itself is modified by
the medium. This is clearly seen in neutrons, which undergo $\beta$-decay as
free particles, but are stable in nuclei. The parton model interpretation of
the EMC effect is that the medium reduces the number of high-momentum quarks.
This momentum reduction leads, via the uncertainty principle, to the notion
that quarks in nuclei are confined in a larger volume than that of a free
nucleon. There are two general ways to realize this simple idea: either scalar
and vector mean-field effects cause bound nucleons to be larger than free
ones, or short-range correlations (SRCs) cause the nucleon structure to be
modified by including either $NN^{*}$ configurations or deconfined six-quark
configurations that are orthogonal to the two-nucleon wave functions.
Interestingly, two-nucleon SRCs might explain the observed linear correlation
between the magnitude of the EMC effect at $0.3\leq x_{N}\leq 0.7$ and the
size of the plateau observed in quasi-elastic scattering at $1.5\leq x_{N}\leq
2$, which would solve the single-nucleon sum rule problem [34, 35].

A similar compensation mechanism could be at work in the antishadowing region,
which is imposed by shadowing through the momentum sum rule. In the Breit
frame, small momentum quarks and gluons, because of the uncertainty principle,
spread over a distance comparable to the nucleon-nucleon separation. Quarks
and gluons from different nucleons can then overlap spatially and fuse, thus
increasing the density of high momentum partons (antishadowing) at the expense
of that of lower momentum ones (shadowing) [36]. In perturbative QCD, this
process is flavor-dependent and $q\bar{q}\to g$ fusion results, e.g., in
shadowing for antiquarks and antishadowing for gluons. The fact that there is
no clear evidence of antishadowing in the DY process can be interpreted either
with an important role of valence ($v$) quarks or as a consequence of the
evolution in $Q^{2}$ [2, 3].

The non-perturbative nature of nuclear interactions, need for phenomenological
models and incongruous nuclear and partonic interpretations of the effects
described above are strong motivations to parameterize and fit nuclear PDFs to
the available data in a model-independent way. Improving on Eq. 2, Eskola
parameterized in 1992 the ratio of heavier nuclear structure functions over
deuterons separately in each region at the starting scale $Q_{0}^{2}$, matched
it at the transition points (whose definition depended on $A$), and evolved it
in $Q^{2}$ [37]. As observed experimentally [17] and predicted theoretically
[38], shadowing then vanished only very slowly towards larger $Q^{2}$, in
particular for quarks and antiquarks. When the parameterized ratio was fitted
to DIS and $pA$ DY data, while imposing baryon number and momentum
conservation, the nuclear data could be described fairly independently of the
underlying proton PDFs [39]. A rigorous statistical analysis of DIS data
initially led to rather large values of $\chi^{2}$ per degree of freedom (dof)
of $1.82-1.93$ [40], which could, however, be reduced to 1.35 in leading order
(LO) and 1.21 in next-to-leading order (NLO) QCD using more precise DIS and
$pA$ DY data [41, 42]. Collider (RHIC) data on $\pi^{0}$ production from 2006
introduced sensitivity to the gluon density beyond scaling violations with a
resulting $\chi^{2}/$dof $=0.79$ in LO and NLO [43]. An equally good value of
$\chi^{2}/$dof $=0.83$ was obtained in NLO with a similar data set, but with a
direct $A$-dependent parameterization of nuclear PDFs at the starting scale
$Q_{0}^{2}$ [44]. The inclusion of neutrino DIS proved to be more difficult
[45], which triggered a discussion about the universality of nuclear effects
in charged-lepton and neutrino scattering [46, 47, 48, 49, 50]. In the last
decade, a wealth of LHC data on electroweak boson, photon, light and heavy
hadron, and jet production have been taken in proton-lead ($p$Pb) collisions,
which have already had a significant impact on the determination of nuclear
PDFs [51, 52, 53]. These modern developments will be reviewed thoroughly in
the following, whereas other recent reviews have mostly focused on proton PDFs
[54, 55, 56] or sketched a larger multi-dimensional picture of the nucleus
[57].

The remainder of this article is organized as follows: in Sec. 2, we briefly
review the theoretical foundations of nuclear DIS and its factorization. In
Sec. 3, we describe the different methodological approaches in the global fits
of nuclear PDFs. In Secs. 4 and 5, we discuss the impact of the different
experimental data in roughly chronological order, the main focus being of
course on the LHC. The impact of the future Electron Ion Collider and
connections to other fields in nuclear, particle and astroparticle physics
(lattice QCD, the search for gluon saturation, the QGP and astrophysical
phenomena) are briefly addressed in Sec. 6, before we conclude this review in
Sec. 7.

##  2 Theoretical foundations

We start our discussion of nuclear PDFs by reviewing their theoretical
foundations both within the operator product expansion (OPE) and the QCD-
improved parton model, including target-mass and other higher-twist effects.

###  2.1 Factorization in the OPE and QCD-improved parton model

Figure 1: Leading-order diagrams for (a) neutral-current (NC) DIS with charged
leptons, (b) charged-current (CC) DIS with neutrinos or antineutrinos, and (c)
charm dimuon production. The charm quark is understood to hadronize before the
semileptonic decay.

Deep-inelastic scattering of high-energy leptons is the key process for
studying the hadronic structure of nucleons $N$ or nuclei $A$ with mass
$M_{N}$ or $M_{A}=A\,M_{N}$ and four-momentum $P_{N}$ or $P_{A}$ in terms of
their partonic (quark and gluon) degrees of freedom (cf. Fig. 1). The charged
lepton $l$ or (anti-)neutrino $\nu$ has incoming (outgoing) four-momentum $k$
($k^{\prime}$), the squared center-of-mass energy is
$s_{N,A}=(k+P_{N,A})^{2}$, and $X$ represents all final-state hadrons with
total four-momentum $P_{X}$ and squared mass $W_{N,A}^{2}=(q+P_{N,A})^{2}$. In
the laboratory frame, the lepton energy loss is $\nu=q\cdot
P_{N,A}/M_{N,A}=E-E^{\prime}$, and the exchanged vector boson $V$ has squared
momentum transfer $Q^{2}=-q^{2}>0$.

The inclusive differential cross section $d\tilde{\sigma}^{A}\sim
L^{\mu\nu}\tilde{W}^{A}_{\mu\nu}$ can be written as a combination of a
pointlike leptonic tensor $L^{\mu\nu}$ and the hadronic tensor

$\tilde{W}^{A}_{\mu\nu}(P_{A},q)=\frac{1}{4\pi}\int d^{4}z~{}e^{iq\cdot
z}\langle A|J_{\mu}^{\dagger}(z)\ J_{\nu}(0)|A\rangle=\frac{1}{4\pi}\ {\rm
disc}\,\tilde{T}_{\mu\nu}(P_{A},q),$ (4)

where the latter is given in terms of a product of hadronic currents and can
be related to the discontinuity of the virtual forward Compton scattering
amplitude

$\tilde{T}^{A}_{\mu\nu}(P_{A},q)=\int d^{4}z\ e^{iq\cdot z}\ \langle
A|\mathcal{T}J_{\mu}^{\dagger}(z)\ J_{\nu}(0)|A\rangle.$ (5)

The Operator Product Expansion (OPE) then allows one to expand the hadronic
matrix element of the forward scattering amplitude in a complete set of local
operators [58]

$\tilde{T}^{A}_{\mu\nu}(P_{A},q)=-2i\sum_{j,\tau,n}c_{\tau,\mu\nu}^{j,\mu_{1}\cdots\mu_{n}}\langle
A|O^{j,\tau}_{\mu_{1}\cdots\mu_{n}}|A\rangle\
=-2i\sum_{j,k}\frac{2^{2k}}{Q^{4k}}C_{j}^{2k}A_{2k}\tilde{\Pi}^{j,k}_{\mu\nu}\
+\ \mathcal{O}(\tau>2),$ (6)

where $c_{\tau,\mu\nu}^{j,\mu_{1}\cdots\mu_{n}}$ denote the hard scattering,
$\tau$ the twist of the operator $O$, defined as its mass dimension minus its
spin, and $j$ different operators with the same twist. Up to power
corrections, one can identify the product of the perturbative Wilson
coefficients $C_{i}^{2k}$ and reduced hadronic matrix elements $A_{2k}$ as
integer Mellin moments of structure functions

$\int_{0}^{1}dy~{}y^{2k-1}~{}\tilde{F}^{A}_{i}(y,Q^{2})=C_{i}^{2k}A_{2k}+\mathcal{O}(\tau>2)\,,$
(7)

with $y^{2k-1}\to y^{2k-2}$ for $i=2$. The Lorentz structure in terms of
metric tensors and momenta is encoded in $\tilde{\Pi}^{j,k}_{\mu\nu}$.

In the QCD-improved parton model, the nuclear structure functions

${\tilde{F}}^{A}_{i}(x_{A},Q^{2})=\sum_{j=q,g}\int_{x_{A}}^{1}\frac{dy_{A}}{y_{A}}\
C_{i,j}\tilde{f}_{j}^{A}(y_{A},Q^{2})\ +\mathcal{O}(\tau>2)$ (8)

depend on the Bjorken scaling variable $x_{N,A}=Q^{2}/(2q\cdot
P_{N,A})=Q^{2}/(2\nu M_{N,A})$ with $x_{N}\in[0,A]$ ($x_{A}=x_{N}/A\in[0,1]$)
with logarithmic scaling violation in $Q^{2}$ (see below). They are given as
convolutions of target-independent short-distance Wilson coefficients
$C_{i,j}$ with universal nuclear parton distribution functions (PDFs)
$\tilde{f}_{j}^{A}$. Inspection of Eqs. 7 and 8 shows that nuclear PDFs can be
understood as moments of matrix elements of local twist-two operators composed
of quark and gluon fields [59]. The nuclear PDFs
$\tilde{f}_{i}^{A}(x_{A},Q^{2})$ above are related to the more familiar
average-nucleon nuclear PDFs $f_{i}^{A}(x_{N},Q^{2})$ through
$f_{i}^{A}(x_{N},Q^{2})=\tilde{f}_{i}^{A}(x_{A},Q^{2})/A$. This rescaling is a
key step that allows us to compare structure functions across different
nuclei, including the free nucleon. The evolution of the PDFs with the scale
$Q^{2}$ is perturbatively calculable and given by the Dokshitzer-Gribov-
Lipatov-Altarelli-Parisi (DGLAP) evolution equations [7, 8, 9, 10]

$\displaystyle\frac{df_{i}^{A}(x_{N},Q^{2})}{d\ln Q^{2}}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{s}(Q^{2})}{2\pi}\int_{x_{N}}^{A}\frac{dy_{N}}{y_{N}}\
P_{ij}\left(\frac{x_{N}}{y_{N}}\right)\ f_{j}^{A}(y_{N},Q^{2})\,,$ (9)

where $\alpha_{s}$ is the QCD coupling and $P_{ij}$ are the partonic splitting
functions. Furthermore, the PDFs satisfy sum rules due to charge, baryon
number and momentum conservation,

$\displaystyle\int_{0}^{A}dx_{N}\\{f_{u_{v}}^{A},f_{d_{v}}^{A}\\}(x_{N},Q^{2})=\\{2Z+N,Z+2N\\}$
, $\displaystyle\int_{0}^{A}dx_{N}x_{N}\sum_{i}f_{i}^{A}(x_{N},Q^{2})=A,$ (10)

where $Z$ is the electric charge of the nucleus with baryon number $A=Z+N$. It
is therefore common (but not necessary) to decompose the nuclear PDFs as

$f_{i}^{A}(x_{N},Q^{2})=\frac{Z}{A}f_{i}^{p/A}(x_{N},Q^{2})+\frac{A-Z}{A}f_{i}^{n/A}(x_{N},Q^{2})\,,$
(11)

where the bound neutron PDFs $f_{i}^{n/A}(x_{N},Q^{2})$ are commonly obtained
from those of the bound proton $f_{i}^{p/A}(x_{N},Q^{2})$ by assuming isospin
symmetry,

$f_{u,\bar{u}}^{n/A}(x_{N},Q^{2})=f_{d,\bar{d}}^{p/A}(x_{N},Q^{2})\,,\ \
f_{d,\bar{d}}^{n/A}(x_{N},Q^{2})=f_{u,\bar{u}}^{p/A}(x_{N},Q^{2})\,.$ (12)

In principle the above integrations extend to $A$, although the dominant
support of the PDFs is expected to be in the region $x_{A}\leq 1/A$, or
$x_{N}\leq 1$. One therefore usually assumes $f_{i}^{A}(x_{N},Q^{2})=0$ for
$x_{N}>1$, which has the advantage that the same evolution equations can be
used for all nuclei in the interval $x\in[0,1]$.

Accounting for the heavy-quark masses is essential in an accurate description
of the free proton data [60]. Standard methods to handle the quark masses are
nowadays General Mass (GM) Variable Flavor Number Schemes (VFNS) [61] such as
the simplified Aivazis-Collins-Olness-Tung (SACOT) schemes [62, 63, 64] or the
FONLL schemes [65, 66], which provide systematic ways to interpolate between
the Fixed Flavor Number Scheme (FFNS), in which heavy quarks are not
considered as partons, and the Zero Mass (ZM) VFNS, in which heavy quarks are
treated as massless partons.

###  2.2 Target mass corrections

Target mass corrections (TMCs) can be discussed in terms of collinear
factorization [67] or be obtained from the OPE by inverting moments of
structure functions, cf. Eq. 7. They can be written in the general form [68]

$F_{i}^{A,{\rm
TMC}}(x_{N},Q^{2})=\sum_{j}A_{i}^{j}F_{j}^{A}(\xi_{N},Q^{2})+B_{i}^{j}h_{j}^{A}(\xi_{N},Q^{2})+C_{i}g_{2}^{A}(\xi_{N},Q^{2}),$
(13)

where $\xi_{N}\\!=\\!2x_{N}/(1+r_{N})$ is the Nachtmann variable [69] with
$r_{N}\\!=\\!\sqrt{1+4x_{N}^{2}M_{N}^{2}/Q^{2}}$ and

$\displaystyle
h_{j}^{A}(\xi_{N},Q^{2})\sim\int_{\xi_{N}}^{A}d\xi_{N}^{\prime}\frac{F_{j}^{A}(\xi_{N}^{\prime},Q^{2})}{\xi_{N}^{\prime}}$
$\displaystyle\quad,\quad$ $\displaystyle
g_{2}^{A}(\xi_{N},Q^{2})=\int_{\xi_{N}}^{A}d\xi_{N}^{\prime}h_{2}^{A}(\xi_{N}^{\prime},Q^{2})$
(14)

are auxiliary functions with $1/\xi_{N}^{\prime}\to 1/{\xi_{N}^{\prime}}^{2}$
for $j=2$. Specifically, the proportionality factor in $h_{2}^{A}$ is unity
and we have

$F_{2}^{A,{\rm
TMC}}(x_{N},Q^{2})=\left(\frac{x_{N}^{2}}{\xi_{N}^{2}r_{N}^{3}}\right)F_{2}^{A}(\xi_{N},Q^{2})+\left(\frac{6M_{N}^{2}x_{N}^{3}}{Q^{2}r_{N}^{4}}\right)h_{2}^{A}(\xi_{N},Q^{2})+\left(\frac{12M_{N}^{4}x_{N}^{4}}{Q^{4}r_{N}^{5}}\right)g_{2}^{A}(\xi_{N},Q^{2}).$
(15)

Quark masses modify $\xi_{N}$ by $\xi_{N}\rightarrow R_{ij}\xi_{N}$, where the
factor $R_{ij}$ depends on the incoming and outgoing quark masses $m_{i,j}$.
In the case $m_{i}=0$, one obtains the slow-rescaling limit
$R_{ij}=1+(nm_{j})^{2}/Q^{2}$ with $n=1$ for CC [70] and $n=2$ for NC DIS
[60].

###  2.3 Higher twist corrections

Figure 2: Classification of multiple parton scattering in nuclear medium: (a)
interactions internal to the nucleus, (b) initial-state (IS) and (c) final-
state (FS) interactions [71].

Interactions internal to the nucleus as in Fig. 2 (a) change the nuclear PDFs
with respect to those of the free nucleon. However, since only a single parton
participates in the hard scattering, the structure functions can still be
factorized as in Eq. 8, and the leading twist nuclear PDFs can be
parameterized at an initial scale $Q_{0}$, evolved with (in principle
$A$-dependent) evolution equations and fitted to experimental data or modeled
theoretically.

Next-to-leading power corrections to Eq. 8 of ${\cal O}(r_{T}^{2}\sim
1/p_{T}^{2})$, ${\cal O}(m_{J}^{2}/p_{T}^{2})$ and ${\cal
O}(\alpha_{s}(Q^{2})\Lambda^{2}/Q^{2})$ arise from the transverse size $r_{T}$
of the initial nucleus, the non-vanishing invariant mass $m_{J}$ of the final
jet, and IS and FS interactions involving more than one parton as shown in
Figs. 2 (b) and (c). In hadron-nucleus collisions, both are enhanced by
$A^{1/3}$ due to the large density of soft partons in the nucleus. For IS
interactions, $\Lambda^{2}\sim 0.01$ GeV2 is the squared scale of the twist-
four correlation function of single parton pairs and proportional to the
transverse field strength. For FS interactions, long-range soft parton
interactions must also be considered. It can be shown that the $A^{1/3}$
enhancement can be factorized to all powers in hadron-nucleus, but not
nucleus-nucleus collisions and then involves correlation functions of multiple
parton pairs. In general, however, even the hadron-nucleus Drell-Yan cross
section cannot be factorized beyond next-to-leading power [71].

##  3 Global analyses of nuclear PDFs

The inverse problem of extracting nuclear PDFs from experimental data is
approached in a similar way as the determination of (free) proton PDFs [56],
i.e. all global analyses are based on optimizing the correspondence between
theoretical calculations and experimental measurements by minimizing a figure-
of-merit function that is typically of the form

$\chi^{2}=\sum_{i,j}\left(D_{i}-T_{i}\right)C^{-1}_{ij}\left(D_{j}-T_{j}\right)\,.$
(16)

Here, $D_{i}$ denote the experimental values for observables in the fit, and
$T_{i}$ are the corresponding theoretical values, which depend on the PDFs.
The covariance matrix is defined as
$C_{ij}=\sigma_{i}^{2}\delta_{ij}+\sum_{\alpha}\bar{\sigma}_{i\alpha}\bar{\sigma}_{j\alpha}$,
where $\sigma_{i}$ is the total uncorrelated uncertainty added in quadrature
and $\bar{\sigma}_{i\alpha}$ is the correlated systematic uncertainty from
source $\alpha$. How to exactly assign values for $\sigma_{i}$ and
$\bar{\sigma}_{i\alpha}$ varies from one analysis to another and depends on
whether uncertainties are multiplicative or additive [72, 73]. In the case of
fitting nuclear PDFs, only a handful of the data sets provide the full
information on the correlated systematic uncertainties – in most of the cases
only the overall normalization uncertainty is given.

An essential part of global PDF analyses is the propagation of experimental
uncertainties into the PDFs. From the practical point of view, the two
principal methods are the Hessian [74, 75] and Monte Carlo methods [76, 77].
The Hessian uncertainty analysis [74, 75] is based on expanding the $\chi^{2}$
function in the vicinity of its minimum value $\chi^{2}_{0}$,

$\chi^{2}\approx\chi^{2}_{0}+\sum_{i,j}H_{ij}\delta a_{i}\delta
a_{j}=\chi^{2}_{0}+\sum_{i}z_{i}^{2}\,,$ (17)

where $\delta a_{i}$ are deviations from the best-fit parameters and $H_{ij}$
is the second-derivative matrix, the Hessian matrix. In the second step one
diagonalizes the Hessian matrix by finding its eigendirections. The PDF error
sets $S_{k}^{\pm}$ are then defined as deviations along these eigendirections
to positive/negative directions such that $\chi^{2}$ increases by a fixed
amount $\Delta\chi^{2}$. The value for $\Delta\chi^{2}$ can be defined in
various ways. A common feature in the current global analyses is that
$\Delta\chi^{2}$ is of the order of the number of fit parameters for 68%
confidence-level (CL) uncertainties and somewhat higher for 90% CL
uncertainties. In the plots of this review, the uncertainties $\big{(}\delta
X\big{)}^{\pm}$ for a given PDF-dependent quantity $X$ are calculated by the
asymmetric prescription [78],

$\displaystyle\big{(}\delta X\big{)}^{\pm}$
$\displaystyle=\sqrt{\sum_{k}\bigg{[}\begin{subarray}{c}\max\\\
\min\end{subarray}\left[X(S_{k}^{+})-X(S_{0}),X(S_{k}^{-})-X(S_{0}),0\right]\bigg{]}^{2}}\,,$
(18)

where $S_{0}$ denotes the best fit. Out of the fits discussed in this review,
nCTEQ15HQ [51], EPPS21 [52], TUJU21 [79] and KSASG20 [80] make use of the
Hessian method. The Monte Carlo method is based on preparing several fits in
which the central values of the experimental data have been randomly shifted
within the uncertainties. In the case of uncorrelated experimental
uncertainties, the fitted data points are obtained through

$D_{i}\rightarrow D_{i}\left(1+\sigma_{i}R_{i}\right)\,,$ (19)

where $R_{i}$ is a taken from a Gaussian distribution centered around $0$ and
with a unit standard deviation. The above formula can also be generalized to
the case of correlated uncertainties [81, 82]. In the nNNPDF3.0 fit [53], the
100% uncertainty band for a given quantity $X$ is then defined as the
minimum/maximum value obtained by calculating $X$ with all PDF replicas. The
90% uncertainty is defined by disregarding the highest/lowest 10%.

###  3.1 The nCTEQ framework

In the nCTEQ15 NLO analysis [44] and its sequels, the ansatz for the nuclear
PDFs at $Q_{0}=1.3$ GeV follows the CTEQ6M parameterization [83]

$xf_{i}^{p/A}(x,Q_{0}^{2})=c_{0i}x^{c_{1i}}(1-x)^{c_{2i}}e^{c_{3i}x}(1+e^{c_{4i}}x)^{c_{5i}},$
(20)

where $i=u_{v},d_{v},g,\bar{u}+\bar{d},s+\bar{s}$ and we have dropped the
index $N$ of $x$, while

$\frac{f^{p/A}_{\bar{d}}(x,Q_{0}^{2})}{f^{p/A}_{\bar{u}}(x,Q_{0}^{2})}=c^{\prime}_{0}\,x^{c^{\prime}_{1}}(1-x)^{c^{\prime}_{2}}+(1+c^{\prime}_{3}x)(1-x)^{c^{\prime}_{4}}.$
(21)

The normalization coefficients $c_{0i}$ are constrained by the momentum and
valence quark sum rules. The proton baseline is similar to the fit CTEQ6.1M
[84], but has minimal influence from nuclear data [85]. There are currently no
uncertainties associated with this proton baseline PDF. The $A$-dependence of
the parameterization is directly included in the coefficients

$c_{ki}\longrightarrow c_{ki}(A)\equiv p_{ki}+a_{ki}(1-A^{-b_{ki}}),\quad
k=\\{1,...,5\\}.$ (22)

The 16 free parameters in nCTEQ15 describe the $x$-dependence in $u_{v}$,
$d_{v}$, $g$, $\bar{u}+\bar{d}$, while the parameters in $\bar{d}/\bar{u}$
were fixed, as was $s=\bar{s}=\kappa(\bar{u}+\bar{d})/2$. Starting with the
fits nCTEQ15WZ [86] and nCTEQ15WZ+SIH [87] including weak ($W,Z$) boson and
single-inclusive hadron (SIH) production at the LHC, three free parameters
were added for $s+\bar{s}$. The 38 error PDFs are obtained with the Hessian
method and a tolerance of $\Delta\chi^{2}=35$.

In the latest fit nCTEQ15HQ [51], the weakly fragmentation-function dependent
SIH data were complemented by open heavy quark and quarkonium (HQ) production
data from the LHC that had shown great potential to constrain the gluon in
preceding reweighting studies [88, 89]. These data were fitted with a data-
driven method [90], in which the cross sections of hadrons $A_{1}$ and $A_{2}$
are taken to be dominated by the gluon-gluon subprocesses,

$\sigma(A_{1}A_{2}\rightarrow\mathcal{Q}+X)=\int
dx_{1}dx_{2}f_{g}^{A_{1}}\left(x_{1},\mu^{2}\right)f_{g}^{A_{2}}\left(x_{2},\mu^{2}\right)\frac{1}{2\hat{s}}\overline{\left|\mathcal{A}_{gg\rightarrow\mathcal{Q}+X}\right|^{2}}\mathrm{dPS}\,,$
(23)

where
$\mathcal{Q}=D^{0},J/\psi,B\\!\to\\!J/\psi,\Upsilon,\psi^{\prime},B\\!\to\\!\psi^{\prime}$,
the squared factorization scale $\mu^{2}$ is related to the geometric mean of
$M_{\mathcal{Q}}^{2}$ and $p_{T}^{2}$, and $\mathrm{dPS}$ denotes the two-
particle phase space. The effective matrix elements
$\overline{\left|\mathcal{A}_{gg\rightarrow\mathcal{Q}+X}\right|^{2}}$ are
parameterized by a generalized Crystal Ball function,

$\displaystyle\begin{split}\overline{\left|\mathcal{A}_{gg\rightarrow\mathcal{Q}+X}\right|^{2}}=\frac{\lambda^{2}\kappa\hat{s}}{M_{\mathcal{Q}}^{2}}e^{a|y|}&\times\begin{cases}e^{-\kappa\frac{p_{\rm
T}^{2}}{M_{\mathcal{Q}}^{2}}}&\text{ if }p_{\rm T}\leq\left\langle p_{\rm
T}\right\rangle\\\ e^{-\kappa\frac{\left\langle p_{\rm
T}\right\rangle^{2}}{M_{\mathcal{Q}}^{2}}}\left(1+\frac{\kappa}{n}\frac{p_{\rm
T}^{2}-\left\langle p_{\rm
T}\right\rangle^{2}}{M_{\mathcal{Q}}^{2}}\right)^{-n}&\text{ if }p_{\rm
T}>\left\langle p_{\rm T}\right\rangle\end{cases},\end{split}$ (24)

where $M_{\mathcal{Q}}$ denotes the mass of particle $\mathcal{Q}$,
$\hat{s}=x_{1}x_{2}s$, and $p_{\rm T}$ and $y$ correspond to the transverse
momentum and rapidity of $\mathcal{Q}$. The free parameters $\lambda$,
$\kappa$, $\left<p_{T}\right>$, $n$ and $a$ are fitted for each final state to
$pp$ data. The fits agree with NLO GM-VFNS [91] and NRQCD [92] calculations
within their scale uncertainties. In the fit nCTEQ15HQ, the LHC $p$Pb data are
consistently included in the fit as absolute cross sections in the case of all
observables.

###  3.2 The EPPS framework

The latest EPPS analysis, EPPS21 [52], is rooted in a series of global fits
[43, 93, 94, 95, 96], which parameterize the bound proton PDFs at the starting
scale $Q_{0}=1.3$ GeV as

$f^{p/A}_{i}(x,Q_{0}^{2})=R^{p/A}_{i}(x,Q_{0}^{2})f^{p}_{i}(x,Q_{0}^{2})\,.$
(25)

Here, $i=u_{v},d_{v},g,\bar{u},\bar{d},s=\bar{s}$, and the free proton PDFs
$f^{p}_{i}(x,Q^{2})$ are taken from the fit CT18A [97], which includes more
LHC $pp$ data sensitive to strange quarks than the default CT18 fit. This
reduces the dependence of the proton baseline on the $\nu$Fe DIS data, which
is also part of the CT18A analysis. The nuclear modifications
$R^{p/A}_{i}(x,Q_{0}^{2})$ are parameterized through 24 free parameters. The
parameterization is piecewise smooth in $x$, so that parameters controlling
different $x$ regions mix as little as possible, i.e.

$\displaystyle
R_{i}^{p/A}(x,Q^{2}_{0})=\left\\{\begin{array}[]{lr}a_{0i}+a_{1i}\big{(}x-x_{ai}\big{)}\Big{[}e^{-xa_{2i}/x_{ai}}-e^{-a_{2i}}\Big{]}\,,&x\leq
x_{ai}\\\\[5.0pt]
b_{0i}x^{b_{1i}}\big{(}1-x\big{)}^{b_{2i}}e^{xb_{3i}}\,,&x_{ai}\leq x\leq
x_{ei}\\\\[5.0pt]
c_{0i}+c_{1i}\left(c_{2i}-x\right)\left(1-x\right)^{-\beta_{i}}\,,&x_{ei}\leq
x\leq 1\,,\end{array}\right.$ (29)

where $x_{ai}$ and $x_{ei}$ are the locations of the anticipated antishadowing
maximum and minimum of the EMC effect, respectively. The $A$-dependence is
encoded in such a way that larger nuclei tend to have larger nuclear effects
at $x=0,x_{ai},x_{ei}$ through

$\displaystyle R_{i}^{p/A}(x,Q^{2}_{0})=1+\Big{[}R_{i}^{p/A_{\rm
ref}}(x,Q^{2}_{0})-1\Big{]}\left(\frac{A}{A_{\rm ref}}\right)^{\gamma_{i}}\,,\
\gamma_{i}>0\,,\ A_{\rm ref}=12\,.$ (30)

However, for very small nuclei a monotonic $A$ scaling is not necessarily a
justified assumption (e.g. certain small nuclei are more tightly bound), and
such deviations are also allowed in the EPPS parameterization through
$R^{p/A}_{i}(x,Q_{0}^{2})\longrightarrow
1+f_{A}\Big{[}R^{p/A}_{i}(x,Q_{0}^{2})-1\Big{]}$, where $f_{A}=1$ by default.
The parameterization is applied for $A\geq 3$, while for smaller nuclei the
nuclear modification is set to unity. This is in line with the baseline CT18A
proton PDFs [97] which include DIS data on deuteron targets with no nuclear
corrections [98]. The parameterization of $R^{p/A}_{i}(x,Q_{0}^{2})$, and
thereby the nuclear PDFs, are not restricted to be positive at small $x$. Even
if the positivity was imposed at $Q_{0}$, the backward evolution to smaller
$Q^{2}$ would result in negative values, especially for gluons. Requiring the
positivity at $Q_{0}$ appears thus too restrictive and would induce an
increased dependence on the parameterization scale $Q_{0}$.

The EPPS analyses use ratios of cross sections or structure functions whenever
possible in order to remove the dependence on free proton PDFs as much as
possible. Ratios of cross sections have the additional advantage that they are
perturbatively much more stable than absolute cross sections, reducing the
risk of fitting missing higher order effects into the nuclear modifications.
Furthermore, experimental uncertainties – known and unknown – can be expected
to cancel, e.g. the one from the luminosity. The uncertainties in the EPPS21
analysis are evaluated through the Hessian method with a global tolerance of
$\Delta\chi^{2}=33$. In addition, the dependence of nuclear modifications on
the free proton PDFs is mapped by repeating the fit with each of the 58 CT18A
error sets as well. As a result, the EPPS21 fit comes with 106 error sets,
which are correlated with the CT18A error sets.

###  3.3 The nNNPDF framework

In the nNNPDF3.0 NLO analysis [53], six independent combinations of nuclear
PDFs are parameterized at $Q_{0}=1\,{\rm GeV}$ in an evolution basis,

$\displaystyle xf_{i}^{p/A}(x,Q_{0}^{2})$ $\displaystyle=$ $\displaystyle
B_{i}x^{\alpha_{i}}(1-x)^{\beta_{i}}{\rm NN}_{i}(x,A)\,,\
i=\Sigma,T_{3},T_{8},V,V_{3},g\,,$ (31)

where $\Sigma,T_{3},T_{8},V,V_{3}$ label certain combinations of quark PDFs
[54], $g$ is the gluon and ${\rm NN}_{f}(x,A)$ represents the value of the
neuron in the output layer of the neural network associated to each
independent nuclear PDF. The normalisation coefficients
$B_{\Sigma}=B_{T_{3}}=B_{T_{8}}=1$, while $B_{V}$, $B_{V_{3}}$ and $B_{g}$
enforce the momentum and valence sum rules and are determined at $Q_{0}$ for
each value of $A$. The preprocessing exponents $\alpha_{i}$ and $\beta_{i}$
are required to control the small and large $x$ behavior of the nuclear PDFs.
They are fitted simultaneously with the network parameters. The exponents
$\alpha_{V}$ and $\alpha_{V_{3}}$ are restricted to lie in the range $[0,5]$
during the fit to ensure integrability of the valence distributions. The other
exponents $\alpha_{i}$ are restricted to the range $[-1,5]$, consistent with
momentum sum rule requirements, while the exponents $\beta_{i}$ lie in the
range $[1,10]$. The figure of merit is defined as

$\chi_{\rm fit}^{2}=\chi_{\rm t_{0}}^{2}+\kappa_{\rm pos}^{2}+\kappa_{\rm
BC}^{2}\,,$ (32)

where the first term is the contribution from experimental data with a
covariance matrix that takes into account the normalization uncertainties of
the different data sets with a self-consistent iterative $t_{0}$ procedure
[73]. The second term $\kappa_{\rm pos}^{2}$ imposes the positivity of
physical cross sections, and the third term $\kappa_{\rm BC}^{2}$ ensures that
in the limit $A\to 1$ the nNNPDF predictions reduce to those of the free-
proton boundary condition using a grid of 100 points, half of which are
distributed logarithmically between $x=10^{-6}$ and $0.1$ and the remaining
half are linearly distributed between $0.1$ and $0.7$. Ratios of open heavy
quark production (i.e. $D^{0}$ meson) data are included with the reweighting
method.

The free proton baseline PDFs are fitted separately with the NNPDF3.1
methodology [99] to all data sets included in the NNPDF4.0 NLO analysis [100]
except those involving nuclei with $A\geq$ 2\. Each replica of nNNPDF3.0 has
been fitted using a randomly chosen replica of proton PDFs coming from these
fits. As a result, the uncertainties of nNNPDF3.0 reduce to those of the
proton baseline in the $A\to 1$ limit, and the correlations between free-
proton and nuclear PDFs are charted and available for the users. The $A\geq$ 2
data that was removed from the proton fit are included in the fit of nuclear
PDFs. Consequently, the nuclear effects in deuteron are fitted in a model-
independent way.

A common challenge in training neural network models is the choice of the
hyperparameters such as their architecture and activation functions, the
optimization algorithm and learning rates. Here, Kernel Density Estimators
[101, 102] have been proven to outperform random or grid searches in selecting
the most promising sets. One observes, e.g., that a network with one hidden
layer and $25$ nodes beats a network with two hidden layers.

###  3.4 NNLO and model-dependent approaches

A few analyses have been performed in next-to-next-to-leading order (NNLO)
QCD, albeit with a restricted data set. For example, the nNNPDF1.0 analysis
[103] fitted only the light-quark singlet, octet and gluon PDFs to the
charged-lepton DIS data. The TUJU19 analysis [104] took a step further by
complementing the charged-lepton DIS data with neutrino DIS data in the
framework of the open source tool xFitter [105]. Using an ansatz

$xf^{p/A}_{i}\left(x,Q_{0}^{2}\right)=c_{0i}\,x^{c_{1i}}(1-x)^{c_{2i}}\left(1+c_{3i}\,x+c_{4i}\,x^{2}\right)$
(33)

and a definition of $\chi^{2}$ that followed the HERAPDF2.0 analysis [106],
the free proton baseline was first fitted to the HERA data. Due to the limited
data set, $s=\bar{s}=\bar{u}=\bar{d}$ had to be assumed, but with 13 free
parameters and a tolerance of $\Delta\chi^{2}=20$, a very similar proton
baseline as HERAPDF2.0 was obtained. These proton PDFs were then utilized in
the fit to the heavy nuclei data, in which the $A$-dependence was directly
encoded in the fit parameters $c_{k}$ as in the nCTEQ fits [44], including the
case of the deuteron. In the updated TUJU21 analysis [79], also LHC data for
weak boson production were included both in the proton and nuclear fits.
However, it was still not possible to consider the flavor decomposition in the
sea quark sector. The TUJU19 and TUJU21 analyses were the first to consider
both the proton and heavier nuclei in the same framework, even though they did
not provide a common Hessian matrix.

The KSASG20 analysis [80] is based on CT18 [97] free-proton PDFs and uses a
nuclear modification factor of cubic HKN form [40, 42]

${\cal
W}_{i}(x,A)=1+\left(1-\frac{1}{A^{\alpha}}\right)\frac{a_{i}(A)+b_{i}(A)\,x+c_{i}(A)\,x^{2}+d_{i}(A)\,x^{3}}{(1-x)^{\beta_{i}}},$
(34)

which is flexible enough to accomodate both shadowing and antishadowing
effects. From nuclear volume and surface contributions, one has $\alpha=1/3$
[107], the $a_{i}$ control shadowing, and $\beta_{i}$ are related to Fermi
motion. Only charged-lepton and neutrino DIS data, supplemented with fixed-
target DY data, were fitted with a very restricted flavor decomposition.

Model-dependent nuclear PDFs have been proposed based on the leading twist
approximation for nuclear shadowing with antishadowing constrained by the
momentum sum rule [20], SRCs of nucleon pairs motivated by the EMC effect
[108], and the four-component Kulagin-Petti model [109, 110], which can
reproduce a variety of DIS, fixed-target DY [111] and LHC weak-boson data
[112]. As the viewpoint of the present review is a data-based global analysis
of nuclear PDFs, we will not discuss the modeling of nuclear effects further.

###  3.5 Overview of global nuclear PDF analyses

The key features of the nuclear PDF analyses described above are summarized in
Tab. 1. In addition to the main methodological assumptions and input
parameters, we also list the fitted experimental data types, grouped into
fixed-target and collider data. Since TUJU21 and KSASG20 include (almost) no
collider data, we focus on nCTEQ15HQ, EPPS21 and nNNPDF3.0 in the following.
Figure 3 shows the regions of $x$ and $Q^{2}$ covered by the world data and
included in these analyses. As can be appreciated from the plot, the LHC data
taken during the first decade of $pA$ runs have radically expanded the
available range both in $x$ and $Q^{2}$ as well as diversified the global
analysis. With this large coverage, the question of process independence of
nuclear PDFs can now be addressed much more convincingly.

Table 1: Key features of recent global analyses of nuclear PDFs.

Analysis | nCTEQ15HQ [51] | EPPS21 [52] | nNNPDF3.0 [53] | TUJU21 [79] | KSASG20 [80]
---|---|---|---|---|---
Theoretical input: |  |  |  |  |
Perturbative order | NLO | NLO | NLO | NNLO | NNLO
Heavy-quark scheme | SACOT$-\chi$ | SACOT$-\chi$ | FONLL | FONLL | FONLL
Value of $\alpha_{s}(M_{Z})$ | 0.118 | 0.118 | 0.118 | 0.118 | 0.118
Charm mass $m_{c}$ | $1.3\,{\rm GeV}$ | $1.3\,{\rm GeV}$ | $1.51\,{\rm GeV}$ | $1.43\,{\rm GeV}$ | $1.3\,{\rm GeV}$
Bottom mass $m_{b}$ | $4.5\,{\rm GeV}$ | $4.75\,{\rm GeV}$ | $4.92\,{\rm GeV}$ | $4.5\,{\rm GeV}$ | $4.75\,{\rm GeV}$
Input scale $Q_{0}$ | 1.3 GeV | 1.3 GeV | 1.0 GeV | 1.3 GeV | 1.3 GeV
Data points | 1484 | 2077 | 2188 | 2410 | 4353
Independent flavors | 5 | 6 | 6 | 4 | 3
Parameterization | Analytic | Analytic | Neural network | Analytic | Analytic
Free parameters | 19 | 24 | 256 | 16 | 18
Error analysis | Hessian | Hessian | Monte Carlo | Hessian | Hessian
Tolerance | $\Delta\chi^{2}=35$ | $\Delta\chi^{2}=33$ | N/A | $\Delta\chi^{2}=50$ | $\Delta\chi^{2}=20$
Proton PDF |  $\sim$CTEQ6.1 | CT18A |  $\sim$NNPDF4.0 |  $\sim$HERAPDF2.0 | CT18
Proton PDF correlations |  | ✓ | ✓ |  |
Deuteron corrections | (✓)a,b | ✓c | ✓ | ✓ | ✓
Fixed-target data: |  |  |  |  |
SLAC/EMC/NMC NC DIS | ✓ | ✓ | ✓ | ✓ | ✓
– Cut on $Q^{2}$ | 4 GeV2 | 1.69 GeV2 | 3.5 GeV2 | 3.5 GeV2 | 1.2 GeV2
– Cut on $W^{2}$ | 12.25 GeV2 | 3.24 GeV2 | 12.5 GeV2 | 12.0 GeV2 |
JLab NC DIS | (✓)a | ✓ |  |  | ✓
CHORUS/CDHSW CC DIS | (✓/-)b | ✓/- | ✓/- | ✓/✓ | ✓/✓
NuTeV/CCFR $2\mu$ CC DIS | (✓/✓)b |  | ✓/- |  |
$pA$ DY | ✓ | ✓ | ✓ |  | ✓
$\pi A$ DY |  | ✓ |  |  |
Collider data: |  |  |  |  |
$Z$ bosons | ✓ | ✓ | ✓ | ✓ |
$W^{\pm}$ bosons | ✓ | ✓ | ✓ | ✓ |
Light hadrons | ✓ | ✓d |  |  |
– Cut on $p_{T}$ | 3 GeV | 3 GeV |  |  |
Jets |  | ✓ | ✓ |  |
Prompt photons |  |  | ✓ |  |
Prompt D0 | ✓ | ✓ | ✓e |  |
– Cut on $p_{T}$ | 3 GeV | 3 GeV | 0 GeV |  |
Quarkonia ($J/\psi$, $\psi^{\prime}$, $\Upsilon$) | ✓ |  |  |  |

a nCTEQ15HIX [27]; b nCTEQ15$\nu$ [113]; c through CT18A; d only $\pi^{0}$ in
DAu; e only forward ($y>0$).

Figure 3: An illustration of the $x$ and $Q^{2}$ regions probed by the current
lepton-$A$, pion-$A$ and proton-$A$ data included in the global analyses of
nuclear PDFs. Figure 4: Comparison of the 208Pb nuclear modifications
resulting from the EPPS21 (full, blue) [52], nCTEQ15HQ (dashed, red) [51] and
nNNPDF3.0 (dot-dashed, green) [53] global analyses of nuclear PDFs, i.e. the
PDFs of lead divided by the summed PDFs of 82 free protons and 126 free
neutrons. Uncertainty bands correspond to 90% CL.

Figure 4 compares the nuclear modifications of the lead nucleus PDFs at
$Q^{2}=10$ GeV2 from EPPS21 (full, blue) [52], nCTEQ15HQ (dashed, red) [51]
and nNNPDF3.0 (dot-dashed, green) [53]. Qualitatively, there is good overall
agreement between all three within the 90% CL uncertainty bands (shaded
areas). Closer inspection reveals nevertheless still significant differences
both among the central values and the widths of the uncertainty bands in
several distributions and $x$ regions. Due to the precise fixed-target NC DIS
data, the nuclear modifications of $u$ and $d$ quarks (and to a lesser extent
antiquarks) are well constrained at $x\sim 0.1$, though the widths of the
error bands differ by at least a factor of two. Larger uncertainties remain in
the Fermi motion region, in particular in nNNPDF3.0, where the $x$ dependence
at the starting scale is not parametric. In the case of $\bar{u}$ and
$\bar{d}$, the fit nCTEQ15HQ comes with the smallest uncertainties, which
could, however, be due to not fitting $\bar{u}$ and $\bar{d}$ separately.
Thanks to the LHC data, the gluon uncertainties are now much better
constrained than in the previous rounds of global fits down to $x\sim
10^{-5}$, and they also impact the sea (anti)quarks and their uncertainties at
higher $Q^{2}$. Below $x\sim 10^{-3}$, the nCTEQ15HQ and nNNPDF3.0 gluon
uncertainty bands do not overlap, which will be reflected in some of the plots
comparing theoretical predictions and LHC data in Sec. 5. The largest
uncertainties are seen for the strange quark distributions, which are
constrained only by – to some extent problematic – neutrino data and by LHC
weak boson data, where the strange quark originates, however, mostly from
gluon splittings. In Supplemental Material we provide also a comparison of the
absolute nuclear PDFs.

##  4 Fixed-target data

We now turn to the detailed discussion of the available experimental data and
their impact on global nuclear PDF analyses in roughly chronological order. We
first focus on the early fixed-target NC and CC DIS and DY data, but review
also the more recent JLab DIS data.

###  4.1 Early DIS data and constraints on quarks

Measurements of fixed-target electron and muon NC DIS on various nuclei from
SLAC, EMC and NMC, but also from the BCDMS, FNAL and HERMES experiments form
the backbone of global nuclear PDF determinations [3]. Most of these data sets
can be fitted with an excellent $\chi^{2}/$dof, but a few of them such as the
1988 EMC measurement of $F_{2}^{\rm Sn}/F_{2}^{\rm D}$ [114] are difficult to
describe [42, 51, 53, 79, 80, 115]. In contrast, the NMC $F_{2}^{\rm
Sn}/F_{2}^{\rm C}$ data [15, 116] can be fitted well [51, 52, 53, 79, 80].
Other outliers are e.g. the E665 data [117] on $F_{2}^{\rm C}/F_{2}^{\rm D}$,
$F_{2}^{\rm Ca}/F_{2}^{\rm D}$ and $F_{2}^{\rm Pb}/F_{2}^{\rm D}$, though e.g.
the ratio $(F_{2}^{\rm Pb}/F_{2}^{\rm D})/(F_{2}^{\rm C}/F_{2}^{\rm
D})=F_{2}^{\rm Pb}/F_{2}^{\rm C}$ is consistent with the NMC data [116].

The kinematic reach $x\gtrsim 5\times 10^{-3}$ and $Q^{2}\lesssim 140\,{\rm
GeV}^{2}$ of these fixed-target data is naturally more limited than the one of
the HERA experiments H1 and ZEUS, which form the bulk of the data in free
proton analyses [55, 56]. In addition, cuts are often applied in order to
limit effects of target-mass and other higher-twist corrections, which could
be larger in nuclear reactions [71]. At low $Q^{2}$, NC DIS is governed by
virtual photon exchange. In the kinematic region of fixed-target experiments,
the cross section

$\frac{d^{2}\sigma^{lA}}{dxdQ^{2}}=\frac{4\pi\alpha^{2}}{Q^{4}}\left[F^{A}_{2}(x,Q^{2})\left(\frac{y^{2}}{2}+1-y-\frac{xyM^{2}}{s-M^{2}}\right)-xy^{2}F^{A}_{\rm
L}(x,Q^{2})\right]$ (35)

is dominated by the structure function $F^{A}_{2}(x,Q^{2})$, which in LO is
sensitive only to the squared charge-weighted sum of quarks and antiquarks
(cf. Sec. 1). As the $lA$ DIS data that enter the global fits are given in
terms of ratios

$\frac{d\sigma^{lA_{1}}}{d\sigma^{lA_{2}}}\approx\frac{F_{2}^{A_{1}}}{F_{2}^{A_{2}}}\,,$
(36)

where $A_{2}$ is typically D or carbon (C), the DIS data can mainly directly
constrain the overall nuclear modification of valence and sea quarks. The
contributions of gluons enter the cross section only at order $\alpha_{s}$,
and the direct constraints for the gluon densities are therefore weak.
However, the gluons drive the $Q^{2}$ dependence of $F_{2}^{A}(x,Q^{2})$ at
small values of $x$ [118],

$\frac{dF^{A}_{2}(x,Q^{2})}{d\log
Q^{2}}\approx\frac{10\alpha_{s}(Q^{2})}{27\pi}xf^{A}_{g}(2x,Q^{2})\,,\
x\rightarrow 0\,.$ (37)

Through this relation, it was understood early on that the $Q^{2}$ dependence
of the ratios $F_{2}^{\rm Sn}/F_{2}^{\rm C}$ measured by NMC [15] around
$x\approx 0.01\ldots 0.02$ and $Q^{2}\approx 1\ldots 10\,{\rm GeV}^{2}$ can
constrain the $A$ dependence of the gluon nuclear modifications [119]. In
particular, a very strong $A$ dependence of gluons would contradict the
measured positive $Q^{2}$ slopes of $F_{2}^{\rm Sn}/F_{2}^{\rm C}$. There are
also similar HERMES data for $F_{2}^{\rm Kr}/F_{2}^{\rm D}$ [120], but the
$Q^{2}$ lever arm is not as long in the perturbative regime. In principle, the
longitudinal structure function $F^{A}_{L}$ carries a direct sensitivity to
the gluon [121], but data are scarce [122]. Likewise, the cross sections for
charm production would provide more direct information on the gluons – modulo
a possible intrinsic charm PDF [123] – but not much data are available [124,
125].

Since $F_{2}^{A}$ probes predominantly the squared charge-weighted sum of
quark PDFs, the flavor decomposition is also difficult to pin down. To
understand this, let us write e.g. the valence quark distributions as

$\displaystyle f_{u_{v}}^{A}$
$\displaystyle=R^{A}_{v}\left(\frac{Z}{A}f_{u_{v}}^{p}+\frac{A-Z}{A}f_{d_{v}}^{p}\right)+\delta
R^{A}_{v}\left(\frac{2Z}{A}-1\right)\frac{f_{u_{v}}^{p}f_{d_{v}}^{p}}{f_{u_{v}}^{p}+f_{d_{v}}^{p}}\,,$
(38) $\displaystyle f_{d_{v}}^{A}$
$\displaystyle=R^{A}_{v}\left(\frac{Z}{A}f_{d_{v}}^{p}+\frac{A-Z}{A}f_{u_{v}}^{p}\right)-\delta
R^{A}_{v}\left(\frac{2Z}{A}-1\right)\frac{f_{u_{v}}^{p}f_{d_{v}}^{p}}{f_{u_{v}}^{p}+f_{d_{v}}^{p}}\,,$
(39)

where
$R^{A}_{v}\equiv({R^{p/A}_{u_{v}}f_{u_{v}}^{p}+R^{p/A}_{d_{v}}d_{d_{v}}^{p}})/({f_{u_{v}}^{p}+f_{d_{v}}^{p}})$
is an average nuclear modification of the valence quarks and the difference is
$\delta R^{A}_{v}\equiv R^{p/A}_{u_{v}}-R^{p/A}_{d_{v}}$. The first terms
dominate and thus the large-$x$ data constrain very tightly the average
modification $R^{A}_{v}$. Having data for several different combinations of
$Z$ and $A$, combined with the fact that the nuclear effects are expected to
scale with $A$, will give also constraints on $\delta R^{A}_{v}$ and thereby
to the mutual differences of nuclear effects in up and down valence quarks.
The same reasoning naturally applies for the up and down sea quarks. Eqs. 38
and 39 also clearly demonstrate the anticorrelation of the nuclear effects
between up and down quarks. For an isoscalar nucleus $Z=A/2$, the up and down
quark distributions are always equal.

To facilitate the interpretation of nuclear effects, many early NC DIS
experiments corrected their data for isospin effects, i.e. for the unequal
numbers of protons and neutrons in heavy nuclei compared to the deuteron.
These data were long taken at face value – and are still done so e.g. in
nNNPDF3.0 – and fitted by setting $Z=N=A/2$. However, this is not necessary in
global fits and has in the past even caused some confusion about the nuclear
valence quark modification [96]. Comparing
$F_{2}^{A}=[ZF_{2}^{p/A}+(A-Z)F_{2}^{n/A}]/A$ (cf. Eq. 11) with the isoscalar
expression $\hat{F}_{2}^{A}=[F_{2}^{p/A}+F_{2}^{n/A}]/2$ leads to
$\hat{F}_{2}^{A}=\beta F_{2}^{A}$ with

$\beta=\frac{A}{2}\left(1+\frac{F_{2}^{n/A}}{F_{2}^{p/A}}\right)/\left(Z+(A-Z)\frac{F_{2}^{n/A}}{F_{2}^{p/A}}\right)\,.$
(40)

The experiments then assumed $F_{2}^{n/A}/F_{2}^{p/A}=F_{2}^{n}/F_{2}^{p}$ and
parameterized this ratio from DIS data on protons and deuterons with, e.g.,
$1-0.8x$ [12] or $0.92-0.86x$ [126]. These functions then allow to calculate
$\beta$ and either apply it also to the theoretical calculations [104] or to
remove the isoscalar correction from the data and fit the true nuclei [27,
96].

###  4.2 DIS at high $x$ and nuclear effects in the deuteron

Recent precise JLab measurements taken with $6-10$ GeV electrons on various
nuclear targets [127, 128, 129, 130] at low to intermediate $Q^{2}$ and
$W^{2}$ can considerably reduce the nuclear PDF uncertainties in the high $x$
region [27, 131]. This region requires, however, a good control over potential
target-mass (cf. Sec. 2.2) and other higher-twist corrections (cf. Sec. 2.3),
hadronic resonances and (similarly to the other fixed-target DIS data) nuclear
effects in the deuteron. Target-mass corrections scale with powers of
$M_{N}^{2}/Q^{2}$ and are thus suppressed even for heavy nuclei. Their leading
effect is a shift in the probed momentum fraction $x_{N}$ to the Nachtmann
variable $\xi_{N}$ [68]. When the kinematic cuts are relaxed from $Q^{2}>4$
GeV2, $W^{2}>12.25$ GeV2 to $Q^{2}>1.69$ GeV2, $W^{2}>2.89$ GeV2, the
subleading target mass effects provide a uniform shift of below one percent
for all nuclei, leaving the ratios $F_{2}^{A}/F_{2}^{\rm D}$ unaffected [27].
Some global analyses therefore include the JLab data with lower kinematic cuts
[27, 52, 80, 131]. The compatibility of the JLab data with global fits can be
taken as evidence that other higher-twist effects can be neglected [52]. They
can, however, also be parameterized as

$\displaystyle F_{2}^{A}(x,Q)$ $\displaystyle\to
F_{2}^{A}(x,Q)\left[1+\frac{A^{1/3}h_{0}x^{h_{1}}(1+h_{2}x)}{Q^{2}}\right]\,,$
(41)

where the values $\\{h_{0},h_{1},h_{2}\\}=\\{-3.3~{}{\rm GeV}^{2},1.9,-2.1\\}$
come from the CJ15 proton PDF analysis [132] and $A^{1/3}$ scaling is assumed
(cf. Sec. 2.3). This form of higher-twist corrections leads to a slight
reduction of $F^{A}_{2}$ at intermediate $x\,{\sim}\,0.3$ and a substantial
enhancement at high $x$ and $Q^{2}<16\,{\rm GeV}^{2}$, and it improves the
global $\chi^{2}$ by around 3% within the nCTEQ15HIX global analysis [27].

In the resonance region at very low $W^{2}\in[1.21;2.89]$ GeV2, the nuclear
effects at large $x$ are surprisingly similar to those in DIS, which may
signal the applicability of quark hadron duality due to the averaging over
nuclear resonances [133]. As described in the introduction, Fermi motion can
be accounted for by a convolution of the nucleon structure function with the
nucleon momentum distribution or effectively by a rescaling of the variable
$x$. The rise of $F_{2}^{A}/F_{2}^{\rm D}$ at large $x$ can be well described
by the parameterization [27]

$x^{\prime}=x-\varepsilon\,x^{\kappa}\log_{10}A\ .$ (42)

As the PDFs decrease with $x$, the negative shift ensures that the transformed
function is larger than the unmodified one and non-vanishing as $x\to 1$. The
overall size of the rescaling effect is controlled by $\varepsilon$,
$\kappa>0$ ensures that only the large-$x$ region is modified, and the
$\log_{10}A$ term implies an increasing modification across the full range of
nuclear $A$ values from the proton $(A=1)$ to lead $(A=208)$. A good
description of the JLab data is obtained with $\kappa=10$ and $\epsilon\sim
0.03$. While this suggests that it may be possible to expand the kinematic
reach to $W^{2}<2.89$ GeV2, the resonance region is currently avoided in all
global fits (cf. Tab. 1).

Given that the nuclear DIS data are usually presented as ratios
$F_{2}^{A}/F_{2}^{\rm D}$ and that most global fits of proton PDFs use
deuteron data as well, a good control over the nuclear effects in the deuteron
is required. The deuteron is much more loosely bound than heavier nuclei and
therefore often approximated as an isoscalar combination of a free proton and
neutron. However, its structure at large $x$ is still modified by Fermi
motion, nuclear binding and off-shell effects, while at small $x$ rescattering
still induces some shadowing. These nuclear effects are of the order of a few
percent in the available DIS data and below 1% in the available DY data. They
can be accounted for in different ways: one option is to use a free-proton
baseline that fitted deuteron DIS data without nuclear corrections [52]. In
this case the deuteron nuclear effects are, to some extent, fitted into the
$u$ vs. $d$ quark flavor separation. Then no additional nuclear correction
should be applied, although also this has been done, and in this case the
effect of double counting should be quantified [80]. Alternatively, one can
rescale the fitted $F_{2}^{A}/F_{2}^{\rm D}$ data by a ratio $F_{2}^{\rm
D}/F_{2}^{p}$ [27] as modeled, e.g., within the CJ15 global proton analysis
[132] (cf. also [134]). If this is done with a free proton PDF that already
includes deuteron data [27], the same deuteron correction should be applied in
both cases [85]. A third possibility is to fit the free proton without
deuteron data and then the deuteron in the same way as the other nuclei [79].
The theoretical treatment of the deuteron affects the description of all NC
DIS data and therefore also the question of the compatibility of CC DIS data
with NC DIS and electroweak boson production at the LHC (cf. Secs. 4.4 and
5.1).

###  4.3 Drell-Yan process and constraints on antiquarks

The Drell-Yan process, i.e. the inclusive production of electroweak gauge
bosons in hadron collisions, followed by a leptonic decay of the gauge boson,
has been of enormous historical importance for the quark flavor separation in
protons [54]. For heavier nuclei, fixed-target measurements have been made by
the FNAL E605 [135], E772 [33] and E866 [136] experiments in $pA$ collisions
covering several nuclei from carbon to tungsten in the kinematic range
$x>10^{-2}$ and dilepton mass $M_{ll}<15$ GeV. At such low values of $M_{ll}$,
well below the $Z$-boson peak, the DY process is dominated by an off-shell
intermediate photon with the cross section differential in the lepton pair
rapidity $y_{ll}$ and $M_{ll}$ given by

$\displaystyle\frac{d^{2}\sigma^{pA}_{\rm
DY}}{dM_{ll}dy_{ll}}\sim\sum_{q}e_{q}^{2}\left[f^{p}_{q}(x_{1})f^{A}_{\bar{q}}(x_{2})+f^{p}_{\bar{q}}(x_{1})f^{A}_{q}(x_{2})\right]$
$\displaystyle{\rm with}$ $\displaystyle x_{1,2}=\frac{M_{ll}e^{\pm
y_{ll}}}{\sqrt{s}}\,,$ (43)

which tests a squared charge-weighted combination of quarks and antiquarks.
Fixed-target experiments generally have larger acceptance in the $x_{1}\gg
x_{2}$ region, where $x_{1}$ is defined with respect to the proton beam, and
in this case the first term in Eq. 43 is the dominant one, where
$f^{p}_{q}(x_{1})$ is mostly determined by the valence quark content of the
proton. For isoscalar nuclei, $f^{A}_{\bar{u}}=f^{A}_{\bar{d}}$, and the cross
section ratios between $pA$ and $p$D collisions, measured in the E772 and E866
experiments, become

$\frac{d\sigma_{\rm DY}^{pA}}{d\sigma_{\rm DY}^{p{\rm D}}}_{\Big{|}{\rm
isoscalar}\ A}\approx\frac{f^{A}_{\bar{u}}(x_{2})}{f^{\rm
D}_{\bar{u}}(x_{2})}=\frac{f^{A}_{\bar{d}}(x_{2})}{f^{\rm
D}_{\bar{d}}(x_{2})}\,.$ (44)

As a result, the measured DY ratios are sensitive to the nuclear modifications
of sea-quark distributions at $x_{2}\sim 0.03...0.3$ [137, 138]. In principle,
the dependence on $M_{ll}$ should also retain sensitivity to the gluon through
the DGLAP evolution. Combining the photon-mediated DIS and DY measurements
gave historically the first chance to disentangle the nuclear effects in
valence and sea quarks leading to the conclusion that there was not such a
clear antishadowing effect in sea quarks as there was for valence quarks. The
E605 data [135] for $p$Cu collisions is given in terms of absolute cross
sections and is often used in fits of proton PDFs. In the future, fixed-target
$pA$ data from the FNAL E906/SeaQuest experiment [139] are expected to improve
the precision of the available data. The renewed facilities at RHIC should
also be able to provide new measurements on the DY process [140], and similar
measurements are planned at the LHCb experiment at the LHC [141].

In principle, the pion-nucleus DY process has the potential to constrain the
flavor decomposition of the valence quarks [142]. It depends on the pion PDFs,
but this dependence cancels largely in ratios of nuclear cross sections.
Unfortunately it turns out that the precision of the $\pi A$ DY data from the
CERN NA3 [143], NA10 [144] and FNAL E615 [145] experiments is not high enough
to provide significant discrimination power on top of the DIS data. However,
the new CERN-based facility AMBER [146] may be able to improve upon the
current precision.

###  4.4 Neutrino DIS data and flavor separation

Due to the weak nature of neutrino interactions, heavy nuclear targets such as
iron or lead have traditionally been used to obtain CC DIS data with
sufficient statistics. Despite the nuclear targets, these data have routinely
been included in global analyses of proton PDFs, with or without nuclear
corrections [147], and have formed the principal constraint for a possible $s$
vs. $\bar{s}$ asymmetry [148]. In addition, determinations of the weak mixing
angle in neutrino DIS [149] have relied on a sufficient understanding of the
nuclear structure [150, 151, 152]. In the case of CC neutrino DIS, the
differential cross section is

$\frac{d^{2}\sigma^{\nu,\bar{\nu}A}}{dxdy}\\!=\\!\frac{G_{F}^{2}M_{W}^{4}}{2\pi
xyQ^{2}}\left(\\!\frac{Q^{2}}{Q^{2}\\!+\\!M_{W}^{2}}\\!\right)^{2}\\!\left[\\!\left(1\\!-\\!y\\!-\\!\frac{x^{2}y^{2}M^{2}_{N}}{Q^{2}}\right)\\!F^{\nu,\bar{\nu}A}_{2}\\!\\!+\\!y^{2}xF^{\nu,\bar{\nu}A}_{1}\pm\\!\left(y\\!-\\!\frac{y^{2}}{2}\right)xF^{\nu,\bar{\nu}A}_{3}\right]\\!,$
(45)

where $G_{F}$ is the Fermi constant, $M_{W}$ the $W$-boson mass and the $+$
($-$) sign is taken for incoming (anti-)neutrinos. At LO and high $Q^{2}$,

$\displaystyle d^{2}\sigma^{\nu A}$
$\displaystyle\propto\left(f_{d}^{A}+f_{s}^{A}+f_{b}^{A}\right)+(1-y)^{2}\left(f_{\bar{u}}^{A}+f_{\bar{c}}^{A}\right)\,,$
(46) $\displaystyle d^{2}\sigma^{\bar{\nu}A}$
$\displaystyle\propto\left(f_{\bar{d}}^{A}+f_{\bar{s}}^{A}+f_{\bar{b}}^{A}\right)+(1-y)^{2}\left(f_{u}^{A}+f_{c}^{A}\right)\,.$
(47)

Due to the suppressing factor $(1-y)^{2}$, there is an increased sensitivity
to the strange quark distribution in comparison to NC charged-lepton DIS,
particularly for antineutrinos. In addition, the up and down quarks enter the
cross sections with different weights than in the case of NC charged-lepton
DIS. Adding the neutrino data thus helps in constraining the differences
between nuclear effects in up and down quarks (cf. Eqs. 38 and 39). However,
in a global fit these data are also sensitive to the assumed proton PDFs.

Data on inclusive neutrino DIS have been taken by e.g. by CDHSW [153], CCFR
[154] and NuTeV [14] on iron and by CHORUS on lead [155]. Also the charm
production has been measured through muonic decays of produced charmed hadrons
(cf. Fig. 1) [156, 157]. In principle, all of these data should be relevant
for nuclear PDFs, but are only partially included in global analyses due to
concerns about possible mutual tensions between the neutrino data sets,
tensions with the charged-lepton DIS data, and due to the fact that part of
these neutrino data are in some cases already used in the proton PDF fits that
are used as baselines in the fits of nuclear PDFs. One of the difficulties is
also that there are no references from $\nu p$ or $\nu$D scattering, so that
the data are reported as absolute cross sections.

Figure 5: Average ratios of neutrino (left) and antineutrino (right) cross
sections as measured by CHORUS [155], NuTeV [14] and CDHSW [153] to a
theoretical NLO prediction with CT18A PDFs within the kinematic range
$Q^{2}>4$ GeV2 and $W^{2}>12.25\,{\rm GeV}^{2}$. The data are compared with
EPPS21 [52], nCTEQ15HQ [113] and nNNPDF3.0 [53] predictions. The SLAC/NMC NC
DIS parameterization in Eq. 2 is also shown as a reference.

Figure 5 shows neutrino (left) and antineutrino (right) cross sections divided
by the theoretical NLO predictions with CT18A proton PDFs in the SACOT-$\chi$
scheme, including approximate target mass and electroweak corrections as used
in Ref. [46]. The ratios are evaluated as weighted averages over
$Q^{2}>4\,{\rm GeV}^{2}$ and $W^{2}>12.25\,{\rm GeV}^{2}$ as in Ref. [113].
The obtained ratios are compared with predictions for NuTeV data using EPPS21
[52], nCTEQ15HQ [113] and nNNPDF3.0 [53]. Also the SLAC/NMC NC DIS
parameterization (Eq. 2) is shown for comparison. Tensions between different
data sets, nuclear PDFs, and data and theory are clearly visible. The largest
differences between nuclear PDFs occur in the case of $\bar{\nu}$ DIS at
$x\gtrsim 0.2$, where the nNNPDF3.0 values are significantly above those of
EPPS21 or nCTEQ15HQ. This can be explained by the large enhancement of
$\bar{d}$ and $s$ densities of nNNPDF3.0 in comparison to EPPS21 or nCTEQ15HQ
(cf. Supplemental Material). From the data sets, in particular the NuTeV
neutrino data stand out from the others, but also larger deviations between
the CHORUS and CDHSW data have been observed, if less restrictive kinematic
cuts are imposed and electroweak corrections are neglected [113]. To some
degree the observed tensions can be alleviated by normalizing the data by the
cross sections integrated over $x$ and $y$ [46, 48], by neglecting the NuTeV
systematic error correlations, or by introducing additional cuts in $x$ [113].
The current consensus seems to be that at least the CHORUS data can be
included in global analyses without significant tensions. In addition, the
charm dimuon data are used in nCTEQ15$\nu$ [113] and nNNPDF3.0 [53] and the
CDHSW data in TUJU21 [79] and KSASG20 [80]. There have been speculations about
differences in nuclear shadowing in CC and NC processes [21], at least at low
$Q^{2}$ [20], but $W^{\pm}$ and $Z$ production at the LHC probing nuclear PDFs
at significantly higher $Q^{2}$ can be fitted well in global analyses.

In the future, novel neutrino-nucleus DIS data may become available through
dedicated experiments measuring neutrinos produced in high-luminosity $pp$
collisions at the LHC. Indeed, the first observations of such collider
neutrinos have already been made by the FASER [158] and SND@LHC [159]
collaborations. The impact of such future measurements on nuclear PDFs has
been considered recently in Ref. [160].

##  5 Collider data

We now turn to the discussion of the LHC ($p$Pb), but also RHIC (DAu)
measurements used in global fits of nuclear PDFs. Within collinear
factorization, the $p$A (or DA) cross sections

$d\sigma(pA\rightarrow\mathcal{O}+X)=\sum_{\rm i,j[,k]}f_{i}^{p}\otimes
f_{j}^{A}\otimes d\hat{\sigma}(ij\rightarrow\mathcal{O},[k]+X)\
[\,\otimes\,D_{k}^{\mathcal{O}}\,]$ (48)

for the observable $\mathcal{O}$ involve convolutions of (nuclear) PDFs
$f_{i,j}^{p,A}$ with perturbative partonic cross sections $d\hat{\sigma}$ and
in the case of inclusive light or heavy-flavored hadron ($h$) production non-
perturbative fragmentation functions (FFs) $D_{k}^{h}$. These are numerically
costly and must in many cases be evaluated with precomputed grids [161]. Since
one of the colliding objects is a proton and also the nuclear PDFs depend on
the proton, absolute LHC cross sections depend directly and indirectly on the
proton PDFs. To reduce this dependence and also cancel other theoretical and
experimental uncertainties, the nuclear modification ratio

$\displaystyle R_{pA}$ $\displaystyle\equiv
d\sigma(pA\rightarrow\mathcal{O}+X)\big{/}d\sigma(pp\rightarrow\mathcal{O}+X)$
(49)

and forward-to-backward ratio

$\displaystyle R_{\rm FB}$ $\displaystyle\equiv
d\sigma(pA\rightarrow\mathcal{O}+X)_{\big{|}y>0}\,\big{/}d\sigma(pA\rightarrow\mathcal{O}+X)_{\big{|}y<0}$
(50)

are often introduced. Here $y$ refers to the rapidity of the observable
$\mathcal{O}$.

In fits of collider data, normalization uncertainties originating from the
luminosity measurements play a special role. For example, the measured and
calculated $R_{p{\rm Pb}}$ for hadron production in the $y\gg 0$ region (small
$x_{N}$) at the LHC is often rather flat, and changes in the nuclear PDFs can
be compensated by treating the normalization uncertainty as a correlated
systematic uncertainty, see e.g. Ref. [113]. If the luminosity uncertainty is
common for $y>0$ and $y<0$, the $R_{\rm FB}$ is free from this additional
freedom. Today, all global fits of nuclear PDFs account for the systematic
normalization uncertainties.

###  5.1 Electroweak bosons

Table 2: Summary of $Z$, $W^{\pm}$, and low invariant mass $Z/\gamma^{*}$
rapidity distributions available from LHC $p$Pb collisions.

Data set | nCTEQ15HQ [51] | EPPS21 [52] | nNNPDF3.0 [53] | TUJU21 [79] | KP16 [112]
---|---|---|---|---|---
Run-I: |  |  |  |  |
ATLAS $Z$ [162] | ✓ | ✓ | ✓ | ✓ | ✓
CMS $Z$ [163] | ✓ | ✓ | ✓ | ✓ | ✓
ALICE $Z$ [164] |  |  | ✓b |  |
LHCb $Z$ [165] | ✓ |  | ✓b |  |
ATLAS $W^{\pm}$ [166]c | ✓ |  |  |  | ✓
CMS $W^{\pm}$ [167] | ✓ | ✓ | ✓ |  |
ALICE $W^{\pm}$ [164] | ✓ |  | ✓b |  |
Run-II: |  |  |  |  |
CMS $Z$ [168] |  |  | ✓b |  |
CMS $Z/\gamma^{*}$ [168] |  |  | ✓b |  |
ALICE $Z$ [169] |  |  | ✓b |  |
LHCb $Z$ [170] |  |  |  |  |
CMS $W^{\pm}$ [171, 172] | ✓ | ✓a | ✓ | ✓ |
ALICE $W^{\pm}$ [173] |  |  |  |  |

a added in EPPS21 [52]; b added in nNNPDF3.0 [53]; c unpublished.

The first global analysis of nuclear PDFs to include $W$ and $Z$ boson data
from $p$Pb collisions was EPPS16 [96]. However, the impact of the Run-I ATLAS
[162] and CMS [163, 167] data was still rather limited due to low statistics.
Thereafter, measurements for these electroweak processes have been published
by all four LHC experiments at $\sqrt{s}=5.02$ TeV (Run-I) and by ALICE, CMS
and LHCb at $\sqrt{s}=8.16$ TeV (Run-II). They are now, in different
combinations, used in all recent global analyses [51, 52, 53] and also in the
NNLO and model-dependent fits TUJU21 [79] and KP16 [112]. The available data
are summarized in Tab. 2.

At the moment, the most stringent constraints come from the Run-II CMS $W$
data [172]. Fig. 6 compares the nuclear modification ratio $R_{p{\rm Pb}}$
constructed from the CMS Run-II [171, 172] measurements with NLO calculations
using the EPPS21 [52], nCTEQ15HQ [51] and nNNPDF3.0 [53] nuclear PDFs. While
the EPPS21 analysis included the shown $R_{p{\rm Pb}}$ data, the nCTEQ15HQ and
nNNPDF3.0 analyses fitted absolute $p$Pb cross sections. As one can see, the
spread between the different predictions is still rather significant. In
comparison to a calculation with no nuclear effects, i.e. 82 free protons and
126 free neutrons, the data indicates a clear sign of shadowing at forward
rapidities or $x\ll 1$. The relative ordering of EPPS21, nCTEQ15HQ and
nNNPDF3.0 values follows the one of the corresponding gluon shadowing in Fig.
4. Also, as can be seen from Fig. 4, even after inclusion of these electroweak
data, the overall variation in the strange quark PDF is still quite
significant, which indicates that the constraints for the nuclear strange
quark PDFs are still not very strong.

Figure 6: Nuclear modification ratios for $W^{+}$ (left) and $W^{-}$ bosons
(right) at CMS Run-II [171, 172] compared with EPPS21 [52], nCTEQ15HQ [51],
nNNPDF3.0 [53] and a calculation with 82 free protons and 126 free neutrons.

The CMS Run-II measurement for $Z$ boson production [168] reports similarly
small uncertainties to the $W$ boson measurement. However, it is not possible
to obtain a good quantitative description of these data with any nuclear PDFs
due to large fluctuations of the data around midrapidity ($y_{ll}=0$), which
lead e.g. to a $R_{\rm FB}$ that does not tend to unity towards
$y_{ll}\rightarrow 0$ as one would expect. Along with the on-shell $Z$
production, CMS measured also low-mass cross sections in the window $15\,{\rm
GeV}<M_{ll}<60$ GeV. Within the TUJU21 analysis [79] it was noticed that to
simultaneously reproduce the normalization of the CMS low-mass and $Z$ cross
sections, the NNLO QCD corrections appear to be necessary. This is the first
time the necessity of NNLO corrections is seen in the case of $p$Pb
collisions.

Dielectron pairs have also been measured by ALICE at Run-I in the low-mass
region $M_{ll}<3$ GeV and with $0<p_{T,ll}<8$ GeV [174], which is in principle
very sensitive to the gluon density and avoids the fragmentation contribution
present for real photons [175, 176, 177]. Currently the data are unfortunately
still dominated by the heavy-flavor ($c$, $b$) decay background, but the
statistics should be improved in Run-II and the background reducible with
heavy-flavor tagging, in particular in LHCb [141].

###  5.2 Photons

Another electroweak probe of nuclear PDFs is the prompt production of real
photons with finite transverse momentum $p_{T}$ [178, 179, 180, 181]. It
proceeds directly through quark-antiquark annihilation ($q\bar{q}\to g\gamma$)
and the QCD Compton process ($qg\to q\gamma$), which dominates at large
$p_{T}$ [182]. The radiation of massless photons from quarks in pure QCD
processes gives rise to a photon fragmentation contribution [183, 184, 185],
which is important at small $p_{T}$. Isolating the photon and thus reducing
the surrounding hadronic energy (e.g. in a cone) suppresses the fragmentation
component and also the non-prompt photon background from hadronic (in
particular pion) decays [186].

Direct photon production in proton (and pion) nucleus collisions was first
measured in fixed-target mode at the Fermilab E706 experiment [187]. Also
PHENIX [188] and STAR [189] measurements in DAu collisions at RHIC are
available. At the LHC, high-$p_{T}$ isolated photons in $p$Pb [190] collisions
have been measured by the ATLAS collaboration. Ratios of these ATLAS cross
section measurements to the corresponding $pp$ data [191] have a reduced
sensitivity to missing higher-order effects, fragmentation functions and
proton PDFs. They can be reasonably described at NLO QCD and are therefore
included in the nNNPDF3.0 [53] analysis. The absolute $pp$ and $p$Pb cross
sections are, however, underestimated by NLO QCD by up to 30% at the lowest
values of $p_{\rm T}\sim 20\,{\rm GeV}$. This could indicate the necessity to
include NNLO corrections [192]. The impact of the ATLAS prompt photon data in
the global fit is small in comparison to dijet and heavy flavor production due
to larger uncertainties. Looking into the future, the ALICE collaboration has
proposed to build a new forward calorimeter (FoCal) for LHC Run-IV [193, 194],
optimized for direct photons in the rapidity region $3.2<\eta<5.8$.

###  5.3 Light hadrons

As discussed in Sec. 4, before the beginning of the LHC era not much was known
about the nuclear gluon PDFs. The first direct evidence for the presence of
shadowing, antishadowing and the EMC effect in gluons came from inclusive
hadron production in D$A$ collisions at RHIC. This process involves a gluon
contribution already at LO and is therefore a candidate to constrain the
nuclear gluons in the perturbative region, i.e. when the transverse momentum
$p_{T}$ of the hadron is sufficiently large. This possibility was first
discussed in Refs. [94, 41, 195] in the light of early RHIC data [196, 197,
198, 199]. The first global analysis to fit this type of data was EPS08 [95],
which included negatively-charged hadron BRAHMS [197] as well as PHENIX [198,
200] and STAR [201] pion data. It was, however, noticed that the rapidity
dependence of the BRAHMS data at low values of $p_{T}\gtrsim 2\,{\rm GeV}$ was
too strong to be optimally reproduced within a global fit, inducing tensions
with the NMC data for $F_{2}^{\rm Sn}/F_{2}^{\rm C}$ [15]. These negatively-
charged hadron BRAHMS data were eventually dropped from the EPS09 analysis
[43], when it was noticed that it was difficult to reproduce even the $pp$
reference data. The later EP(P)S analyses [52, 96] have retained only PHENIX
$\pi^{0}$ data [200], while nCTEQ15 [44] included these and updated STAR
$\pi^{0}$ data [189].

Inclusive hadron production is not only sensitive PDFs, but also to the final-
state hadronization encoded in the parton-to-hadron FFs. In the EP(P)S and
nCTEQ fits, the FFs are taken from global fits of hadron production in
$e^{+}e^{-}$, $eN$, and $pp$ collisions [202, 203, 204, 205, 206]. The
sensitivity to FFs was recently studied within the nCTEQ15WZ+SIH [87]
analysis, which also propagated the FF uncertainties, when available, into the
fit. Combining the latest RHIC DAu data on pions, kaons and $\eta$ mesons
[189, 200, 201, 207] with the corresponding ALICE $p$Pb measurements [208,
209, 210] led to a consistent description of the data and to a considerable
reduction in the nuclear gluon uncertainty. The fact that consistent global
fits down to $p_{T}>3\,{\rm GeV}$ are possible can be also taken as an
indication that higher-twist final-state rescattering (cf. Sec. 2.3) is indeed
a subleading effect. The inclusive hadron data have therefore also been
retained in the latest nCTEQ15HQ analysis [51]. Alternatively, the RHIC data
in DAu collisions have also been interpreted in terms of nuclear-modified FFs
[211, 212], which were used in the DSSZ [115] global analysis of nuclear PDFs
and resulted in reduced nuclear effects for the gluon PDF.

The latest LHC measurements of $R_{p{\rm Pb}}$ for high-$p_{T}$ neutral pion
production come from LHCb [213], complementing the ALICE midrapidity data with
forward/backward measurements. While the forward-rapidity data agree with the
NLO predictions with nuclear PDFs, there appears to be a slight normalization
difference between the predictions and the LHCb data at negative rapidities.
The preliminary LHCb data for $\eta$ mesons looks consistent with nuclear PDFs
[214]. These data are not yet included in global fits. The LHCb forward pion
data also agree well with the LHCb forward charged-hadron ($h^{\pm}$) data
[215]. The corresponding backward data are, however, in disagreement with the
nuclear PDF predictions, which hints that the baryon production in the lead-
going direction cannot be described solely within the factorization. The same
issue is visible at midrapidity as well [216, 217, 218, 219] and is more
pronounced at RHIC [198, 199]. Moreover, it has been noticed that even in
simpler $pp$ collisions at LHC energies the collinear factorization around
midrapidity appears to be applicable only at $p_{\rm T}\gtrsim 10\,{\rm GeV}$
for $h^{\pm}$ production [220]. As a result, only the production of mesons is
considered in global fits of nuclear PDFs.

###  5.4 Jets

Jet measurements in $p$Pb collisions probe the intermediate- to large-$x$
regime of nuclear PDFs at large interaction scales ($Q^{2}\gtrsim 10^{3}$
GeV2). A complication in $p$Pb compared to $pp$ collisions is the
significantly larger background from the underlying event. Indeed, in Glauber-
type models an average $p$Pb collision contains around $7\pm 5$ $pN$
interactions [221]. To reduce the model dependence in how the multi-parton
interactions (MPIs) are dealt with, the jet $p_{T}$ must therefore be large
enough, or the jet cone must be small enough, to reduce the probability of
particles from MPIs to occupy the same phase space. However, at too small cone
sizes the jet cross sections become unstable due to an incomplete cancellation
of infrared divergences. Also the hadronization corrections, which tend to
widen the partonic jets, grow. While in $pp$ collisions all non-perturbative
corrections are applied to the theoretical predictions, the $p$Pb data have
already been subtracted for the backgrounds from MPIs. This works rather well,
i.e. the obtained ratios $R_{p{\rm Pb}}$ are broadly consistent with the
expectations from nuclear PDFs.

Figure 7: Comparison of the normalized CMS dijet nuclear modification ratios
[222] with NLO calculations using EPPS21, nCTEQ15HQ and nNNPDF3.0 nuclear
PDFs.

Currently, the most constraining data are the Run-I CMS dijet data
differential in the average $p_{T}$ ($p^{\rm ave}_{\rm T}$) and rapidity
($\eta_{\rm dijet}$) of the two jets [222]. They supersede the earlier dijet
data [223], which were included in the EPPS16 analysis. The cross sections are
normalized to the rapidity-integrated cross section, so that most of the
systematic uncertainties cancel. The resulting spectra in $pp$ are then so
precise that they challenge the theoretical description, the NLO perturbative
QCD calculations being in tension with the data [224]. The impact of the NNLO
corrections [225] is still unclear. Despite this tension, ratios of normalized
cross sections between $p$Pb and $pp$ collisions, $R_{\rm pPb}^{\rm norm.}$,
are broadly consistent with nuclear PDFs and included in the EPPS21 and
nNNPDF3.0 global fits, where they have a large impact on the gluon. However,
the most forward data points at the edge of the detector acceptance indicate a
suppression, which cannot be fitted. This questions the reliability of
extracting nuclear PDFs from these data.

Part of the CMS dijet data are compared in Fig. 7 with NLO calculations using
EPPS21, nCTEQ15HQ and nNNPDF3.0. The depletion in the backward direction
$\eta_{\rm dijet}\ll 0$ indicates a presence of an EMC effect for gluons
(though it is partially obscured due to the contribution of valence quarks),
whereas the depletion in the forward direction $\eta_{\rm dijet}\gg 0$ is
consistent with the expected gluon shadowing. At $\eta_{\rm dijet}\gg 0$, the
less shadowed gluons of nCTEQ15HQ seen in Fig. 4 lead to the observed higher
prediction for $R_{\rm pPb}^{\rm norm.}$. The nCTEQ15HQ error bands are the
widest, as these data were not included in the fit. The EPPS21 and nNNPDF3.0
uncertainties are smaller, and near $\eta_{\rm dijet}\sim 1$ there is even a
discrepancy among them. Data for single inclusive jets are also available
[226, 227, 228], but the uncertainties are clearly larger than in the dijet
measurements.

###  5.5 Heavy quarks and quarkonia

Table 3: Heavy quark production data available from LHC $p$Pb collisions.

Observable ${\cal O}$ | $D^{0}$ | $J/\psi$ | $\Upsilon(1S)$ | $\psi(2S)$ | B0, B± |  $c$ jet |  $b$ jet
---|---|---|---|---|---|---|---
Run-I: |  |  |  |  |  |  |
ATLAS |  |  [241, 242]a |  [242]a |  [242]a |  |  |
CMS |  |  [243]a | [244] |  [245]a |  | [246] | [247]
ALICE |  [248, 249, 250]a |  [251, 252]a, [253] | [254] |  [255]a |  |  | [256]
LHCb |  [257]a,b,c |  [258]a | [259] |  |  |  |
Run-II: |  |  |  |  |  |  |
ALICE |  |  [260]a, [261] |  [262]a |  [263]a |  |  |
LHCb | [264] |  [265]a |  [266]a |  | [267] |  |
Fixed target: |  |  |  |  |  |  |
LHCb | [268, 269] | [268, 270] |  | [270] |  |  |

a included in nCTEQ15HQ [51]; b included in EPPS21 [52]; c included in
nNNPDF3.0 [53].

The possibilities of constraining the gluon PDF with inclusive heavy-flavor
production at the LHC have been actively investigated in $pp$ [229, 230, 231,
232] as well as in $p$Pb [51, 88, 89, 233, 234] collisions. However, the
theoretical approaches vary from one analysis to another, and this is also the
case in the fits of nuclear PDFs that include heavy-quark data. The EPPS group
uses GM-VFNS calculations [91, 235, 236], in which heavy quarks are active
partons above the mass thresholds, resumming collinear logarithms from the
initial- and final-state radiation and thereby matching with the variable-
flavor structure of the nuclear PDFs. The nNNPDF group employs an approach
which supplements fixed-order calculations with a similar, though less
complete resummation of collinear logarithms through parton showers [237, 238,
239]. The nCTEQ15HQ group relies on effective matrix-element fitting [240]
introduced in Sec. 3.1. It uses the fact that in fixed-order calculations the
$gg$ initial state dominates at low $p_{T}$ and allows fitting both open heavy
quark and quarkonium production including the hadronization process, which for
quarkonia remains to be fully understood [271].

The four LHC collaborations have collected a vast data set on $D^{0}$,
$B^{0}$, $B^{\pm}$, $J/\psi$, $\Upsilon$ and $\psi^{\prime}$ mesons (cf. Tab.
3), which allow to extend the range in $x_{N}$ to below $10^{-5}$, i.e. more
than one (two) order(s) of magnitude lower than LHC electroweak boson (jet)
production at scales from $m_{c}^{2}$ to $10^{3}$ GeV2 (cf. Fig. 3). Including
these data even partially, the gluon uncertainties of nCTEQ15HQ, EPPS21, and
nNNPDF3.0 have shrunk considerably below $x_{N}=10^{-2}$ in comparison to
their respective predecessors nCTEQ15WZ+SIH [87], EPPS16 [96], and nNNPDF2.0
[272]. While not included in the current global fits, the CMS collaboration
has also measured $c$ [246] and $b$ jets [247], ALICE $b$ jets [256], LHCb
inclusive $B$-meson [267] production discussed e.g. in Ref. [236], and ALICE
heavy-flavor decay electrons [273, 274]. First heavy-flavor measurements have
also been carried out by LHCb in the fixed-target mode with different nuclei
(He, Ar, Ne) [268, 269, 270]. This may eventually allow to study the
$A$-dependence of nuclear PDFs.

Figure 8:  Left: The LHCb Run-II $R_{\rm FB}$ [264] compared to NLO GM-VFNS
calculations [235] using EPPS21, nCTEQ15HQ and nNNPDF3.0 nuclear PDFs. Right:
The LHC Run-II data [277, 278, 279, 280] for exclusive $J/\psi$ production in
PbPb collisions compared to NLO calculations using EPPS21, nCTEQ15HQ and
nNNPDF3.0 nuclear PDFs. The factorization scales $\mu\sim 2.2$ GeV have been
chosen to match the ALICE data at $y=0$.

An important data set in the current global fits is the LHCb Run-I $D^{0}$
measurement [257], which is included in all three fits. In the forward
direction ($y\gg 0$, small $x$), the nuclear modification ratio $R_{p{\rm
Pb}}$ shows a clear suppression consistent with shadowing. In the backward
direction ($y\ll 0$, larger $x$) at the intersection between shadowing and
antishadowing, $R_{p{\rm Pb}}$ is closer to unity. This behavior is consistent
with the CMS dijet and $W^{\pm}$ data. The ALICE $D$-meson data [250] lie at
midrapidity in between the LHCb acceptance and have a somewhat different
normalization. The recent LHCb Run-II $D^{0}$ data [264] are consistent with
nuclear-PDF predictions in the forward direction (shadowing), but indicate a
stronger suppression than expected in the backward direction. Given that these
$R_{p{\rm Pb}}$ data use a $pp$ reference interpolated between $5\,{\rm TeV}$
and $13\,{\rm TeV}$, $R_{\rm FB}$ could arguably be more accurate. Figure 8
(left) compares the new LHCb Run-II measurement with the predictions obtained
using EPPS21, nCTEQ15HQ and nNNPDF3.0 PDFs in a NLO GM-VFNS calculation [235].
Despite the fact that all three use the $5\,{\rm TeV}$ $p$Pb $D^{0}$ data as
an input, there are still significant differences among the predictions.
Recently, preliminary LHCb Run-I data on the $R_{p{\rm Pb}}$ of $D^{+}$ and
$D^{+}_{s}$ have also appeared [275]. They are consistent with the $D^{0}$
results at $y\gg 0$, but the $D^{+}$ data deviate from the $D^{0}$ results at
$y\ll 0$.

The prospects of using top quark production in $p$Pb and PbPb collisions to
understand nuclear PDFs were first quantitatively discussed in Ref. [276].
While the large mass of the top quark renders the production cross sections
small in comparison to charm or beauty production, the process was predicted
to be visible at the LHC. Total top quark cross sections have thereafter been
measured by CMS [281] and ATLAS [282] and also in PbPb collisions by CMS
[283]. The ATLAS measurement in $p$Pb is consistent with the nCTEQ15HQ, EPPS21
and TUJU21, but not with the nNNPDF3.0 nuclear PDFs.

###  5.6 Exclusive and inclusive observables in ultraperipheral collisions

Ultraperipheral collisions (UPCs) of ions are interactions in which the
approaching nuclei do not touch. Instead, they interact at a distance due to
their strong electromagnetic fields [284, 285]. In comparison to typical
minimum-bias $p$Pb (let alone PbPb) collisions, much fewer background
processes take place, and the signal processes are thus easier to isolate.

The exclusive production of $J/\psi$ mesons in UPCs has triggered particular
interest. The process is dominated by the exchange of an almost real photon.
In photon-nucleus collisions, the PDFs appear already at the level of the
matrix element

$\mathcal{M}(\gamma+A\rightarrow J/\psi+A)\sim T_{g}\otimes
f_{g}^{A}+\sum_{q}T_{q}\otimes f_{q}^{A}\,.$ (51)

When squared to obtain a cross section, the latter becomes extremely sensitive
to PDFs. Several LO studies have been performed in the past [286, 287, 288,
289], but the first NLO calculations for PbPb collisions have appeared only
very recently [290, 291] despite the fact that the NLO coefficient functions
$T_{g,q}$ have been known for some time [292]. Figure 8 (right) compares NLO
calculations with several recent nuclear PDFs with the combined experimental
data from the LHC [277, 278, 279, 280]. The factorization scales have been
chosen to match the ALICE data at mid-rapidity. While the central theory
values do not reproduce the behavior of the data particularly well, the
nuclear PDF error bands are much wider than the data uncertainties. This
indicates that these data should further constrain the nuclear PDFs.
Unfortunately, the process is perturbatively unstable: at LO only gluons
contribute, while at NLO there is also a contribution from the quark singlet,
which can even dominate at NLO. The reason is that the LO and NLO gluon
contributions enter $T_{g}$ with opposite signs and there is a significant
cancellation between the two. It has been argued that the theoretical
uncertainties could be brought under better control by summing logarithmically
enhanced contributions at small $x$ [293] and by considering power corrections
in the coefficient functions [294]. There are also other theoretical
uncertainties associated e.g. with modeling of the photon flux and how nuclear
generalized PDFs (GPDs) and collinear PDFs are related [295].

In addition to exclusive observables, also inclusive processes such as dijet
photoproduction are sensitive to nuclear PDFs. Here a photon emitted from one
nucleus breaks up the other to produce a hadronic final state that contains
two hard jets. This process has been calculated in LO with parton showers
[296] and in NLO [297]. Preliminary measurements by ATLAS also exist [298,
299], and the 2017 measurement qualitatively agrees with the NLO calculation.
However, the ATLAS data make use of forward neutrons and rapidity gaps to
resolve the photon-going direction on an event-by-event basis, and this
excludes the diffractive component of UPCs, which is included in NLO
calculations with standard PDFs [300]. The imposed experimental conditions
also require further modeling associated with preventing a Coulomb break-up of
the photon-emitting nucleus and with finite-size effects.

##  6 Other developments

###  6.1 Electron Ion Collider

The Electron Ion Collider (EIC) currently under construction at BNL will
extend the kinematic region of lepton-nucleus scattering compared to fixed
target experiments by one order of magnitude in $x$ (to a few times 10-4) and
$Q^{2}$ (to $10^{3}$ GeV2). This region is of course still considerably
smaller than the one accessible at the LHC (cf. Fig. 3), but the environment
will be much cleaner than in $pA$ collisions, making it easier to disentangle
cold nuclear matter (i.e. leading twist, factorizable) effects from other
possible (higher twist, non-factorizable) contributions. In addition, the EIC
will allow for measurements on a large variety of nuclei in NC and CC DIS as
well as in photoproduction and with inclusive, semi-inclusive (e.g. identified
hadron) and exclusive (e.g. diffractive) final states [301, 302].

While $F_{2}^{A}$ is sensitive to the momentum distributions of gluons mainly
through scaling violations (cf. Eq. 37), the EIC will also allow for
measurements of the longitudinal structure function [303, 304]

$F_{L}^{A}(x,Q^{2})=\frac{\alpha_{s}(Q^{2})}{2\pi}x^{2}\int_{x}^{1}\frac{dz}{z^{3}}\left[\frac{8}{3}F_{2}(x,Q^{2})+4\sum_{q}e_{q}^{2}\left(1-\frac{x}{z}\right)zf_{g}^{A}(x,Q^{2})\right],$
(52)

which has a direct contribution from gluons. Measuring $F_{L}^{A}$ will,
however, require operation at different center-of-mass energies. The impact of
inclusive NC DIS measurements at the EIC on global fits of nuclear PDFs has
been investigated in Refs. [302, 305, 306] in the frameworks of nCTEQ15WZ,
EPPS16 and nNNPDF2.0 with several nuclei and beam-energy configurations (see
also Ref. [103]). At low $Q^{2}$, the predicted impact is significant for all
partonic flavors other than the strange quark. At higher $Q^{2}$, the better
constrained gluon also leads to a better constrained strange quark PDF. Charm
tagging allows to access the reduced charm cross section

$\sigma^{c\bar{c}}_{\rm
red}=\frac{d^{2}\sigma^{c\bar{c}}}{dxdQ^{2}}\frac{xQ^{4}}{2\pi\alpha^{2}[1+(1-y)^{2}]}=F_{2}^{c\bar{c}}-\frac{y^{2}}{1+(1-y)^{2}}F_{L}^{c\bar{c}}$
(53)

related to the charm structure functions $F_{2,L}^{c\bar{c}}$, sensitive to
the gluon and a potential intrinsic charm content in the nucleon, and thus to
further reduce the uncertainties [305, 307]. Dijets in DIS [308, 309] and
photoproduction [310, 311, 312] as well as charm jets [313] provide further
information, the latter in particular on the strange quark in CC DIS.
Similarly to the LHC, exclusive vector meson production is highly sensitive to
the gluon [314, 315]. Diffractive final states with a large rapidity gap or
identified hadrons in the forward direction will allow to access the
completely unknown territory of diffractive nuclear PDFs [316, 317].

Taken together, the EIC measurements will allow for a greater parametric
flexibility in the $x$-dependence at the starting scale $Q_{0}^{2}$ similar to
proton PDFs, lead to a more reliable $A$-dependence, make parameterizations
possible not just in $A$, but also in $Z$, and allow to move away from the
nuclear stability line and to study mirror nuclei. The EIC should also help to
answer the question of potential different nuclear effects in CC and NC DIS,
shed light on shadowing, gluon saturation [318, 319, 320], transverse momentum
distributions (TMDs) [321, 322], the transition to the color glass condensate
(CGC) [323, 324] as well as the EMC effect across a wide range of $A$ and
energy scales. With polarized beams of light nuclei (2H, 3H, 3He), even the
polarized EMC effect could be investigated [301, 302].

###  6.2 Lattice QCD

In the non-perturbative lattice QCD approach, the four-dimensional space-time
is discretized and QCD regularized on a finite Euclidean lattice. Correlation
functions are then computed numerically in the path integral formalism using
methods adapted from statistical mechanics, and the results are extrapolated
to the continuum and infinite volume limits. To make contact with experimental
data, lattice QCD calculations must demonstrate control over all sources of
systematic uncertainty including discretization effects, extrapolation from
unphysical pion masses, finite-volume effects, and renormalization of
composite operators.

Light-cone quantities like PDFs cannot be calculated directly on a Euclidean
lattice. Instead, the traditional approach has been to determine the matrix
elements of local twist-two operators that can be related to the Mellin
moments of PDFs. In principle, given a sufficient number of Mellin moments,
PDFs can be reconstructed from the inverse Mellin transform. In practice,
however, the calculation is limited to the lowest three moments, since power-
divergent mixing occurs between twist-two operators. Three moments are
insufficient to fully reconstruct the momentum dependence of the PDFs without
significant model dependence. The lowest three moments do provide, however,
useful information both as benchmarks of lattice calculations and as
constraints in global extractions of PDFs [325].

Direct extractions of the $x$ dependence of PDFs have been attempted based on
quasi-PDFs in large-momentum effective theory [326], pseudo-PDFs [327] and
other methods. Quasi-PDFs are defined as Fourier transforms of the matrix
elements, whereas pseudo-PDFs are transforms in Ioffe time. For positive
($u-d$) and negative ($\bar{u}-\bar{d}$) isovector quark combinations, the
quasi-PDFs at the physical pion mass were found to agree with global fits at
large $x>0.1$ and $x>0.4$, respectively [328]. Gluon quantities are much
noisier than quark disconnected loops and require calculations with very high
statistics. Up to perturbative matching and power corrections, the Fourier
transform of the gluon quasi-PDF was found to be compatible with the one of
global fits within the statistical uncertainty [329]. Calculating the
small-$x$ behavior requires larger boost momenta, as this results in a faster
decay of the matrix elements, so that truncations in the Fourier transform
matter less [330]. Fourier transforms of strange and charm quark PDFs have
also been obtained, the former being about five times larger and both being
smaller than those in global fits, possibly due to missing contributions from
other flavor distributions. A full analysis of lattice QCD systematics must
still be performed. Nevertheless, the strangeness asymmetry ($s-\bar{s}$) in
the region of $0.3<x<0.8$ was found to be very small with high precision
compared to the uncertainty in global fits [331]. Including these lattice data
in a global fit therefore greatly reduces the size of the $s-\bar{s}$ error
band in the large-$x$ region [332].

Lattice QCD studies of nuclear structure are currently restricted to low $A$
and unphysical pion masses. In particular, the ratio of the longitudinal
momentum fraction carried by the positive isovector quark combination in 3He
to the one in the free nucleon was found to be consistent with unity at the
few-percent level. This is in agreement with, but more precise than, current
determinations from global fits. Including this lattice result in the
nNNPDF2.0 global fit framework reduces the uncertainty on the isovector
momentum fraction ratio by a factor of 2.5 and leads to a more precise
extraction of the $u$ and $d$ quark distributions in 3He [333]. Previously,
the nuclear modification of the gluon momentum fraction for $A\leq 3$ was
found to be less than $\sim$10%. This is consistent with expectations from
phenomenological quark distributions and the momentum sum rule [334].

###  6.3 Relations to other phenomena

The global analysis of nuclear PDFs outlined in this review indicates that, at
sufficiently high interaction scales, collinear factorization is a consistent
way to describe lepton-nucleus and proton-nucleus collisions. However, towards
low interaction scales there are theoretical reasons to believe that the role
of higher-twist effects, suppressed by inverse powers of the interaction
scale, become increasingly important in comparison to simpler lepton-proton
and proton-proton collisions. In the parton model, these effects arise from
processes where two initial-state partons recombine [335, 336, 337, 338]. In
general, these processes modify the linear DGLAP evolution by slowing down the
evolution at small $x$ and increasing it at intermediate $x$. Since the
spatial density of partons is higher in large nuclei, the effect should be
more pronounced there and lead to dynamically generated shadowing and
antishadowing [339]. At very low $Q^{2}$, the recombination of an even larger
number of partons becomes eventually important and, when resummed, can be
interpreted as saturation [4, 340]. Finding a conclusive signature of
saturation has proven rather difficult, calculations based on nuclear PDFs and
the saturation picture giving very similar results. Even forward $D$-meson
production in $p$Pb collisions, which probes the nucleus down to $x\sim
10^{-5}$ at nearly non-perturbative interaction scales, has not revealed clear
deviations from the linear DGLAP dynamics. In addition, the resummation of
small-$x$ BFKL logarithms in the language of PDFs also have a tendency to slow
down the linear DGLAP evolution [341, 342], which further complicates the
search for true saturation effects.

There is a heated discussion revolving around the question whether in very
central, high-multiplicity $p$Pb collisions – and even in smaller systems like
$pp$ or $\gamma A$ collisions – one creates a state of matter that has
collective liquid-type properties – a droplet of Quark-Gluon plasma (QGP)
[343]. Some characteristic features have been experimentally observed [344],
which typically consist of correlations between particles. However, the
measurements can also be interpreted in terms of the initial-state geometry
[345]. Furthermore, the way in which final-state QCD particles hadronize into
a color-neutral state has been shown to display features that could
accidentally be attributed to a liquid-like behavior [346]. It should also be
kept in mind that jet quenching has not been observed in $p$Pb collisions,
which challenges the QGP picture of small systems and is a concrete difference
between observations in $p$Pb and PbPb collisions. Multiplicity-integrated
$p$Pb cross sections appear to be consistent with collinear factorization and
process-independent nuclear PDFs. The relations between the observations
discussed above remain open questions at this moment.

In the case of heavy-ion collisions, the formation of a QGP is nowadays a
generally accepted phenomenon. Nevertheless, even in heavy-ion collisions the
LHC data for electroweak boson [169, 173, 347, 348] and high-$p_{T}$ direct
photon production [349, 350, 351] are consistent with collinear factorization
and process-independent nuclear PDFs [352]. Thus, there is no reason to
believe that the initial state of heavy-ion collisions would not be dictated
by nuclear PDFs. This idea has been pursued in the Eskola-Kolhinen-Ruuskanen-
Tuominen model of heavy-ion collisions, whose most recent versions [353] apply
NLO perturbative QCD calculations and impact-parameter dependent nuclear PDFs
[354] to compute the inital conditions for the subsequent fluid-dynamical
evolution of the system.

The nuclear PDFs also find use in the field of neutrino astronomy and cosmic-
ray physics [355, 356, 357]. Interactions of neutrinos coming from outer space
can be measured in large neutrino telescopes such as IceCube, KM3NeT and
Baikal, where the neutrinos interact with water or ice. Precise theoretical
calculations of the cross sections require nuclear PDFs as an input. In
addition, protons from astrophysical sources can collide with the air
molecules in the atmosphere, which can also produce neutrinos. Precise
calculations of the cross sections for these secondary neutrinos require
nuclear PDFs as well.

##  7 Conclusion

During the last 25 years of research in nuclear PDFs, the field has undergone
an enormous development. Methodologically, simple fits, performed by eye at
LO, have matured into rigorous statistical analyses, including machine-
learning techniques, at NLO and NNLO with full error estimates. Nonetheless,
global nuclear-PDF analyses are still driven by experimental measurements, and
in this respect the $p$Pb collisions carried out during the past decade at the
LHC have opened up a wide, previously unexplored regime both in terms of
kinematics and processes. The following items summarize our take-home
messages:

[SUMMARY POINTS]

  1. 1.

Despite the theoretical and experimental advances, there are still significant
differences among the independent global analyses of nuclear PDFs, both in
terms of the extracted nuclear modifications of PDFs as well as the absolute
nuclear PDFs. In several places the central values of a given analysis can be
outside the error bands of the others, and in some cases even the error bands
of two given analyses do not overlap. The widths of the error bands can be
also very different. The most significant factors behind the observed
differences can be attributed to (i) the assumed form of the non-perturbative
parameterization of nuclear PDFs at low $Q^{2}$, (ii) the data selection,
(iii) fitting absolute cross sections vs. ratios of cross sections, (iv) the
theoretical treatment of heavy flavors, and (v) different baseline free-proton
PDFs. As highlighted in Figs. 4–8, these differences also lead to visible
effects in observables. This underscores the need to carry out the global
analysis independently in several groups. To faithfully chart the theoretical
uncertainties in quantities that depend on nuclear PDFs, it is thus
recommended to use more than one set of nuclear PDFs.

  2. 2.

At the moment, NLO accuracy is the standard in the field of nuclear PDFs, the
full NNLO accuracy being limited by the existence or public availability of
cross-section codes. Most $pA$ data can be well described within the NLO
calculations, but in the case of some observables there is a confirmed (DY
below the $Z$ peak) or conjectured (prompt photons, jets) need for NNLO QCD.
The NNLO accuracy also reduces the theoretical uncertainties in observables
sensitive to effects of partonic saturation, offers a standard candle to
Glauber modeling of heavy-ion collisions through precise predictions for
electroweak observables, and should lead to more precise predictions for
astrophysical applications.

  3. 3.

In the long run, the global analysis of nuclear PDFs should be extended to
include the proton and the deuteron. As of today, most of the free-proton fits
still utilize heavy-target – particulary neutrino DIS – data to constrain the
full flavor decomposition. These same data can be (and are), however, also
taken as constraints on the nuclear modifications of PDFs. To fully chart the
interplay between the two requires a simultaneous extraction of the free-
proton and nuclear PDFs.

## DISCLOSURE STATEMENT

The authors are not aware of any affiliations, memberships, funding, or
financial holdings that might be perceived as affecting the objectivity of
this review.

## ACKNOWLEDGMENTS

The authors thank their nCTEQ and EPPS colleagues for their collaboration and
useful discussions, E. Nocera for providing the nNNPDF3.0 values for Fig. 7,
and V. Guzey for providing the nCTEQ15HQ values for the right-hand panel of
Fig. 8. MK thanks his ALICE colleagues for their collaboration and
acknowledges funding by the BMBF through project 05P21PMCAA and DFG through
GRK 2149 and SFB 1225 “Isoquant,” project-id 273811115. HP acknowledges
funding by the Academy of Finland through the Center of Excellence in Quark
Matter, project 346326. The results shown in Fig. 6 and in the left-hand panel
of Fig. 8 have been computed for this review using computing resources of the
Finnish IT Center for Science (CSC), project jyy2580.

## References

  * [1]
  * [2] Frankfurt LL, Strikman MI. Phys. Rept. 160:235 (1988)
  * [3] Arneodo M. Phys. Rept. 240:301 (1994)
  * [4] Gelis F, Iancu E, Jalilian-Marian J, Venugopalan R. Ann. Rev. Nucl. Part. Sci. 60:463 (2010)
  * [5] Klasen M, Klein-Bösing C, König F, Wessels JP. JHEP 10:119 (2013)
  * [6] Andronic A, Braun-Munzinger P, Redlich K, Stachel J. Nature 561:321 (2018)
  * [7] Gribov VN, Lipatov LN. Sov. J. Nucl. Phys. 15:438 (1972)
  * [8] Gribov VN, Lipatov LN. Sov. J. Nucl. Phys. 15:675 (1972)
  * [9] Dokshitzer YL. Sov. Phys. JETP 46:641 (1977)
  * [10] Altarelli G, Parisi G. Nucl. Phys. B 126:298 (1977)
  * [11] Arnold RG et al. Phys. Rev. Lett. 52:727 (1984)
  * [12] Gomez J et al. Phys. Rev. D 49:4348 (1994)
  * [13] Amaudruz P et al. [New Muon Collaboration]. Nucl. Phys. B 441:3 (1995)
  * [14] Tzanov M et al. [NuTeV Collaboration]. Phys. Rev. D 74:012008 (2006)
  * [15] Arneodo M et al. [New Muon Collaboration]. Nucl. Phys. B 481:23 (1996)
  * [16] Caldwell DO et al. Phys. Rev. Lett. 42:553 (1979)
  * [17] Goodman MS et al. Phys. Rev. Lett. 47:293 (1981)
  * [18] Bauer TH, Spital RD, Yennie DR, Pipkin FM. Rev. Mod. Phys. 50:261 (1978)
  * [19] Armesto N. J. Phys. G 32:R367 (2006)
  * [20] Frankfurt L, Guzey V, Strikman M. Phys. Rept. 512:255 (2012)
  * [21] Kopeliovich BZ, Morfin JG, Schmidt I. Prog. Part. Nucl. Phys. 68:314 (2013)
  * [22] Atwood WB, West GB. Phys. Rev. D 7:773 (1973)
  * [23] Bodek A, Ritchie JL. Phys. Rev. D 23:1070 (1981)
  * [24] Frankfurt LL, Strikman MI. Nucl. Phys. B 181:22 (1981)
  * [25] Saito K, Uchiyama T. Z. Phys. A 322:299 (1985)
  * [26] Ciofi degli Atti C, Simula Y. Phys. Rev. C 53:1689 (1996)
  * [27] Segarra EP et al. Phys. Rev. D 103:114015 (2021)
  * [28] Aubert JJ et al. [European Muon Collaboration]. Phys. Lett. B 123:275 (1983)
  * [29] Bodek A et al. Phys. Rev. Lett. 50:1431 (1983)
  * [30] Bodek A et al. Phys. Rev. Lett. 51:534 (1983)
  * [31] Geesaman DF, Saito K, Thomas AW. Ann. Rev. Nucl. Part. Sci. 45:337 (1995)
  * [32] Thomas AW. Int. J. Mod. Phys. E 27:1840001 (2019)
  * [33] Alde DM et al. [E772 Collaboration]. Phys. Rev. Lett. 64:2479 (1990)
  * [34] Malace S, Gaskell D, Higinbotham DW, Cloet I. Int. J. Mod. Phys. E 23:1430013 (2014)
  * [35] Hen O, Miller GA, Piasetzky E, Weinstein LB. Rev. Mod. Phys. 89:045002 (2017)
  * [36] Nikolaev NN, Zakharov VI. Phys. Lett. B 55:397 (1975)
  * [37] Eskola KJ. Nucl. Phys. B 400:240 (1993)
  * [38] Qiu JW. Nucl. Phys. B 291:746 (1987)
  * [39] Eskola KJ, Kolhinen VJ, Salgado CA. Eur. Phys. J. C 9:61 (1999)
  * [40] Hirai M, Kumano S, Miyama M. Phys. Rev. D 64:034003 (2001)
  * [41] de Florian D, Sassot R. Phys. Rev. D 69:074028 (2004)
  * [42] Hirai M, Kumano S, Nagai TH. Phys. Rev. C 76:065207 (2007)
  * [43] Eskola KJ, Paukkunen H, Salgado CA. JHEP 04:065 (2009)
  * [44] Kovarik K et al. Phys. Rev. D 93:085037 (2016)
  * [45] Schienbein I et al. Phys. Rev. D 77:054013 (2008)
  * [46] Paukkunen H, Salgado CA. JHEP 07:032 (2010)
  * [47] Kovarik K et al. Phys. Rev. Lett. 106:122301 (2011)
  * [48] Paukkunen H, Salgado CA. Phys. Rev. Lett. 110:212301 (2013)
  * [49] Nakamura SX et al. Rept. Prog. Phys. 80:056301 (2017)
  * [50] Kalantarians N, Keppel C, Christy ME. Phys. Rev. C 96:032201 (2017)
  * [51] Duwentäster P et al. Phys. Rev. D 105:114043 (2022)
  * [52] Eskola KJ, Paakkinen P, Paukkunen H, Salgado CA. Eur. Phys. J. C 82:413 (2022)
  * [53] Abdul Khalek R et al. Eur. Phys. J. C 82:507 (2022)
  * [54] Gao J, Harland-Lang L, Rojo J. Phys. Rept. 742:1 (2018)
  * [55] Kovařík K, Nadolsky PM, Soper DE. Rev. Mod. Phys. 92:045003 (2020)
  * [56] Ethier JJ, Nocera ER. Ann. Rev. Nucl. Part. Sci. 70:43 (2020)
  * [57] Klein SR, Mäntysaari H. Nature Rev. Phys. 1:662 (2019)
  * [58] Georgi H, Politzer HD. Phys. Rev. D 14:1829 (1976)
  * [59] Collins JC, Soper DE, Sterman GF. Adv. Ser. Direct. High Energy Phys. 5:1 (1989)
  * [60] Nadolsky PM, Tung WK. Phys. Rev. D 79:113014 (2009)
  * [61] Thorne RS, Tung WK. arXiv:0809.0714 [hep-ph] (2008)
  * [62] Aivazis MAG, Collins JC, Olness FI, Tung WK. Phys. Rev. D 50:3102 (1994)
  * [63] Krämer M, Olness FI, Soper DE. Phys. Rev. D 62:096007 (2000)
  * [64] Tung WK, Kretzer S, Schmidt C. J. Phys. G 28:983 (2002)
  * [65] Cacciari M, Greco M, Nason P. JHEP 05:007 (1998)
  * [66] S. Forte, E. Laenen, P. Nason and J. Rojo, Nucl. Phys. B 834 (2010), 116-162
  * [67] Accardi A, Qiu JW. JHEP 07:090 (2008)
  * [68] Ruiz R et al. arXiv:2301.07715 [hep-ph] (2023)
  * [69] Nachtmann O. Nucl. Phys. B 63:237 (1973)
  * [70] Barnett RM. Phys. Rev. Lett. 36:1163 (1976)
  * [71] Accardi A et al. arXiv:hep-ph/0308248 (2003)
  * [72] d’Agostini G. Nucl. Instrum. Meth. A 346:306 (1994)
  * [73] Ball RD et al. [NNPDF Collaboration]. JHEP 05:075 (2010)
  * [74] Pumplin J, Stump DR, Tung WK. Phys. Rev. D 65:014011 (2001)
  * [75] Pumplin J et al. Phys. Rev. D 65:014013 (2001)
  * [76] Giele WT, Keller S. Phys. Rev. D 58:094023 (1998)
  * [77] Giele WT, Keller SA, Kosower DA. arXiv:hep-ph/0104052 [hep-ph] (2001)
  * [78] Nadolsky PM, Sullivan Z. eConf C010630:P510 (2001)
  * [79] Helenius I, Walt M, Vogelsang W. Phys. Rev. D 105:9 (2022)
  * [80] Khanpour H et al. Phys. Rev. D 104:034010 (2021)
  * [81] Forte S, Garrido L, Latorre JI, Piccione A. JHEP 05:062 (2002)
  * [82] Ball RD et al. [NNPDF Collaboration]. Nucl. Phys. B 809:1-63 (2009) [erratum: Nucl. Phys. B 816:293 (2009)]
  * [83] Pumplin J et al. JHEP 07:012 (2002)
  * [84] Stump D et al. JHEP 10:046 (2003)
  * [85] Owens JF et al. Phys. Rev. D 75:054030 (2007)
  * [86] Kusina A et al. Eur. Phys. J. C 80:968 (2020)
  * [87] Duwentäster P et al. Phys. Rev. D 104:094005 (2021)
  * [88] Eskola KJ, Helenius I, Paakkinen P, Paukkunen H. JHEP 05:037 (2020)
  * [89] Kusina A, Lansberg JP, Schienbein I, Shao HS. Phys. Rev. D 104:014010 (2021)
  * [90] Kom CH, Kulesza A, Stirling WJ. Phys. Rev. Lett. 107:082002 (2011)
  * [91] Kniehl BA, Kramer G, Schienbein I, Spiesberger H. Phys. Rev. D 71:014018 (2005)
  * [92] Butenschoen M, Kniehl BA. Phys. Rev. Lett. 106:022003 (2011)
  * [93] Eskola KJ, Kolhinen VJ, Ruuskanen PV. Nucl. Phys. B 535:351 (1998)
  * [94] Eskola KJ, Kolhinen VJ, Paukkunen H, Salgado CA. JHEP 05:002 (2007)
  * [95] Eskola KJ, Paukkunen H, Salgado CA. JHEP 07:102 (2008)
  * [96] Eskola KJ, Paakkinen P, Paukkunen H, Salgado CA. Eur. Phys. J. C 77:163 (2017)
  * [97] Hou TJ et al. Phys. Rev. D 103:014013 (2021)
  * [98] Accardi A, Hobbs TJ, Jing X, Nadolsky PM. Eur. Phys. J. C 81:603 (2021)
  * [99] Ball RD et al. [NNPDF Collaboration]. Eur. Phys. J. C 77:663 (2017)
  * [100] Ball RD et al. [NNPDF Collaboration]. Eur. Phys. J. C 82:428 (2022)
  * [101] Rosenblatt M. The Annals of Mathematical Statistics 27:832 (1956)
  * [102] Parzen E. The Annals of Mathematical Statistics 33:1065 (1962)
  * [103] Abdul Khalek R et al. [NNPDF Collaboration]. Eur. Phys. J. C 79:471 (2019)
  * [104] M. Walt, I. Helenius and W. Vogelsang, Phys. Rev. D 100:096015 (2019)
  * [105] Bertone V et al. [xFitter Developers’ Team]. PoS DIS2017:203 (2018)
  * [106] Abramowicz H et al. [H1 and ZEUS Collaborations]. Eur. Phys. J. C 75:580 (2015)
  * [107] Sick I, Day D. Phys. Lett. B 274:16 (1992)
  * [108] Denniston AW et al. MS-TP-22-13
  * [109] Kulagin SA, Petti R. Nucl. Phys. A 765:126 (2006)
  * [110] Kulagin SA, Petti R. Phys. Rev. D 76:094023 (2007)
  * [111] Kulagin SA, Petti R. Phys. Rev. C 90:045204 (2014)
  * [112] Ru P, Kulagin SA, Petti R, Zhang BW. Phys. Rev. D 94:113013 (2016)
  * [113] Muzakka KF et al. Phys. Rev. D 106:074004 (2022)
  * [114] Ashman J et al. [European Muon Collaboration]. Phys. Lett. B 202:603 (1988)
  * [115] de Florian D, Sassot R, Zurita P, Stratmann M. Phys. Rev. D 85:074028 (2012)
  * [116] Arneodo M et al. [New Muon Collaboration]. Nucl. Phys. B 481:3 (1996)
  * [117] Adams MR et al. [E665 Collaboration]. Z. Phys. C 67:403 (1995)
  * [118] Prytz K. Phys. Lett. B 311:286 (1993)
  * [119] Eskola KJ, Honkanen H, Kolhinen VJ, Salgado CA. Phys. Lett. B 532:222 (2002)
  * [120] Airapetian A et al. [HERMES Collaboration]. arXiv:hep-ex/0210068 [hep-ex] (2002)
  * [121] Armesto N, Paukkunen H, Salgado CA, Tywoniuk K. Phys. Lett. B 694:38 (2011)
  * [122] Ackerstaff K et al. [HERMES Collaboration]. Phys. Lett. B 475:386 (2000) [erratum: Phys. Lett. B 567:339 (2003)]
  * [123] Ball RD et al. [NNPDF Collaboration]. Nature 608:483 (2022)
  * [124] Aubert JJ et al. [European Muon Collaboration]. Phys. Lett. B 94:96 (1980)
  * [125] Aubert JJ et al. [European Muon Collaboration]. Nucl. Phys. B 213:31 (1983)
  * [126] Ashman J et al. [European Muon Collaboration]. Z. Phys. C 57:211 (1993)
  * [127] Seely J et al. [Hall C Collaboration]. Phys. Rev. Lett. 103:202301 (2009)
  * [128] Arrington J et al. [Hall C Collaboration]. Phys. Rev. C 104:065203 (2021)
  * [129] Karki A et al. [Hall C Collaboration]. arXiv:2207.03850 [nucl-ex] (2022)
  * [130] Schmookler B et al. [CLAS Collaboration]. Nature 566:354 (2019)
  * [131] Paukkunen H, Zurita P. Eur. Phys. J. C 80:381 (2020)
  * [132] Accardi A et al. Phys. Rev. D 93:114017 (2016)
  * [133] Arrington J et al. Phys. Rev. C 73:035205 (2006)
  * [134] Martin AD et al. Eur. Phys. J. C 73:2318 (2013)
  * [135] Moreno G et al. [E605 Collaboration]. Phys. Rev. D 43:2815 (1991)
  * [136] Vasilev MA et al. [E866 Collaboration]. Phys. Rev. Lett. 83:2304 (1999)
  * [137] Ellis SD, Stirling WJ. Phys. Lett. B 256:258 (1991)
  * [138] Martin AD, Roberts RG, Stirling WJ, Thorne RS. Eur. Phys. J. C 4:463 (1998)
  * [139] Lin PJ. doi:10.2172/1398791
  * [140] Helenius I et al. Phys. Rev. D 100:014004 (2019)
  * [141] Aaij R et al. [LHCb Collaboration]. LHCb-CONF-2018-005
  * [142] Paakkinen P, Eskola KJ, Paukkunen H. Phys. Lett. B 768:7 (2017)
  * [143] Badier J et al. [NA3 Collaboration]. Phys. Lett. B 104:335 (1981)
  * [144] Bordalo P et al. [NA10 Collaboration]. Phys. Lett. B 193:368 (1987)
  * [145] Heinrich JG et al. [E615 Collaboration]. Phys. Rev. Lett. 63:356 (1989)
  * [146] Adams B et al. [COMPASS++/AMBER Collaboration]. arXiv:1808.00848 [hep-ex] (2018)
  * [147] Ball RD et al. [NNPDF Collaboration]. Eur. Phys. J. C 79 (2019) no.3, 282
  * [148] Mason D et al. [NuTeV Collaboration]. Phys. Rev. Lett. 99:192001 (2007)
  * [149] Zeller GP et al. [NuTeV Collaboration]. Phys. Rev. Lett. 88:091802 (2002) [erratum: Phys. Rev. Lett. 90:239902 (2003)]
  * [150] Hirai M, Kumano S, Nagai TH. Phys. Rev. D 71:113007 (2005)
  * [151] Eskola KJ, Paukkunen H. JHEP 06:008 (2006)
  * [152] Workman RL et al. [Particle Data Group]. PTEP 2022:083C01 (2022)
  * [153] Berge JP et al. Z. Phys. C 49:187 (1991)
  * [154] Yang UK et al. [CCFR Collaboration]. Phys. Rev. Lett. 86:2742 (2001)
  * [155] Onengut G et al. [CHORUS Collaboration]. Phys. Lett. B 632:65 (2006)
  * [156] Goncharov M et al. [CCFR/NuTeV Collaboration]. Phys. Rev. D 64:112006 (2001)
  * [157] Samoylov O et al. [NOMAD Collaboration]. Nucl. Phys. B 876:339 (2013)
  * [158] Abreu H et al. [FASER Collaboration]. Phys. Rev. Lett. 131:3 (2023)
  * [159] Albanese R et al. [SND@LHC Collaboration]. Phys. Rev. Lett. 131:031802 (2023)
  * [160] Cruz-Martinez JM et al. arXiv:2309.09581 [hep-ph] (2023)
  * [161] Carli T et al. Eur. Phys. J. C 66:503 (2010)
  * [162] Aad G et al. [ATLAS Collaboration]. Phys. Rev. C 92:044915 (2015)
  * [163] Khachatryan V et al. [CMS Collaboration]. Phys. Lett. B 759:36 (2016)
  * [164] Adam J et al. [ALICE Collaboration]. JHEP 02:077 (2017)
  * [165] Aaij R et al. [LHCb Collaboration]. JHEP 09:030 (2014)
  * [166] Aad G et al. [ATLAS Collaboration]. ATLAS-CONF-2015-056
  * [167] Khachatryan V et al. [CMS Collaboration]. Phys. Lett. B 750:565 (2015)
  * [168] Sirunyan AM et al. [CMS Collaboration]. JHEP 05:182 (2021)
  * [169] Acharya S et al. [ALICE Collaboration]. JHEP 09:076 (2020)
  * [170] Aaij R et al. [LHCb Collaboration]. JHEP 06:022 (2023)
  * [171] Khachatryan V et al. [CMS Collaboration]. Eur. Phys. J. C 76:469 (2016)
  * [172] Sirunyan AM et al. [CMS Collaboration]. Phys. Lett. B 800:135048 (2020)
  * [173] Acharya S et al. [ALICE Collaboration]. JHEP 05:036 (2023)
  * [174] Acharya S et al. [ALICE Collaboration]. Phys. Rev. C 102:055204 (2020)
  * [175] Berger EL, Gordon LE, Klasen M. Phys. Rev. D 58:074012 (1998)
  * [176] Klasen M, Brandt M. Phys. Rev. D 88:054002 (2013)
  * [177] Brandt M, Klasen M, König F. Nucl. Phys. A 927:78 (2014)
  * [178] Arleo F, Gousset T. Phys. Lett. B 660:181 (2008)
  * [179] Brenner Mariotto C, Goncalves VP. Phys. Rev. C 78:037901 (2008)
  * [180] Arleo F, Eskola KJ, Paukkunen H, Salgado CA. JHEP 04:055 (2011)
  * [181] Helenius I, Eskola KJ, Paukkunen H. JHEP 09:138 (2014)
  * [182] Klasen M. Rev. Mod. Phys. 74:1221 (2002)
  * [183] Glück M, Reya E, Vogt A. Phys. Rev. D 48:116 (1993) [erratum: Phys. Rev. D 51:1427 (1995)]
  * [184] Bourhis L, Fontannaz M, Guillet JP. Eur. Phys. J. C 2:529 (1998)
  * [185] Klasen M, König F. Eur. Phys. J. C 74:3009 (2014)
  * [186] Frixione S. Phys. Lett. B 429:369 (1998)
  * [187] Apanasevich L et al. [Fermilab E706 Collaboration]. Phys. Rev. D 70:092009 (2004)
  * [188] Adare A et al. [PHENIX Collaboration]. Phys. Rev. C 87:054907 (2013)
  * [189] Abelev BI et al. [STAR Collaboration]. Phys. Rev. C 81:064904 (2010)
  * [190] Aaboud et M al. [ATLAS Collaboration]. Phys. Lett. B 796:230 (2019)
  * [191] Aad G et al. [ATLAS Collaboration]. JHEP 08:005 (2016)
  * [192] Campbell JM, Ellis RK, Williams C. Phys. Rev. Lett. 118 (2017) no.22, 222001 [erratum: Phys. Rev. Lett. 124 (2020) no.25, 259901]
  * [193] Acharya S et al. [ALICE Collaboration]. CERN-LHCC-2020-009
  * [194] Acharya S et al. [ALICE Collaboration]. ALICE-PUBLIC-2023-001
  * [195] Vogt R. Phys. Rev. C 70:064902 (2004)
  * [196] Arsene I et al. [BRAHMS Collaboration]. Phys. Rev. Lett. 91:072305 (2003)
  * [197] Arsene I et al. [BRAHMS Collaboration]. Phys. Rev. Lett. 93:242303 (2004)
  * [198] Adler SS et al. [PHENIX Collaboration]. Phys. Rev. Lett. 91:072303 (2003)
  * [199] Adams J et al. [STAR Collaboration]. Phys. Rev. Lett. 91:072304 (2003)
  * [200] Adler SS et al. [PHENIX Collaboration]. Phys. Rev. Lett. 98:172302 (2007)
  * [201] Adams J et al. [STAR Collaboration]. Phys. Lett. B 637:161 (2006)
  * [202] Kniehl BA, Kramer G, Pötter B. Nucl. Phys. B 582:514 (2000)
  * [203] Aidala CA et al. Phys. Rev. D 83:034002 (2011)
  * [204] de Florian D et al. Phys. Rev. D 91:014035 (2015)
  * [205] de Florian D et al. Phys. Rev. D 95:094019 (2017)
  * [206] Moffat E et al. [Jefferson Lab Angular Momentum Collaboration]. Phys. Rev. D 104:016015 (2021)
  * [207] Adare A et al. [PHENIX Collaboration]. Phys. Rev. C 88:024906 (2013)
  * [208] Acharya S et al. [ALICE Collaboration]. Eur. Phys. J. C 78:624 (2018)
  * [209] Adam J et al. [ALICE Collaboration]. Phys. Lett. B 760:720 (2016)
  * [210] Acharya S et al. [ALICE Collaboration]. Phys. Lett. B 827:136943 (2022)
  * [211] Sassot R, Stratmann M, Zurita P. Phys. Rev. D 81:054001 (2010)
  * [212] Zurita P. arXiv:2101.01088 [hep-ph] (2021)
  * [213] Aaij R et al. [LHCb Collaboration]. Phys. Rev. Lett. 131:042302 (2023)
  * [214] Aaij R et al. [LHCb Collaboration]. arXiv:2310.17326 [nucl-ex] (2023)
  * [215] Aaij R et al. [LHCb Collaboration]. Phys. Rev. Lett. 128:142004 (2022)
  * [216] Aad G et al. [ATLAS Collaboration]. JHEP 07:074 (2023)
  * [217] Khachatryan V et al. [CMS Collaboration]. Eur. Phys. J. C 75:237 (2015)
  * [218] Khachatryan V et al. [CMS Collaboration]. JHEP 04:039 (2017)
  * [219] Acharya S et al. [ALICE Collaboration]. JHEP 11:013 (2018)
  * [220] d’Enterria D, Eskola KJ, Helenius I, Paukkunen H. Nucl. Phys. B 883:615 (2014)
  * [221] Loizides C, Kamin J, d’Enterria D. Phys. Rev. C 97:054910 (2018) [erratum: Phys. Rev. C 99:019901 (2019)]
  * [222] Sirunyan AM et al. [CMS Collaboration]. Phys. Rev. Lett. 121:062002 (2018)
  * [223] Chatrchyan S et al. [CMS Collaboration]. Eur. Phys. J. C 74:2951 (2014)
  * [224] Eskola KJ, Paakkinen P, Paukkunen H. Eur. Phys. J. C 79:511 (2019)
  * [225] Currie J et al. Phys. Rev. Lett. 119:152001 (2017)
  * [226] Aad G et al. [ATLAS Collaboration]. Phys. Lett. B 748:392 (2015)
  * [227] Khachatryan V et al. [CMS Collaboration]. Eur. Phys. J. C 76:372 (2016)
  * [228] Acharya S et al. [ALICE Collaboration]. arXiv:2307.10860 [nucl-ex] (2023)
  * [229] Cacciari M, Mangano ML, Nason P. Eur. Phys. J. C 75:610 (2015)
  * [230] Gauld R, Rojo J, Rottoli L, Talbert J. JHEP 11:009 (2015)
  * [231] Gauld R, Rojo J. Phys. Rev. Lett. 118:072001 (2017)
  * [232] Zenaiev O et al. [PROSA Collaboration]. Eur. Phys. J. C 75:396 (2015)
  * [233] Gauld R. Phys. Rev. D 93:014001 (2016)
  * [234] Kusina A, Lansberg JP, Schienbein I, Shao HS. Phys. Rev. Lett. 121:052004 (2018)
  * [235] Helenius I, Paukkunen H. JHEP 05:196 (2018)
  * [236] Helenius I, Paukkunen H. JHEP 07:054 (2023)
  * [237] Nason P. JHEP 11:040 (2004)
  * [238] Frixione S, Nason P, Oleari C. JHEP 11:070 (2007)
  * [239] Alioli S, Nason P, Oleari C, Re E. JHEP 06:043 (2010)
  * [240] Lansberg JP, Shao HS. Eur. Phys. J. C 77:1 (2017)
  * [241] Aad G et al. [ATLAS Collaboration]. Phys. Rev. C 92:034904 (2015)
  * [242] Aaboud M et al. [ATLAS Collaboration]. Eur. Phys. J. C 78:171 (2018)
  * [243] Sirunyan AM et al. [CMS Collaboration]. Eur. Phys. J. C 77:269 (2017)
  * [244] Tumasyan A et al. [CMS Collaboration]. Phys. Lett. B 835:137397 (2022)
  * [245] Sirunyan AM et al. [CMS Collaboration]. Phys. Lett. B 790:509 (2019)
  * [246] Sirunyan AM et al. [CMS Collaboration]. Phys. Lett. B 772:306 (2017)
  * [247] Khachatryan V et al. [CMS Collaboration]. Phys. Lett. B 754:59 (2016)
  * [248] Abelev BB et al. [ALICE Collaboration]. Phys. Rev. Lett. 113:232301 (2014)
  * [249] Adam J et al. [ALICE Collaboration]. Phys. Rev. C 94:054908 (2016)
  * [250] Acharya S et al. [ALICE Collaboration]. JHEP 12:092 (2019)
  * [251] Abelev BB et al. [ALICE Collaboration]. JHEP 02:073 (2014)
  * [252] Adam J et al. [ALICE Collaboration]. JHEP 06:055 (2015)
  * [253] Acharya S et al. [ALICE Collaboration]. JHEP 06:011 (2022)
  * [254] Abelev BB et al. [ALICE Collaboration]. Phys. Lett. B 740:105 (2015)
  * [255] Abelev BB et al. [ALICE Collaboration]. JHEP 12:073 (2014)
  * [256] Acharya S et al. [ALICE Collaboration]. JHEP 01:178 (2022)
  * [257] Aaij R et al. [LHCb Collaboration]. JHEP 10:090 (2017)
  * [258] Aaij R et al. [LHCb Collaboration]. JHEP 02:072 (2014)
  * [259] Aaij R et al. [LHCb Collaboration]. JHEP 07:094 (2014)
  * [260] Acharya S et al. [ALICE Collaboration]. JHEP 07:160 (2018)
  * [261] Acharya S et al. [ALICE Collaboration]. JHEP 07:137 (2023)
  * [262] Acharya S et al. [ALICE Collaboration]. Phys. Lett. B 806:135486 (2020)
  * [263] Acharya S et al. [ALICE Collaboration]. JHEP 07:237 (2020)
  * [264] Bezshyiko I et al. [LHCb Collaboration]. Phys. Rev. Lett. 131:102301 (2023)
  * [265] Aaij R et al. [LHCb Collaboration]. Phys. Lett. B 774:159 (2017)
  * [266] Aaij R et al. [LHCb Collaboration]. JHEP 11:194 (2018) [erratum: JHEP 02:093 (2020)]
  * [267] Aaij R et al. [LHCb Collaboration]. Phys. Rev. D 99:052011 (2019)
  * [268] Aaij R et al. [LHCb Collaboration]. Phys. Rev. Lett. 122:132002 (2019)