Parton Distributions in Nucleons and Nuclei

We review the current status of Parton Distribution Function (PDF) determinations for unpolarized and longitudinally polarized protons and for unpolarized nuclei, which are probed by high-energy hadronic scattering in perturbative Quantum Chromodynamics (QCD). We present the established theoretical framework, the experimental information, and the methodological aspects inherent to any modern PDF extraction. Furthermore, we summarize the present knowledge of PDFs and discuss their limitations in both accuracy and precision relevant to advance our understanding of QCD proton substructure and pursue our quest for precision in the Standard Model and beyond. In this respect, we highlight various achievements, discuss contemporary issues in PDF analyses, and outline future directions of progress.


A BRIDGE BETWEEN LOW AND HIGH ENERGIES
Nucleons (protons and neutrons) are bound states that make up all nuclei, and hence most of the visible matter in the Universe.Unraveling their fundamental structure in terms of their elementary constituents -quarks and gluons, or collectively called partonsis currently one of the main challenges at the boundary of hadron and particle physics.Such an understanding is rooted in the theoretical framework provided by the Standard Model (SM), part of which contains the field theory that describes the strong interactions of color-charged partons: Quantum Chromodynamics (QCD).Since energy grows with the separation between color charges, one of the defining features of QCD, partons are confined to exist only in neutral color combinations called hadrons, among which are nucleons.
Nucleons are probed by scattering a beam of leptons or protons/antiprotons from them in large-momentum-transfer processes.Since the elementary interactions occur at distances much shorter than the confinement scale (the scale at which partons are bound into nucleons), the measurable cross section of the process can be determined by folding the partonic cross section with Parton Distribution Functions (PDFs).The former encodes the scattering of quasi-free partons in terms of process-dependent kernels that are computed perturbatively in QCD, while the latter detail the momentum distribution of partons that enter the elementary scattering process in terms of universal functions.
Parton distributions are ubiquitous in hadron and particle physics: they are essential tools to interpret experimental data for a variety of hard-scattering processes in light of the underlying theory.Such processes are measured with the greatest precision by different experiments at a number of facilities around the world, among which is the largest proton collider ever built: the Large Hadron Collider (LHC).Collectively, they cover a wide range of momentum-transfer energies and are therefore sensitive to different aspects of the theory that are not yet fully understood.Among these are the following outstanding questions: Are the Higgs boson dynamics the same as those prescribed by the SM, or are there other phenomena at work?What is the real origin of electroweak symmetry breaking and particles' mass?How do nucleon bound states emerge from parton interactions?How does the proton spin emerge from quarks and gluons?How are parton dynamics modified within the bound nucleons of nuclei?What is the interplay with astrophysics?
To successfully address -and possibly answer -these questions, a careful determination of PDFs and their uncertainties is mandatory.Currently, this cannot be reliably achieved from first principles.Instead, PDFs are modeled by means of some parameterization, which is then optimized by comparing the PDF-dependent prediction of one or more physical process to its actual measurement, a procedure that is called (global) QCD analysis (or fit).In this sense, PDF extractions can be labeled generally as a nonlinear regression problem, whereby one has to learn a set of functions from data.The PDFs can then be used to describe any other process that depends on them due to their universal property.
The purpose of this review is twofold.First of all, it aims at providing an accessible, yet complete, introduction to the art of PDF fitting with a focus on the theoretical, experimental, and methodological input that enters modern QCD analyses.In this respect, our review is more concise and pedagogical than those recently presented in (1)(2)(3).Secondly, it aims at presenting altogether various species of collinear PDFs: unpolarized PDFs, relevant for precision physics in the SM; helicity (or polarized) PDFs, relevant for the nucleon spin structure; and nuclear PDFs, relevant for the understanding of cold nuclear matter effects.In this respect, our review updates and extends those in (4,5), and presents a critical snapshot of PDFs as a bridge between low and high energy physics.
The review is organized as follows.In Sect. 2 we present the theoretical formalism used to study PDFs in perturbative QCD.In Sect. 3 we discuss the typical data sets included in a global QCD analysis, and how they constrain the different parton flavors.We review the methodological aspects relevant to PDF fitting and techniques used to represent their uncertainties in Sect. 4. Sect. 5 is an overview of state-of-the-art PDF fits, their main features, and phenomenological implications.This is followed by a discussion on the relevance of PDFs to some of the questions outlined above, and by presenting a critical appraisal on future prospects, in Sect.6.Finally, we conclude with an outlook in Sect.7.

THEORETICAL INPUT
In this section we review the standard theoretical framework used to study PDFs in perturbative QCD.In particular, we describe the factorization of hadronic observables, the details of the perturbative computations, and the theoretical constraints that PDFs must fulfill.

Factorization of hadronic observables
Factorization theorems (6) provide the framework in which perturbative QCD calculations can be performed for a class of sufficiently inclusive hadronic observables that are measured in large-momentum-transfer processes.In the region of asymptotic freedom (7), elementary QCD interactions occur at distances much shorter than the confinement scale -roughly the order of the nucleon's quantum-mechanical wavelength -and short-distance interactions of partons (expressed in terms of process-dependent partonic cross sections calculable in perturbative QCD) can be seperated from their long-distance momentum distributions (given by the non-perturbative, universal PDFs).Factorization realizes the convolution of these two parts to make predictions of experimental observables.While we refer to (8,9) for a detailed treatment of the subject, in the following we briefly discuss how factorization behaves for lepto-and hadro-production, two classes of processes that are key to this review.
Deep-inelastic scattering (DIS), the high-energy inclusive scattering of leptons and nucleons N ( N → X, where X is unobserved), still plays a central role in discussing factorization, both because it was the first process to establish the existence of partons inside the nucleon, and because it comprises a majority of the experimental information used in PDF analyses (see Sect. 3).The unpolarized (polarized) DIS cross sections σ (∆σ) for neutral current (i=NC) interactions involving a photon or Z 0 -boson exchange, and for charged current (i=CC) interactions given by a W ± -boson exchange, can be written as ( 10) where the standard DIS kinematic variables are: the momentum fraction x = Q 2 /2q • P , the inelasticity y = q • P/k • P , and the energy Q 2 = −q 2 , with k, P , and q being the lepton, nucleon, and transferred momentum, respectively.The polarized cross section is defined as the difference ∆σ i = σ i (λN = −1, λ ) − σ i (λN = +1, λ ), where λN (λ ) is the nucleon (lepton) helicity which can have value ±1 corresponding to its orientation parallel (+) or antiparallel (−) to the beam direction.The electroweak factors η i depend on the interaction type and on the lepton helicity λ , where GF is the Fermi constant and MW the W -boson mass.In Eqs. 1. and 2., we have denoted the Quantum Electrodynamics (QED) coupling as α, defined the kinematic factors Y± = 1±(1−y) 2 , and neglected terms proportional to M 2 /Q 2 , with M being the mass of the nucleon.In Eqs. 1. and 3. the ∓ (±) sign refers to either an incoming electron or neutrino (−), or to an incoming positron or antineutrino (+).At leading twist, the (un)polarized structure functions F i 2 , F i 3 , F i L and g i 4 , g i 1 , g i L , collectively labeled F i , factorize as, where the sum runs over all active quarks n f and the gluon, C i f are the partonic scattering coefficients computed as a power series in the strong coupling as = αs/4π, and f N are either the unpolarized (f N ) or the polarized (∆f N ) PDFs in a nucleon N .Finally, the PDFs can be further decomposed in terms of the net momentum density of partons aligned (↑) or anti-aligned (↓) to the parent nucleon's polarization: The PDFs also have a bilocal operator definition in canonical field theory (see Ref. (6) and Sect. 2 in Ref.
(2)).For example, the unpolarized quark PDFs can be written as, where ψ f , ψf are the quark fields, γ are the Dirac matrices, W is the Wilson line ensuring the Gauge invariance of the operator, and all four-vectors have been expressed using light-cone coordinates.Similar expressions also hold for the gluon PDF and for polarized PDFs.Additional particles in the final state can be measured in conjunction with the outgoing lepton in lepton-nucleon scattering.In semi-inclusive DIS (SIDIS), for instance, a single hadron type h is identified ( N → hX).In the case of unpolarized incoming particles, the SIDIS cross section reads (11), The process is sufficiently inclusive to allow the structure functions F h 1,L to factorize as where the hard scattering coefficients are now denoted as C h f,f , and the fragmentation function (FF) D h f (12) is introduced.This is the time-like counterpart of the PDF encoding the hadronization of a parton f into a hadron h.The corresponding momentum fraction is defined as z = P • P h /P • q, where P h is the four-momentum of the outgoing hadron.
The cross section for a generic unpolarized nucleon-nucleon process, which depends on a single scale Q (e.g.pp → AX, where A can be, for example, a single jet of hadrons, a lepton or quark pair, or an electroweak boson) also has a factorized expression, where s is the center-of-mass energy of the hadronic collision, τ = Q 2 /s is the scaling variable of the hadronic process, σ 0 f f is the leading order (LO) partonic cross section, f N and f ,N are the unpolarized PDFs stemming from each nucleon, and C f f are the hard scattering coefficients.Moreover, expressions similar to Eq. 9. exist for processes with one or both colliding nucleons polarized, for factorizable multi-scale processes (e.g.Higgs production in W fusion), for less inclusive processes (e.g.electroweak boson production in association with jets), and for one-particle inclusive production (e.g.neutral pion production).In the last case, an additional convolution with the final particle FF must be included.The factorized result of Eq. 9. can also be extended to total and differential cross sections.
Should hard scattering processes occur off nuclei instead of nucleons, it is custom to assume that they can still be described in terms of factorization theorems (13).If so, the partonic scattering coefficients in Eqs.4., 8., and 9. remain the same but the PDFs are of nuclei, defined as the average of proton and neutron densities bound in a nucleus, In this expression, the pair of atomic numbers (A, Z) identifies the nuclear isotope, and f p,n/A are the proton and neutron bound PDFs, usually related to the nucleon PDFs by, Since the nuclear and nucleon PDFs are defined by the same leading twist operators of Eq. 6.
(though acting on different states), it is natural to assume that nuclear modifications can be absorbed into the PDFs without altering factorization theorems.This assumption might not hold for processes or kinematic regions subject to sizable higher-twist corrections (2).

Perturbative calculations
The partonic scattering coefficients C i f , C h f f and C f f entering Eqs.4., 8. and 9., which we collectively denote as C, can be expressed as a perturbative series in the strong coupling, where the explicit form of the coefficients c (k) depend on the specific process.At the lowest perturbative order, LO, they either vanish or are proportional to a Dirac delta, in which case the hadronic cross sections reduce to a combination of the PDFs that couple to the relevant final state.Each process is therefore sensitive to different partons (see Sect. 3).The computation of partonic cross sections to higher orders, usually next-to leading and next-to-next-to-leading orders (NLO and NNLO), generally entail three classes of singularities: ultraviolet (UV) singularities, which are renormalized through the running of the QCD coupling as(µR), at a scale µR; infrared singularities associated with loop graphs, which cancel corresponding soft singularities from the emission of real gluons; and collinear singularities due to the initial partons emitting gluons at zero angle, which are subtracted by terms arising in the renormalization of the PDF operators (6).The resulting partonic cross sections are therefore free of any infrared type singularities, but remain dependent on a renormalization scale µR as a result of regularizing the UV divergences.Renormalization also implies the choice of a scheme, the one most widely used being the modified minimal subtraction scheme (MS) (14), which includes additional common factors in the counterterms required to make finite predictions.
Since the PDFs defined in factorization are renormalized separately, an additional scale dependence µF arises in their argument, the choice of which can be made independently from the scale choice in the partonic cross section calculation.However, both the renormalization and the factorization scales µR and µF are unphysical, i.e. the hadronic cross section does not depend on them should it be computed to all orders in perturbation theory.At a given fixed order a k s these cancellations will only be approximate, that is, for the hadronic cross section σ of any of the processes in Eqs.1.,2.,7.and 9..This suggests that the dependence of the hadronic cross section on the renormalization and factorization scale choice decreases as one carries out the perturbative calculation to higher orders (unless new kinematic configurations open up), and that scale variations can be used to estimate the accuracy of the perturbative order.Unless otherwise stated, we assume µR = µF = Q in this review (including the previous section).
The dependence of the strong coupling on the renormalization scale and of the PDFs on the factorization scale is completely determined by the invariance of the renormalization group and the corresponding equations: the QCD beta function for the strong coupling (which is currently known up to five loops (15)), and the DGLAP evolution equations (16) for the PDFs.This is a set of 2n f + 1 coupled integro-differential equations of the form, often conveniently rewritten in terms of the gluon, non-singlet f N NS = q − q and singlet f N S = n f q q + (with q + = q + q) quark-antiquark combinations.In this way, Eq. 14. reduces to 2n f − 1 decoupled equations for the non-singlet distributions and two coupled equations for the gluon and singlet distributions.The splitting functions P f f can be either unpolarized (P f f ) or polarized (∆P f f ) and can be expanded as a power series in as, with the coefficients p (k) computed in perturbative QCD up to NNLO (17) and partly to N 3 LO (18) for P f f , and NNLO for ∆P f f (19).Expressions similar to Eqs. 14.-15.also hold for the FF in Eq. 8., with the corresponding splitting functions known up to NNLO (20).
Since the coefficients for partonic cross sections in Eq. 12. are known up to NNLO for most of the unpolarized processes and up to NLO for most of the polarized processes, the corresponding PDFs can be determined only up to these accuracies.Beyond LO, the coefficients of partonic scattering in Eq. 12. and of the splitting functions in Eq. 15. contain terms proportional to ln Q 2 , ln(1/x) and ln(1 − x).While the first appears in the splitting kernels and are summed by the DGLAP equations, resummation techniques (21) and BFKL equations ( 22) may be used to deal with ln(1/x) terms at small x.At large x, threshold resummation techniques (23) handle terms proportional to ln(1−x).
Sums over the number of flavors n f appear in all factorization formulas and in the evolution equations.Decoupling arguments (24) imply that the contribution of heavy quark flavors to any process are power-suppressed at scales which are below the threshold for their production (25).Therefore, when expressing predictions for processes at different scales in terms of the same PDF set, it is necessary to use a variable-flavor number (VFN) scheme in which different numbers of active flavors are adopted consistently.In the vicinity of the threshold for heavy quark production, the quark mass m cannot be ignored.This can be accounted for in a general-mass VFN (GM-VFN) scheme that interpolates, in a modeldependent way, between the fixed-flavor number (FFN) scheme near production threshold and the asymptotic result of the zero-mass VFN (ZM-VFN).In the FFN scheme, heavy quark mass effects are built into the partonic cross sections, but terms proportional to large logarithms of the form ln(Q/m) are not resummed; in the ZM-VFN scheme, heavy quark mass effects are ignored, but ln(Q/m) terms are instead resummed into the heavy quark PDFs.Massive coefficient functions are known up to O(α 2 s ) for both neutral-( 26) and charged-current (27) DIS.A number of partial results also exist at O(α 3 s ), see (28) and references therein.Various GM-VFN schemes have been worked out in the literature up to NNLO, including the exact dependence on the heavy-quark mass up to O(α 2 s ): ACOT (29), TR (30) and FONLL (31).These schemes differ by subleading terms, which may not be entirely negligible at NLO in the vicinity of the quark threshold, but rapidly decrease at NNLO (32).In the case of FONLL, the GM-VFNS was generalized to include mass effects for heavy-quark initiated contributions (33) should non-perturbative charm and bottom quark content in the proton be required.While most of the details have been worked out in the unpolarized case, these can be extended to polarized PDFs and even FFs (34).
We shall finally note that any process involving electroweak interactions also receives higher-order electromagnetic or weak corrections.Both coefficent functions in Eq. 12. and splitting functions in Eq. 15. should be modified, and additional PDFs for the photon (and generally for leptons) should be included in factorization formulas and evolution equations (35).Since the QED coupling at the electroweak scale is α ∼ α 2 s ∼ 0.01, one expects NLO corrections in the electromagnetic interaction to be roughly of the same order of NNLO QCD corrections.Electroweak corrections have never been systematically included in global analyses of PDFs so far, except for unpolarized photon-initiated partonic cross sections.In this case, it was demonstrated (36) that the photon PDF can be almost completely determined (to a precision comparable to quark PDFs) by relating the entire photon contribution to the proton structure function (refer to Sect.7 in (1) for a thorough review).

Theoretical constraints
There are several theoretical constraints that can be used to help determine PDFs.
Positivity of hadronic observables.Hadronic observables must always be positive regardless of the shape or size of the PDFs, which can be negative for unpolarized PDFs beyond the lowest perturbative order, and for polarized PDFs at all orders.Positivity can be enforced by choosing ad hoc PDF parameterizations or imposed by defining a set of control cross sections on a grid of kinematics that captures a sufficiently large region of the phase space.The parameter configurations that lead to negative cross sections are subsequently discarded.In principle, one could also isolate analytically the terms that lead to negative results from collinear subtraction in the DGLAP equations.Such a procedure was used to derive positivity bounds on polarized PDFs from their unpolarized counterparts (37).
Sum rules.Energy conservation implies that unpolarized parton momentum fractions x carried by each parton must sum to unity, i.e.PDFs must fulfill the momentum sum rule, 16.
In addition, valence sum rules ensure that an unpolarized hadron maintains its valence structure, i.e. for a proton made of two up quarks and one down quark, SU f symmetry relations.There are constraints that can be imposed on polarized PDFs using an established relationship between weak baryon decays and the non-singlet combinations of PDF moments.The first follows from the Bjorken sum rule (38), which relates the difference between the lowest moments of the proton and neutron structure functions g p 1 and g n 1 to the isovector axial charge gA assuming exact SU f (2) symmetry.The second follows from relating the non-singlet combination of PDFs to the octet axial charge a8, taken from weak hyperon decay assuming exact SU f (3) symmetry.In terms of PDF moments, The second relationship in particular plays a vital role in constraining the strange polarization, forcing its first moment to have a value ∼ −0.1.While SU f (2) symmetry has been confirmed in a recent QCD analysis to almost 2% (39), there is some uncertainty surrounding the level of SU f (3) symmetry breaking in the value of a8 (40).As a result, imposing the hyperon decay value in a fit can bias the shape of the strange polarization (see Sect. 5).
Nuclear boundary condition.In addition to momentum and valence sum rules, nuclear PDFs for low A nuclei are constrained by unpolarized proton PDFs at A = 1.This boundary condition is usually implemented by defining a nuclear PDF parameterization that exactly reproduces the central value of the proton PDF at A = 1.Alternatively (41), one can fit the A = 1 distribution along with other nuclei, and then discard parameter configurations that deviate from a given proton boundary condition (with its uncertainties and correlations).This approach results in smaller uncertainties for neighboring nuclei.

EXPERIMENTAL INPUT
Theoretical predictions for hadronic cross sections, obtained by the formalism presented in Sect.2, must be compared to experimental data in order to determine PDFs.Table 1 presents a summary of the main hadronic processes used in current global QCD analyses.
For each reaction, we indicate the leading partonic process contributing to the corresponding factorization formulas, the probed parton flavors, and whether there exists available data in the unpolarized, polarized, and nuclear sectors.Furthermore, hadronic processes are grouped into seven categories: fixed target DIS, collider DIS, fixed target SIDIS, fixed target Drell Yan (DY), collider DY, jet and hadron production, and top production.The kinematic coverage in x and Q 2 of the corresponding measurements is displayed in Fig. 1 for each PDF species.Applicability of perturbative QCD implies that only the data above Q 2 = 1 GeV 2 are usually considered (indicated by the dashed horizontal line in Fig. 1).Since each process probes a different parton flavor (or combinations thereof) in different kinematic regions, a global QCD analysis is required to piece the information together.The ability to describe such a variety of hadronic observables simultaneously is in fact a great success of perturbative QCD.Restricted datasets are also sometimes considered in order to maximize consistency across measurements.However, this may limit the precision and accuracy of the resulting PDF distributions.
As is apparent in Fig. 1, the amount of experimental data points and the variety of hadronic measurements is significantly larger for unpolarized proton scattering than for polarized and nuclear reactions.The reason for this is that, historically, accelerating polarized nucleons and heavy ions has been technically more challenging.As a result, the experimental facilities pivotal for the extraction of unpolarized proton PDFs (namely the HERA collider and the LHC) have largely remained unparalleled in the polarized and nuclear cases.However, several new facilities and upgrades are expected in the near future that aim to decrease this discrepancy (see Sect. 6).
Nevertheless, in all cases, the bulk of the experimental information is provided by DIS.In the unpolarized sector, neutral and charged current DIS has been measured by various experiments at the CERN Super-Synchrotron (CERN-SPS), such as NMC, BCDMS and CHORUS, at the Stanford Linear Accelerator (SLAC), at Fermilab (NuTeV) and at HERA.In the polarized case, neutral current DIS has been measured by various experiments at the CERN-SPS (EMC, SMC, COMPASS), at SLAC, and at JLab.Lastly, in the nuclear case, both neutral and charged current DIS has been measured by experiments at the CERN-SPS, SLAC and Fermilab.In principle, DIS is able to constrain all partons: quark-antiquark total and valence combinations for each flavor through neutral and charged current interactions, respectively, and the gluon through scaling violations via DGLAP evolution and higher Jet and hadron production pp(p{p, A}) → (di−)jet(s) + X gg, qg, qq → jet(s) g, q pp(pp) → h + X gg, qg, qq → π, K, D g, q Table 1 The main hadronic processes commonly used to determine PDFs.For each reaction, the leading partonic process, the probed partons, and whether available data exists in the unpolarized (U), polarized (P), and nuclear (N) cases is indicated.
order QCD corrections in the structure functions.Complimentary to DIS are DY and W ± /Z 0 production measurements which are necessary to disentangle the sea quark distributions.Various measurements of both fixed-target and collider DY scattering have been performed: in the unpolarized case by the E665 and E866 experiments at Fermilab, by the D0 and CDF experiments at the Tevatron, and by the ATLAS, CMS and LHCb experiments at the LHC; in the polarized case, by the STAR and PHENIX experiments at RHIC; and in the nuclear case, by the E772, E866, E605 experiments at Fermilab.Sea quark distributions can likewise be probed in SIDIS, although analyses of the corresponding data are complicated by the non-perturbative finalstate FF that enters the theoretical prediction (Eq.8.).Semi-inclusive pion, kaon, and charged hadron production were measured in both the unpolarized and polarized cases by The kinematic coverage, in the (x, Q 2 ) plane, of the hadronic cross section data for the processes commonly included in global QCD analyses of unpolarized, polarized, and nuclear PDFs.The extended kinematic ranges attained by the LHeC and the EIC are also displayed.
the COMPASS experiment at CERN and by the HERMES experiment at HERA.
For the gluon distribution, jet and hadron production hold much of the constraining power.In the unpolarized sector, jet production has been measured both in DIS by the H1 and ZEUS experiments at HERA, in proton-(anti)proton collisions by the D0 and CDF experiments at the Tevatron, and by the ATLAS, CMS and ALICE experiments at the LHC.In the polarized case, jet and pion production has been measured in proton-proton collisions by the STAR and PHENIX experiments at RHIC.In the nuclear case, jet production and hadron production have been measured in proton-ion collisions by the CMS and ALICE experiments at the LHC and by the STAR and PHENIX experiments at RHIC.Finally, single top and top pair production is used to probe the sea quark and gluon distributions, respectively.In this case, measurements are only available for unpolarized PDF analyses from the ATLAS and CMS experiments at the LHC.
The interested reader can obtain further details on all of these measurements (including their references) from Sect. 3 of Ref. (1) for unpolarized PDFs, from Sect. 3 of Ref. (3) for the polarized distributions, and from Sect. 3 of Ref. (42) for the nuclear PDFs.The impact of these observables on PDFs will be discussed in Sect. 5.

METHODOLOGICAL INPUT
Since parton distributions cannot be computed reliably from first principles, they are instead determined by comparing theoretical predictions of hadronic cross sections in the form of Eqs.4.,8.and 9. to experimental data.In other words, PDFs must be modeled by a function of fit parameters and optimized by data to yield matching theoretical predictions.It can be classified more generally as a nonlinear regression problem, in which a confidence interval in the space of PDFs is determined by optimizing a suitable goodness-of-fit measure (43).In this section, we discuss the methodological aspects to perform a global QCD anaylsis, in particular the PDF parameterization, optimization, and uncertainty representation.

Parameterization
Unpolarized and polarized proton PDFs are usually parameterized for each parton f as, at an input scale Q 2 0 .The power-like factors x α f and (1−x) β f describe the low-x and high-x behavior of PDFs, and are inspired by general QCD arguments, namely Regge theory (44) and Brodsky-Farrar quark counting rules (45) for the small-and large-x asymptotic limits, respectively.While various models can predict approximate values for α f and β f (46), they are determined from the data in global QCD analyses.The factor N is an overall normalization that accounts for theoretical constraints given by Eqs.16.-18., which also imply β f > 0 (to ensure PDF integrability).The function I(x; a), which depends on a set of parameters a, is defined to interpolate between the small-and large-x regions.Typically it is chosen to have a simple polynomial form, e.g.
, where a = {γ, δ, . . .} are parameters to be determined by the fit, although some have used more flexible Chebyshev (47,48) or Bernstein polynomials (49).The function I can also be defined as the output of a neural network (50,51), to reduce bias coming from the choice of parameterization as much as possible.In this context, neural networks provide a convenient unbiased set of (nonlinear) basis functions.
For PDFs of nuclei, an additional dependence on the atomic mass A is required.The treatment of such dependence is subject to various choices.For instance, the nuclear modification R A f (x, Q0) can be parameterized directly (42), in which case the PDFs of a proton bound in a nucleus can be constructed according to Eq. 11., given a proton PDF set.Alternatively, the nuclear proton PDF can be parameterized similarly to the free proton PDFs (52) by giving an A-dependent functional form to the parameters in Eq. 19, (α f , β f , a) → (α f (A), β f (A), a(A)).Lastly, one can remain agnostic about the A dependence by modeling it as part of a neural network (41).
The choice of the input scale Q0 is typically set below the charm mass threshold mc if heavier quark PDFs are assumed to be generated by QCD radiation, or slightly above it if the charm PDF is parameterized on the same footing as the light quark PDFs.Should the FFN scheme be chosen (53), various input scales are used above each threshold leading to different PDF sets for every number of active flavors.For unpolarized PDFs, there are usually seven independent PDFs to parameterize, f = u, ū, d, d, s, s, g, which increase to eight or nine if the charm (54) and the photon (55) PDFs are also parameterized.Should a restricted data set be used, or different species of PDFs be fitted, the number of independent PDFs must be adjusted to match the partons to which the data is sensitive.Any linear combination of independent PDFs can be used, for example the sum and difference of the quark and anti-quark distributions, q + ≡ q+q and q − ≡ q−q, or the singlet and non-singlet distributions of the DGLAP evolution basis 14..

Optimization
Once the parameterization is set at an input scale Q0, the PDFs can be evolved to the scale of the data Q by means of the DGLAP equations 14, and then convoluted with partonic cross sections to obtain theoretical predictions for the hadronic cross sections at a given perturbative order.The optimal PDF parameters are then found by optimizing a suitable figure of merit, commonly taken as the log-likelihood χ 2 .Given a set of N dat measurements Di, and a set of corresponding theoretical predictions Ti({a}), the χ 2 is defined as where the covariance matrix covij is composed of all of the uncertainties of the data, namely, uncorrelated (σ unc i ) and ncorr correlated (σ corr,A k,i ) uncertainties.The χ 2 also reads which is equivalent to Eq. 20. upon minimization with respect to the shift parameters, rj.This allows one to study also the behavior of the shift parameters at the minimum, where their distribution ought to be a univariate Gaussian with mean zero.
If correlated multiplicative uncertainties are provided (e.g. in the case of normalization errors such as the luminosity uncertainty), they must be treated carefully in order to avoid the so-called d'Agostini bias (56), which results in the optimal parameters giving a theoretical prediction that underestimates the data.To avoid this, multiplicative uncertainties in Eq.21 can be treated iteratively, in which they are multiplied by central theory predictions from a previous fit instead of by the experimental data.This procedure, called the t0-method (57), was proven to rapidly converge.Alternatively, one can fit an overall normalization parameter and allow it to fluctuate within the multiplicative uncertainty range.This can be naturally included in Eq. 22 as one of the rj parameters.
The optimization of Eq. 20. is performed numerically, which means DGLAP evolution and convolutions between PDFs and partonic cross sections must be performed quickly, accurately, and precisely in order to make global fits viable.In addition, fitting strategies must be optimized to explore a sufficiently large region of PDF parameter space.
Several publicly available codes solve DGLAP evolution equations efficiently (58-60).These programs, as well as most private evolution codes used in PDF fits, have been benchmarked against standard tables (60,61).The calculation of hadronic cross sections is performed via the use of lookup tables, where partonic cross sections convoluted with evolution kernels are precomputed and stored for each data set on a suitable interpolation grid.Hadronic cross sections then reduce to the scalar product between such interpolation tables and the PDFs parameterized at the scale Q0.This method is usually realized as part of each private optimization code (see e.g. ( 62)), except for the APFELgrid program (63).In this case, DGLAP evolution kernels and DIS partonic cross sections are provided by the APFEL package (59), while partonic cross sections for other unpolarized processes are provided in the format of APPLgrid (64) and FastNLO (65) tables.The latter are obtained in turn from multi-purpose Monte Carlo generators (66), which are accurate to a given perturbative order and contain appropriate numerical interfaces (67).
The choice of an efficient optimization strategy largely depends on the dimension of the PDF parameter space.Numerical gradient-based algorithms such as the Newton's method, the Levenberg-Marquardt method (which supplements the Newton's method with the steepest descent method), or variable metric methods (which only rely on gradient information) are used for spaces of moderate dimension (roughly on the order of 40 free parameters).If the dimension of the PDF parameter space is larger, these methods start to become unsuitable due to inefficient numerical inversions of large matrices and the increased possibility of ending in a local minima.Stochastic genetic algorithms (68) or deterministic gradient descent methods, possibly associated with deep-learning and hyperoptimization techniques (69), can be used to efficiently explore the parameter space in this case.To avoid fitting the noise in the data, it is important to devise a suitable stopping criterion.In this respect, a widely used method in the literature is cross-validation, where the data points are randomly divided into two sets: training and validation.The χ 2 is then computed for both sets separately, but optimized only on the training set.The fit terminates when the χ 2 of the validation set starts to increase (while the χ 2 of the training set continues to decrease).To avoid information loss, the procedure should be repeated a sufficiently large number of times starting from different data partitions.

Representation of PDF uncertainties
The confidence interval in the space of PDFs, namely a representation of the PDF probability density, is intrinsically derived from the probability density of the fitted parameters given the data.In Bayesian language we are interested in P(a|D), in which the expectation E and variance V of an observable O depending on a set of PDFs can be computed as, In this section, we review the statistical estimation of E and V with focus on the three components that contribute to V : experimental, procedural, and theoretical uncertainties.
Experimental uncertainties.There exists at least two commonly used approaches for propagating the data uncertainties into the PDFs: the Hessian and Monte Carlo methods.The Hessian technique (70) (including its Lagrange multiplier generalization (71)), adopted when a fairly limited parameterization is used, is based on the standard leastsquares approach, and assumes that P(a|D) is multi-Gaussian.Once the best-fit is determined, the χ 2 is approximated to first nontrivial order about its minimum.The desired confidence level is obtained as the volume in parameter space about the minimum that corresponds to a fixed increase in χ 2 .For Gaussian uncertainties, the 68% (or one-sigma) confidence level corresponds to the volume enclosed by ∆χ 2 = χ 2 − χ 2 min = 1.Given the set of parameters a0 that minimize the χ 2 and the Hessian matrix H in the parameter space, whose elements are Hij = 1 2 ∂ 2 χ 2 (a) ∂a i ∂a j a=a 0 , it follows that Here, npar is the number of free parameters and neig is the number of eigenvectors of the matrix H.The confidence level is therefore entirely determined by the inverse of the Hessian matrix, which is equivalent to the covariance matrix in parameter space.Uncertainties on the PDFs, and on any quantities depending on them, can be compactly represented by providing eigenvectors of the Hessian matrix rescaled by their eigenvalues.They can then be computed by adding in quadrature the variation along the direction of each eigenvector by a fixed amount, ∆O k = O(ai) − (a0) (more sophisticated formulas hold for asymmetric eigenvectors).However, data sets can be incompatible, in which case the prescription ∆χ 2 = 1 for the parameter shifts may lead to unrealistically small uncertainties.To compensate this, the scaling along the eigendirections is inflated by a tolerance factor T = ∆χ 2 > 1.
Its specific value is usually chosen by studying the distribution of best-fit parameters across experiments or by determining a different tolerance along each Hessian eigenvector (72).
In addition to the tolerance criteria, there are two additional disadvantages to the Hessian method.The first is that it becomes rapidly unviable when the parameter space is large.Secondly, the Gaussian assumption, and also a linear approximation for the observable O that is made in deriving Eq. 24, might become inadequate to handle flat directions that arise whenever the data loosely constrain the PDF parameters.
A more robust method is the Monte Carlo technique, which samples the probability distribution P(a|D) by bootstrapping the starting data sample (i.e.smearing the experimental data about their respective uncertainties) and performing a fit to each of these new data replicas.The distribution of the data is therefore mapped into the distribution of the PDFs, provided that the number of data replicas nrep is sufficiently large (typically on the order of a few hundred).The expectation value and uncertainty of any quantity depending on the PDFs can therefore be computed as a mean and variance over the PDF ensemble, where a 0,k are the best-fit parameters for each replica.Despite being computationally expensive, the Monte Carlo method has some obvious advantages with respect to the Hessian method: any distribution can be used to smear the data, should they be non-Gaussian; there is no need to rely on the linear approximation of observables; there are no limitations in how large the parameter space can be; and the probability density P(a|D) can be readily updated with Bayes' theorem to incorporate a new piece of experimental information without performing a new fit.The last of these is realized through reweighting (73), which consists of assigning each PDF replica a weight that quantifies its agreement with the new data.In this case, formulas in Eq. 25. should be replaced by their weighted counterparts.Since replicas with small weights become immaterial, the PDF ensemble looses part of its statistical power, a drawback of the method when such a loss is excessive.In any case, Monte Carlo uncertainties are statistically sound in that they can be rigorously validated (the same is not true for the tolerance criterion).This can be done by performing a closure test (68), where one assumes that the underlying parton distributions are known by using a specific (previously determined) PDF set to generate artificial data.Fits are then performed to check whether the output PDFs provide consistent and unbiased estimators of the underlying truth, and that the confidence levels reproduce the correct coverage, among other important considerations.
Techniques have also been developed to convert Hessian PDF sets to Monte Carlo PDF sets and vice versa (74), and to optimize the number of replicas in a PDF ensemble without statistical loss (75).Such methods have led to an extension of the reweighting procedure to Hessian sets (76), the construction of statistical combinations of different PDF sets (77), to optimized versions of PDF sets (78,79), and to statistical tools that analyze the sensitivity of the data set to PDFs (and to any quantity depending on them (80).Finally, we shall note that additional optimization algorithms, not based on least-squares, exist to determine the probability density P(a|D) of the PDF parameters.Among these are Nested Sampling and Markov chain Monte Carlo methods, which have been pioneered in (81).
Procedural uncertainties.There are three classes of procedural uncertainties, all connected to optimization and naturally incorporated in P(a|D).The first is the methodological uncertainty related to the choices made in a fit, namely the basis functions, the functional form, the number of parameters, and the minimization strategy.The second is the extrapolation uncertainty related to the fact that data points, even when infinitely precise, are not covering the entire phase space.The third is the functional uncertainty related to the fact that a set of functions (which are infinite dimensional objects) is determined from a finite amount of information.Within the Hessian method, the three classes of uncertainties are usually assessed by comparing different fits obtained with varied procedural input, and accounted for by the tolerance criterion.In the Monte Carlo method, the methodological uncertainty can be made immaterial by tuning the fitting procedure.Such a tuning is achieved through a level-0 closure test, in which perfect data (no uncertainties) are generated from an assumed underlying law.The test is successful if a perfect fit to the data can be produced, i.e. a fit with vanishing χ 2 .The extrapolation and functional uncertainties are irreducible, yet they can be completely characterized.The former can be determined by looking at the remaining uncertainty on the final PDF in a level-0 closure test.The latter can be determined in a level-1 closure test, where the data replicas are fluctuated to mimic statistical noise (but not the uncertainty of the data itself).The functional uncertainty then arises from all of the statistically equivalent fits.In QCD analyses, closure tests were realized completely for unpolarized PDFs (68) and for nuclear PDFs (41).
Theoretical uncertainties.There is a variety of theoretical uncertainties related to the assumptions made in the computation of hadronic observables that are more challenging to propagate into PDF uncertainties.These include missing higher-order uncertainties (MHOU) due to the truncation of the perturbative expansion to a given order, uncertainties in the input values of the physical parameters such as the strong coupling and the heavy quark masses, uncertainties due to the neglect of power-suppressed corrections should they not be included in factorization formulas, and uncertainties due to nuclear corrections when data on nuclear targets are used in fits of proton PDFs.The first is usually estimated via scale variations, see Eq. 13..The second is accounted for by performing different fits with varied values of the input physical parameters and combining the results.In this respect, the Hessian method allows these uncertainties to be treated as additional sources of nuisance by summing them in quadrature in Eq. 24..The third can usually be kept under control by removing data points from the fit that are particularly sensitive to power suppressed corrections and by looking at the stability of the fit upon variations of this cut-off (82,83).Alternatively, power corrections can be modeled and fitted along with PDFs (84), as well as nuclear corrections (72,84,85), and uncertainties estimated from PDF variations.
A general procedure to represent theory uncertainties in PDFs has been proposed recently in the framework of the Monte Carlo method (86).Assuming that they are Gaussian, it follows from Bayes' theorem that they can be included by redefining the covariance matrix in Eq. 20. as the sum of an experimental and a theoretical part, covij = cov exp ij + cov th ij .The resulting χ 2 is then used both in the sampling of the data replicas and in the optimization of the fit.The problem of propagating theory uncertainties into PDFs is therefore reduced to estimating the theoretical covariance matrix cov th ij by way of an educated guess.Estimates were formulated so far in two cases of unpolarized PDFs studies.First, for nuclear uncertainties by defining the matrix elements of cov th ij as the difference between theoretical predictions obtained either with a free proton or a nuclear PDF, and then taking an average over replicas (87).Secondly, for MHOU (at NLO) by defining the matrix elements of cov th ij as the difference between theoretical predictions obtained with either central or varied factorization and renormalization scales according to various prescriptions (88).In these studies, correlations across data points induced either by the nuclear target in the first case, or by the structure of higher-order corrections in the partonic cross sections and splitting functions in the second, were accounted for properly.Furthermore, nuclear uncertainties and MHOU were validated against the exact nuclear and NNLO results, respectively.The inclusion of such theoretical uncertainties improves the description of the data, shifts the central value of the PDFs towards the truth, and slightly increases their uncertainties.

STATE-OF-THE-ART RESULTS
As is evident by the previous sections, the determination of unpolarized, polarized, and nuclear PDFs is a particularly involved problem.It is therefore addressed by various collaborations of physicists who regularly produce and update general-purpose PDF sets, many of which have a history as long as two decades (see Sect. 1 in (4) and Sect.2.1 in (1) for an overview).Furthermore, while most collaborations perform their global QCD analyses privately, the xFitter collaboration has developed an open-source fitting framework (89).Most of the recent PDF determinations are summarized with their theoretical, experimental, and methodological features in Table 2.All but LSS15, DSSV14 and JAM17, are publicly available through the LHAPDF library (90) and can be readily visualized on-line (91).In this section, we delineate the current status of unpolarized, polarized, and nuclear PDFs using recent PDF extractions.

Unpolarized PDFs
The most recent global determinations of unpolarized PDFs are CT18 (92), MMHT14 (47), NNPDF3.1 (50), JAM19 (93) and ABMP16 (53).Since their publication, the MMHT14 and the NNPDF3.1 analyses have been updated with new data and improved theoretical frameworks.In particular, the first was extended to include the HERA I-II legacy mea-  surements (98), as well as differential measurements in top-pair production (99) and recent jet production measurements (100) from the LHC.It was also extended to a simultaneous determination of αs (101) and of the photon PDF (102).The NNPDF3.1 analysis has since included direct photon (103) and single-top production (104) measurements at the LHC.Similar to MMHT14, it also extended to a determination of αs (105) and of the photon PDF (55).Lastly, the NNPDF collaboration recently assessed theoretical uncertainties from nuclear corrections (87) and missing higher orders (88), as well as studied the impact of small-x resummation (106).
The differences among these PDF sets are summarized in Table 2 except for ABMP16, which is the only unpolarized PDF set determined in the FFN scheme for n f = 3, 4, 5 active flavors.In this case, αs and the heavy quark masses were also free parameters determined together with the PDFs, and each of the six PDF flavors were parameterized in terms of a simple polynomial for a total of 25 free parameters.In addition, ABMP16 utilized the Hessian method for error propagation, with ∆χ 2 = 1.All of the different PDF sets are based on a fairly similar dataset and are available at NLO and NNLO except for JAM19, which is the only PDF set that achieved a simultaneous NLO determination of unpolarized proton PDFs and FFs from charged hadron production in SIDIS and electron-positron annihilation.Additional unpolarized PDF sets exist in the literature, namely JR14 (107), CJ15 (84) and HERAPDF2.0(108), although these PDF sets are based on a reduced set of measurements and are somewhat more limited in scope.In particular, HERAPDF2.0was an analysis of only the HERA I-II legacy data, JR14 studied the impact of a valence-like input below Q 2 ∼ 1 GeV 2 (see also (109)), and CJ15 focused on the high-x data region via the assessment of power-suppressed corrections.
Representative unpolarized PDFs from NNPDF3.1 are presented as a function of x at Q 2 = 100 GeV 2 in the left panel of Fig. 2 for all flavors, and in Fig. 3 for the up quark and gluon distributions normalized to the NNPDF3.1 result.Also illustrated in Fig. 3 is the sea quark ratio Rs = (s + s)/(ū + d).Results correspond to the NNLO PDF sets except for JAM, which is displayed at NLO, and all uncertainty bands correspond to 68% confidence levels.The general features of unpolarized PDFs can be summarized as follows.
Valence quarks.The relevance of quark valence distributions is twofold.First, they give the global properties of the nucleon across their entire range in x, such as its charge and baryon number.Furthermore, at large x they carry most of the nucleon's momentum and thus are sensitive to its non-perturbative structure and to the production of new particles, e.g.heavy W and Z bosons at high rapidities and invariant masses.Up and down valence distributions are primarily constrained by collider DY measurements at small x (in particular by the W ± asymmetry), by neutral-and charged-current DIS measurements at medium x, and by fixed-target DY measurements at medium to high x.Results obtained from various PDF sets agree very well for the up valence distribution, which has an overall relative uncertainty of a few percent in the data region, while they are a little more widespread for the down valence distribution (see Sect. 6.2 in (1) for more details).
Sea quarks.Anti-up and anti-down sea quarks are probed by parity-violating processes, in particular by charged-current DIS and DY measurements.Unlike the valence distributions, they are strongly suppressed at large x and show a steep rise at small x due to being produced mainly by gluon splitting.Anti-up and anti-down PDFs display larger uncertainties than their valence distribution counterparts, with more marked differences across PDF sets.
Concerning the strange quark PDF, discrepancies have arisen in analyses of inclusive W ± production from ATLAS with respect to charged current neutrino DIS measurements.While the former support a ratio of the strange to non-strange light sea distributions around unity, Rs ∼ 1.13 at x = 0.023 and scale Q 2 = 1.9 GeV 2 , the latter give a result Rs ∼ 0.5 at similar kinematics (50).The tension may be somewhat relieved if a flexible parameterization is used (50), if massive corrections to neutrino DIS structure functions are included in the analysis, or if experimental correlations are relaxed (110).The picture is further complicated by recent analyses of charged kaon production in SIDIS, which can provide crucial strange quark PDF information due to the enhancement from the favored s(s) → K − (K + ) FF.Using a reweighting procedure, a minor suppression of Rs was found in the intermediate x region compared to the result obtained without the inclusion of SIDIS data (111).On the other hand, the JAM19 analysis found a more suppressed Rs result by simultaneously fitting the FFs with the proton PDFs (93).All of the various Rs distributions can be seen in Fig. 3, where uncertainties are rather large.To alleviate this, better understanding of the various systematic uncertainties that affect the data (especially nuclear modifications in neutrinonucleus DIS) and clean processes that are sensitive to strangeness are desperately needed.One such process is W -boson production in association with a charm quark.However, in this case important NNLO corrections to the theoretical prediction remain absent.
Heavy quarks.Heavy quarks produce sizable contributions to inclusive DIS structure functions, typically about 30% (3%) in F2 for the charm (bottom).They are usually determined perturbatively through gluon splitting in quark-antiquark pairs.Under this assumption, the distributions turn out to be consistent across different PDF sets, with a precision similar to that of the gluon PDF (see below).An intrinsic, non-perturbative heavy quark component can also be assumed, however.In this case, heavy quark PDFs must be parameterized and determined together with the light quark PDFs.This approach is default for the treatment of the charm in NNPDF3.1, which leads to a reduced dependence on the charm mass for several high-energy benchmark cross sections, to an overall improvement in the description of the data, and to a more stable gluon PDF (see also Sect.6.2 for non-perturbative implications).
Gluon.The gluon PDF is probed in various ways.At small to medium x, it can be accessed through measurements of the structure function FL in inclusive DIS as well as in scaling violations in DGLAP evolution.At medium to large x, measurements of the Z-boson transverse momentum distributions can provide information.Lastly at large x, the gluon can be constrained by jet and top-pair production measurements.The available measurements for the last two classes of processes have recently raised some controversy, where it was argued that they can be well described only if experimental correlations are relaxed.For jets, such decorrelations do not affect the ensuing PDFs (92, 100), while they might for top production (92,99), especially in relationship to the specific combination of kinematic distributions included in the fit.For this reason, the comparison of various PDF sets displays different shapes for the gluon PDF in Fig. 3, particularly at large x, albeit with large uncertainties.This behavior may challenge the study of high-invariant mass states beyond the SM, possibly in conjunction with threshold resummation (112).At very small x, the gluon remains largely unconstrained; however, small-x resummation was recently demonstrated to point towards evidence of BFKL dynamics (106).
Photon.It has been recently shown that the photon PDF, relevant in several electroweak processes through photon-initiated and higher-order contributions, can be completely determined from the structure functions of the proton F2 and FL (36).The computation, known as LUXqed formalism, made previous approaches based on a model assumption (113) and a loosely constrained parameterization (114) obsolete.Recent analyses (55,102) supplemented the LUXqed formalism with constraints from the data.In these cases, the photon PDF resulted in an uncertainty of few percent, carrying about 0.5% of the proton's momentum, and accounted for corrections of up to 20% in DY, vector-boson, top quark and Higgs production processes at the LHC (see Sect. 7 in (1) for additional details).

Polarized PDFs
The most recent analyses of polarized PDFs are LSS15 (94), DSSV14 (115), NNPDF-pol1.1 (96) and JAM17 (39).Since publication, the DSSV14 analysis has been updated with a Monte Carlo variant (95), which studied the impact of recent di-jet measurements from STAR (116).The NNPDFpol1.1 PDF set was also updated with STAR measurements, including the same di-jet measurements and additional W -boson production data (117).
The similarities and differences among the polarized PDF sets are summarized in Table 2. Here, all of the global QCD fits are performed at NLO accuracy and are more or less consistent in their theoretical and methodological choices.However, each analysis supplements inclusive DIS measurements with different data.The JAM17 and DSSV14 PDF extractions incorporated spin asymmetry measurements from pion and kaon production in SIDIS by HERMES (118) and COMPASS (119,120).Single-jet production asymmetries in polarized proton-proton collisions measured by STAR (121) were used in the DSSV14 and NNPDFpol1.1 analyses.The DSSV14 dataset also included neutral pion production in polarized proton-proton collision measurements by PHENIX (122), while the NNPDFpol1.1 data contained open-charm production asymmetries in DIS from COMPASS (123) and W -boson production asymmetries in polarized proton-proton collisions from STAR (124).Lastly, while not indicated in Table 2, the JAM17 dataset included charged pion and kaon production from single-inclusive electron-positron annihilation to help constrain the final state FFs simultaneously with the polarized PDFs (39).
Representative examples of polarized PDFs at NLO are illustrated as a function of Polarized PDFs as a function of x for the total strange (left), anti-up (middle) and anti-down (right) PDFs at Q 2 = 10 GeV 2 .Displayed are the JAM17 (red), DSSV14 (blue), and NNPDFpol1.1 (green) NLO PDF sets.Uncertainty bands are computed as 68% confidence levels.
x in the central panel of Fig. 2 for all flavors at Q 2 = 100 GeV 2 , and in Fig. 4 for the total strange (x∆s + = x(∆s + ∆s)), anti-up (x∆ū), and anti-down (x∆ d) distributions at Q 2 = 10 GeV 2 .In the latter figure, results are compared between the JAM17, DSSV14, and NNPDFpol1.1 PDF sets.Moreover, all uncertainty bands correspond to 68% confidence levels.The general features of polarized PDFs can be summarized as follows.
Valence quarks.The up quark valence distribution is the most constrained polarized PDF, primarily due to measurements of the proton structure function g1 over a relatively broad range of x and Q 2 .The corresponding down quark distribution, which is opposite in sign, is smaller in magnitude and shows somewhat larger uncertainties.Nevertheless, a fair agreement between various PDF sets has been achieved for the valence polarizations, with differences originating from theoretical, experimental, and methodological choices generally being smaller than the PDF uncertainty, i.e. the uncertainty of the data.

Sea quarks.
As can be seen in Fig. 2, the polarization of the sea quark is significantly smaller than the polarization of the valence quarks.It is also more dependent on the measurements included in each PDF set and on the flavor assumptions with which the data are analyzed.While all PDF sets show a fairly consistent anti-down quark polarization, the anti-up distribution has opposite sign for the DSSV14 and NNPDFpol1.1 sets in the region 0.1 x 0.5.Recall, however, that this is driven by different processes in the two PDF sets, particularly by fixed-target SIDIS multiplicities and by collider W -boson production spin asymmetries, respectively.The former requires the knowledge of FFs 8., which are taken as a fixed input in the DSSV14 analysis.The difference is somewhat relieved by the JAM17 result, where the FFs entering the analysis of the SIDIS data were determined together with the polarized PDFs.
The situation is even more involved in the case of the polarized strange PDF.Here, the distribution is entirely unconstrained in PDF sets that do not include SIDIS asymmetries from the kaon sector unless one imposes the SU f (3) constraint from weak baryon decays (Eq.18.).This is the case, for instance, in the NNPDFpol1.1 analysis, which displays a negative strange polarization peaked at x ∼ 0.1 as a result.In contrast, the DSSV14 PDF set, based on SIDIS data with a fixed kaon FF, finds a sign-changing strange PDF that is incompatible with the NNPDFpol1.1 result in the range 0.02 x 0.2.This discrepancy is Nuclear modification ratios (Eq.11.) as a function of x at Q 2 = 10 GeV 2 for the gluon PDF in again somewhat alleviated within the large uncertainties of the JAM17 result, obtained by removing the SU f (3) constraint from the analysis and by fitting the kaon FFs simultaneously with the polarized PDFs.In any case, more precise SIDIS polarization asymmetries for kaon production are desirable to reduce the uncertainty of the strange PDF and shed light on possible SU f (3) symmetry breaking effects.Moreover, separate measurements for positively and negatively charged kaons could in principle allow for a separation of the strange and anti-strange quark components for the first time.
Gluon.The polarized gluon distribution has been elusive for a long time, since it is only weakly constrained by DIS and SIDIS measurements, where it enters as a higher-order correction, and by the limited Q 2 coverage of the data, in which sensitivity can come from the DGLAP evolution equations.In fact, it was believed to be small until very recently.The availability of precise jet, di-jet, and hadron production spin asymmetry measurements at RHIC were a game changer in this respect, revealing for the first time a sizable polarization of the gluon PDF in the DSSV14 and in the NNPDFpol1.1 analyses.Such evidence, however, is limited to the region 0.02 x 0.4, outside of which the polarized gluon PDF is affected by large extrapolation uncertainties that prevent any definitive conclusion about its role in understanding the proton spin decomposition (see Sect 6.2 and Ref. ( 125)).

Nuclear PDFs
In the nuclear sector, the most recent PDF determinations are nCTEQ15 (52), EPPS16 (42), nNNPDF1.0(41) and TUJU19 (97).The nCTEQ15 analysis has since been updated with measurements of vector boson production (126) in proton-lead and lead-lead collisions, while the EPPS collaboration recently studied the impact of di-jet (127) and D-meson production (128) in proton-lead collisions.
The differences among these PDF sets are summarized in Table 2. Concerning the data set, they are all based on inclusive DIS data except for neutrino DIS data in the case of nNNPDF1.0.Complementary measurements are added only in the nCTEQ15 and EPPS16 global analyses, where various fixed-target and collider DY observables are included with jet and hadron production data.Finally, all of the PDF sets are determined at NLO, with the DIS-only nNNPDF1.0 and TUJU19 analyses extending also to NNLO.
In the right panel of Fig. 2, an example nuclear PDF set for lead is displayed at GeV 2 as a function of x for all flavors.Furthermore, we illustrate in Fig. 5 the ratio R A f (x, Q 2 ) 11. at Q 2 = 10 GeV 2 for the gluon distribution in a light nucleus, 12  6 C, and a heavy nucleus, 208  82 Pb.Also shown is the quark singlet PDF combination for 208 82 Pb.Here the nNNPDF1.0,EPPS16, and TUJU19 NLO PDF sets are compared with uncertainty bands corresponding to 90% confidence levels.The ratio R A f (x, Q 2 ) presented in Fig. 5 provides the most relevant information in nuclear PDF studies since it sheds light on the various nuclear modifications with respect to the free proton in different x regions.In particular, experimental data have indicated a shadowing behavior (R < 1) at low x 0.1, followed by anti-shadowing (R > 1) in the x ∼ 0.1 region, and lastly the EMC effect (R < 1) in the valence region, x ∼ 0.4, the latter being named after the experimental group from which it was discovered (129).While the mechanisms that generate these effects are still under investigation, they are expected to appear as general features in the nuclear PDF results, which are summarized as follows.
Quarks.In Fig. 2, the behavior of the valence and sea quark PDFs mimic the distributions of their proton counterparts, albeit they are affected by larger uncertainties which tend to increase as the atomic number A increases.The reason for this can be traced back to the quantity and quality of the data, which is significantly more limited in the variety of processes and in the breadth of the kinematic coverage for nuclear PDFs than for proton PDFs (refer to Fig. 1).Note that, opposite to the proton case, the up and down quark valence distributions are fairly similar as a result of the nuclear PDF definition (Eq.10.).Assuming exact isospin symmetry, the visible differences then arise from the non-isoscalarity of the nucleus.In general, the agreement between the different nuclear PDF sets is excellent for the total quark singlet distribution presented in Fig. 5, where the shape clearly displays shadowing, anti-shadowing, and EMC effects in their expected x regions.Additional comparisons of individual quark and antiquark flavors are not made here due to the lack of flavor separation in DIS data used in the nNNPDF1.0 and TUJU19 analyses, but can be found in Refs.(42,52).Lastly, despite containing similar datasets, the nNNPDF1.0distributions show significantly larger uncertainties than the TUJU19 result, especially in the small-x extrapolation region.While this is partly due to the different uncertainty estimation techniques, it can also be attributed to the use of a more flexible parameterization in nNNPDF1.0(see Refs. (41) and (97) for more details).
Gluon.The features of the gluon distribution in Fig. 2 are also similar to those in the proton case, apart from having larger uncertainties.However, the nuclear modifications to the gluon distribution in Fig. 5 are much less distinct than those for the quark singlet.While the TUJU19 PDF set displays very clear shadowing behavior for the ratio, the larger uncertainties for the EPPS16 and nNNPDF1.0distributions prevent any definitive conclusion about nuclear shadowing in the low x region.Furthermore, gluon anti-shadowing seems to appear in different regions of x for the different PDF analyses.For EPPS16, the gluon exhibits anti-shadowing in a similar region as the quark singlet, while the TUJU19 and nNNPDF1.0analyses display an anti-shadowing effect at larger x.The latter is likely an artefact of the fitting procedure rather than an actual physical effect, as the two analyses do not include data sensitive to the gluon.However, further investigation is required to make any conclusive statements.Lastly, there is a noticeable difference in uncertainties between the gluon uncertainties in carbon compared to that of heavier lead nuclei.In fact, this is also true for the quark distributions and arises by imposing the boundary condition at A = 1 (see Sect. 2.3).In this sense, the proton PDFs can play a crucial role in constraining the nuclear distributions for light nuclei.

PRESENT AND FUTURE RELEVANCE OF PDFS
While significant progress has been achieved in the determination of PDFs, various issues remain open.In this section we give a concise summary of the ones that we find more relevant in particle and hadronic physics, and how future colliders will aid in solving them.

Precision physics in the SM and beyond
Since the discovery of the Higgs boson (130), to which an accurate determination of unpolarized PDFs was instrumental, the LHC has allowed the SM to be established with unprecedented precision.Nevertheless, PDFs remain one of the main sources of uncertainty in pushing forward precision and discovery physics at the LHC.
Determination of the strong coupling.It is customary to determine the value of the strong coupling αs from a variety of processes involving hadrons in the initial state, which require knowledge of PDFs.In this respect, determinations of αs have historically been of two kinds (see (131) for a review).The first realizes a simultaneous fit of both αs and the PDFs in a global QCD analysis from a (more or less) wide set of data and processes; the second finds the likelihood of some new data for a single process as a function of αs, based on a pool of pre-existing PDF sets determined with various input values of the strong coupling.While it has been recently argued that only an appropriate implementation of the first method, which takes into account correlations between αs and the PDFs, leads to an unbiased result (132), in both cases the PDF uncertainty becomes relevant.For instance, at NNLO and at the Z-boson mass mZ , the analysis of Ref. (105) found αs(mZ ) = 0.1185 ± 0.0005 exp ± 0.0001 meth ± 0.0011 th , where exp, meth, and th uncertainties are the PDF, methodological and missing higher order theoretical uncertainties, respectively.
Determination of SM parameters.Parton distributions represent one of the dominant sources of theoretical uncertainty for the determination of the Higgs boson couplings and cross sections at the LHC (see (133) and references therein for a review).Likewise, they are the most significant limiting factor in the precise determination of electroweak parameters such as the mass of the W -boson, the mass of the top-quark, and the CKM matrix elements (see (134) and references therein, and (135) for a dedicated study).
The search for BSM physics.Beside precise benchmarks of the SM, extensive searches are in progress at the LHC for physics signatures beyond the SM.Since the majority of measurements at the LHC and elsewhere have so far been consistent with SM predictions, many have turned towards analyzing LHC data in the context of a SM effective field theory (SMEFT), where new physics signals may be found in the dynamics of SM fields via higher (mass) dimension operators (see Ref. (136) for a review).A large number of SMEFT analyses have been performed in both the electroweak and Higgs sectors, as well as the top quark sector, to constrain the Wilson coefficients related to the higher mass dimension terms.Here the PDFs play an essential role in the calculations of the SM part of the SMEFT expansion (see e.g.(137) and references therein).In principle, however, effects beyond the SM could be absorbed by the PDFs themselves in a global QCD analysis.To avoid this, one should determine a set of SMEFT coefficients simultaneously with the PDFs.Despite being a daunting task, it was demonstrated to be achievable in a simple case, whereby a small subset of SMEFT coefficients was determined together with the PDFs by analyzing inclusive DIS data (138).In that case, no significant indication of physics beyond the SM were found.However, much work still remains in the field of SMEFT to obtain better constraints on the Wilson coefficients and identify new physics signals.As higher order QCD and electroweak calculations are made available in the SMEFT, and experimental measurements increase in abundance and precision, the role of PDFs in beyond the SM searches will become ever more relevant.

A window to non-perturbative behavior
Because of their nature, PDFs can shed light on the non-perturbative aspects of nucleons and nuclei at perturbative scales.Among these, some of the most compelling are as follows.
Sea quark asymmetries.Perturbative QCD predicts light sea quarks to be equally produced by gluons in virtual fermion loops.However, a sizable asymmetry in the unpolarized anti-up and anti-down quark sea was observed in DY by the NA51 (139) experiment and confirmed by the NuSea/E866 (140) experiment at CERN and Fermilab, respectively.Such an asymmetry therefore has a non-perturbative origin, which is commonly explained in the context of a meson cloud model (141).Significant effort is underway to further constrain the d − ū unpolarized PDF asymmetry and also to clarify its origin, see e.g.Ref. (142) for a recent study of the meson cloud applied to the ∆ baryon.Experimentally, new DY measurements from the SeaQuest/E-906 experiment at Fermilab (143) are soon expected to be finalized.Significant interest in flavor symmetry breaking effects exists also for other sea quark distributions, specifically the unpolarized (polarized) s − s (∆s − ∆s) and polarized ∆ū − ∆ d asymmetries.The former relates to constraining the strange quark PDF uncertainties, as discussed in Sect. 5.The latter can be fairly well constrained by hadron multiplicities in SIDIS and by W ± single-spin asymmetries in polarized proton-proton collisions.Interestingly, these have shown evidence of a positive ∆ū − ∆ d asymmetry in the intermediate x region, opposite to that generated in the unpolarized case (117).Down to up quark ratio at large x.Quarks carrying large momentum fractions can provide additional insight into non-perturbative quark-gluon interactions.Despite progress in PDF determinations, the d/u and ∆d/∆u ratios are affected by large uncertainties that prevent any conclusive comparison with theoretical models (46,144).The reason for this state of affairs is the lack of experimental data with free neutron targets (or beams).Instead, light nuclei such as deuterium or helium are commonly used in (DIS) experiments, and therefore bound nucleon effects must be taken into account in PDF analyses (see e.g.Ref. (145) and references therein).The CJ collaboration recently focused on the determination of the u and d PDFs at large x in an analysis where nuclear effects in electron-nuclei scattering were treated by means of a weak binding approximation and off-shell corrections (84).On the experimental front, the 12 GeV program at JLab ( 146) is well underway with several dedicated experiments aiming to reveal large-x PDF dynamics.In particular, the MARATHON experiment has recently completed measurements of unpolarized DIS from a tritium target, which aims to extract the neutron to proton structure function ratio up to x ∼ 0.85 (147).Results are expected to be reported soon, which can subsequently be used in a global QCD analysis to further reduce the d/u ratio uncertainty.
Non-perturbative charm.Apart from stabilizing predictions for a class of unpolarized hadronic cross sections at high energies, non-perturbative charm in the unpolarized proton is also relevant at low energies.Its role is primarily discussed in the context of a non-zero 5-quark Fock state component of the proton wave function, |uudcc (148).Phenomenologically, world scattering data sensitive to the unpolarized charm quark PDF, most notably those measuring the F c 2 structure function, can be analyzed by separating the charm component into a perturbative and a non-perturbative piece, the latter of which carries some fraction of the proton's momentum, x IC ≡ 1 0 dxx [c(x) + c(x)] = 0. Various phenomenological analyses suggest that this component is small but sizable.The CTEQ-TEA collaboration found x IC 2.5% at Q = 1.65 GeV using various non-perturbative models (149), while the NNPDF collaboration found x IC = (1.6 ± 1.2)% at the same energy scale, in a fit where the charm PDF was parameterized together with the light quark PDFs (54).On the other hand, the JR collaboration found x IC = (0.15 ± 0.09)%, utilizing data at low final-state invariant-mass energies from SLAC (150).While the NNPDF result is compatible with the other two within its large uncertainties, it corresponds to a completely different shape of the charm PDF.Such differences are likely to originate from the various theoretical and methodological choices that are inherent in each analysis.
The proton spin.The size of the contribution of quarks, antiquarks and gluons to the nucleon spin is quantified by the first moments of the corresponding polarized PDFs, according to the canonical decomposition of the proton's total angular momentum (151).Both the NNPDFpol1.1 and the DSSV14 analyses agree on the values of the quark singlet and gluon polarized PDF first moments in the data region, which are found to be 0.23-0.30and 0.20-0.23 at Q 2 = 10 GeV 2 , respectively.However, uncertainties affecting the smallx extrapolation region, as well as those surrounding SU f symmetry constraints, prevent an assessment of the residual contribution to the total angular momentum of the proton coming from the orbital angular momentum of quarks and gluons.

Parton dynamics in the nuclear medium
The relevance of nuclear parton distributions is perhaps the most far-reaching of the collinear distributions, especially in the following domains.
Nuclear modifications and the characterization of quark-gluon plasma.Nuclear PDFs are key to understanding how parton dynamics are modified in the nuclear medium and the mechanisms that generate the shadowing, anti-shadowing, and EMC effects (152).In addition, they play an essential role in the characterization of the quark-gluon plasma, the hot and dense medium created in the early Universe currently being studied by heavy ion collisions at RHIC and LHC (153).In this respect, measurements of proton-ion and ion-ion collisions will help constrain the nuclear PDFs and their uncertainties.In addition to those outlined in Sect.3, the LHCb D-meson production data were demonstrated to constrain the nuclear gluon PDF at low x (128).More recently, measurements of J/ψ production in ultra-peripheral collisions of two lead nuclei were made available by the ALICE experiment at the LHC (154); similar measurements of gold nuclei are expected to be made by the STAR experiment at RHIC (155).Such collisions, where the impact parameter is larger than the combined radii of the colliding nuclei, provide a rather clean way to study photon-nucleus interactions.While coherent photoproduction of vector mesons can provide auxillary constraints on nuclear PDFs in general, J/ψ is of special interest since it is a gluon initiated process and therefore its associated theoretical prediction is sensitive to the nuclear gluon PDF (156).Once included in a PDF determination, the data can potentially provide additional evidence of shadowing and gluon saturation at small x values.
Relationship with proton PDFs.A significant number of observables used in global QCD analyses of free-proton PDFs involve proton-nucleus or lepton-nucleus scattering, namely proton-deuteron (and proton-ion) DY measurements, and charged-current neutrino-nucleus measurements.Being essential to free-proton PDF extractions, in particular to constrain the sea quark PDFs, nuclear corrections must be carefully taken into account using a nuclear model or parameterized and determined from a fit to the data.Whatever approach is adopted, effects from nuclei will consequently increase the uncertainty in proton PDFs since nuclear corrections are known less precisely.A methodology to include such an uncertainty in proton fits was developed in (87), as was discussed in Sect.4.3.
Implications for astroparticle physics.Neutrino telescopes, such as IceCube (157), KM3NeT (158) and Baikal-GVD (159), have been developed in recent years to measure ultra high energy neutrino fluxes (with energy Eν 10 7 GeV) as a way to study, in conjunction with cosmic ray observatories such as Auger (160), sources of cosmic rays in the Universe and QCD at multi-TeV scales, far outside the energy scales probed by available colliders.Theoretical predictions of neutrino-nucleus scattering cross sections can be made in the framework of perturbative QCD for both signal (161) and backround (162) events.Since the detectors measure scattering events from a target source material consisting of nuclei, usually water molecules, a precise knowledge of nuclear PDFs is necessary.Nuclear modifications induced by the medium increase the cross sections by 1-2% at energies below 100 TeV (antishadowing), and reduce it by 3-4% at higher energies (shadowing) (163).Nuclear effects also alter the inelasticity distribution of neutrino interactions in water/ice by increasing the number of low inelasticity interactions, with a larger effect for neutrinos than antineutrinos.These effects are particularly important in the energy range below a few TeV.An accurate control of nuclear PDF uncertainties is pivotal to characterize all of these effects, in particular above the IceCube energy reach, still within the larger reach of Auger (Eν ∼ 10 12 GeV), where the absence of experiment constraints remains the limiting factor to determine the composition of the highest energy hadronic cosmic rays (164).

The role of future colliders
Several new accelerator upgrades and designs are being planned that will improve much of the PDF statuses discussed in the above sections.The LHC is currently upgrading to its high luminosity (HL) phase (165), where various observables will be measured with significantly increased statistical precision.This will extend the LHC physics output for another decade, providing rich precision measurements of SM parameters such as the Higgs couplings and the W -boson mass (133,134).At lower energies, the JLab has just completed its upgrade to 12 GeV (146), which will allow for a careful investigation of aspects related to the non-perturbative structure of the unpolarized and polarized proton.
Beside these upgrades, there are at least two additional accelerator designs being developed that will be particularly beneficial for PDF studies.The first is an electron-proton collider at the LHC (LHeC), which intends to inject roughly 60 GeV electron beams into the LHC proton collider (166).With the point-like nature of the electron probe and the energies of the LHC, high precision DIS measurements will be achieved for Bjorken-x down to ∼ 5 × 10 −6 , a kinematic region that has yet to be explored.The second is an electron-ion collider (EIC) planned to be built on the current site of RHIC, with medium energies of up to 90 GeV (167).The aim of the EIC is to study with high precision the three-dimensional structure of nucleons and nuclei.In terms of collinear structure, however, the new accelera-tor will be able to provide high precision DIS measurements for nuclear PDF determination down to x ∼ 10 −4 .One of the physics aims, in fact, is to provide such measurements in order to pin down the nuclear gluon distribution and reveal gluon saturation effects, similarly to what HERA did for the proton.It is also a machine that can polarize the beams, thus providing important constraints for the polarized PDFs.
The extended kinematic regions attained by both the LHeC and the EIC are displayed as dashed regions in Fig. 1.Impact studies of projected pseudodata foreseen at these facilities have been extensively performed for unpolarized PDFs in Ref. (168,169), for polarized PDFs in Refs.(170,171), and for nuclear PDFs in Refs.(41,172,173).

OUTLOOK
Global QCD analyses of PDFs continue to be an active field in particle and hadronic physics.The upgrade of the LHC, combined with the advent of new colliders, will enhance their relevance even further in the next decades.Now that PDFs are entering a new precision era, they will be critical for calculations of signal and background events in physics searches beyond the SM and in bringing electroweak symmetry breaking under complete control.They will allow us to explore nuclear modifications, the onset of gluon saturation, and aid in understanding high-energy neutrino interactions in astroparticle physics.And they will shed light on the non-perturbative structure of the proton and unravel the proton's spin decomposition.To achieve such ambitious goals, we foresee more sophisticated PDF determinations in the future, all of which should fulfill the following criteria.
1. PDF determinations should be as global as possible, i.e. the range of the input dataset should be extended to cover unexplored kinematic regions and new processes.In this respect, experimental uncertainties should be carefully scrutinized, in particular as correlated systematic uncertainties dominate collider measurements.The observed tensions between data sets, and the different parton sensitivity to different observables across PDF extractions must be carefully understood with appropriate statistical tools.Benchmarking (and, if beneficial, combination) of different PDF sets, including polarized and nuclear PDFs, must become standard.2. Global QCD analyses should be as accurate as possible, i.e. theoretical calculations should be performed at the highest available order in perturbative QCD.This is currently NNLO for unpolarized and nuclear PDFs, and NLO for polarized PDFs.Approximate higher-order PDF fits can actually be performed only to inclusive DIS data.In this respect, a N 3 LO determination of unpolarized PDFs is desirable to study the effect of using a NNLO matching condition, and whether it leads to a charm PDF closer to the fitted one.A NNLO determination of polarized PDFs, instead, will allow us to check the perturbative stability of all of the current polarized PDF sets.Unpolarized PDF determinations should also resum large logarithms: large-x resummation is particularly important for searches of new physics and matching of fixed-order calculations to parton showering in Monte Carlo event generators, while small-x resummation is relevant to reveal the onset of BFKL dynamics and gluon saturation.Finally, the treatment of heavy quark PDFs on the same footing as the light quark distributions, and also the systematic inclusion of electroweak corrections, should become standard for unpolarized PDFs.
3. PDF determinations should represent uncertainties as faithfully as possible, i.e. their statistical meaning should be carefully validated.Furthermore, a complete characterization of theoretical uncertainties should be implemented.This includes the uncertainty from input physical parameters and from the neglect of missing higher orders or other corrections in the calculation of hadronic observables (such as nuclear and higher-twist corrections).In view of related needs for precision physics, this will become mandatory for unpolarized PDFs. 4. PDF analyses should aim at becoming simultaneous, i.e. parameters from different non-perturbative aspects of the theory should be determined at the same time.The simplest example is represented by the determination of the strong coupling along with the PDFs.More ambitiously, we envision three cases where simultaneous fits will become critical.The first is in the extraction of polarized PDFs, where observables are often presented as ratios of spin-dependent to spin-averaged cross sections.In this case, unpolarized and polarized PDFs should be treated simultaneously in a universal QCD analysis.The second is the analysis of semi-inclusive measurements involving the production of an identified hadron in the final state.We believe that these processes can be analyzed accurately only if the final-state FFs are determined together with the initial-state PDFs.While this is primarily significant for QCD analyses of polarized PDFs, it can also provide additional information on the sea distributions in unpolarized PDF fits.The third is analyses of LHC processes in the SMEFT.We believe that searches for physics beyond the SM can be realized in this framework only if SMEFT Wilson coefficients are extracted simultaneously with the PDF parameters.Because such fits will dramatically increase the complexity of the parameter space, development of new optimization techniques such as those based on machine learning and artificial intelligence will likely be required.5. PDF extractions should benefit from input provided by ab-initio lattice QCD calculations (174) and by advancements in the theoretical and phenomenological understanding of non-perturbative functions that generalize collinear PDFs (175).These fields witnessed spectacular progress in recent years, which cannot be ignored.
As can be seen by this review, significant efforts are underway to achieve many of the objectives listed here.While much work remains, global QCD analyses have undoubtably taken major steps forward to address challenging questions in QCD and deepen our understanding of nucleon and nuclei substructure.

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.

Figure 2 A
Figure 2 A representative snapshot of unpolarized and polarized proton PDFs, and on nuclear (lead) PDFs.Parton distributions are taken from the NNPDF3.1 (NNLO), DSSV14 (NLO) and EPPS16 (NLO) analyses, respectively.Uncertainty bands correspond to Monte Carlo 68% confidence levels for NNPDF3.1 and DSSV14, and to Hessian 90% confidence levels for EPPS16.

Figure 3
Figure 3 Representative examples of unpolarized PDFs as a function of x at Q 2 = 100 GeV 2 for the up quark (left) and gluon (middle) PDFs normalized to the NNPDF3.1 central value, and for the ratio Rs = s+s ū+ d (right).Results correspond to NNLO MMHT14, CT18, NNPDF3.1 and NLO JAM19 PDF sets.Uncertainty bands computed as 68% confidence levels.

12 6
C and 208 82 Pb, and the quark singlet PDF in208  82 Pb.Results are shown for the nNNPDF1.0,EPPS16 and TUJU19 analyses at NLO. Uncertainty bands correspond to 90% confidence levels.