Hadron structure in high-energy collisions

TL;DR

This paper reviews the concept of parton distribution functions (PDFs) in hadrons, focusing on their definitions, how they are fitted to data, and methods to test the validity of the assumptions involved.

Contribution

It provides a comprehensive overview of PDFs in protons and nuclei, detailing the statistical fitting process and validation tests, with insights into handling assumption failures.

Findings

Detailed explanation of PDF definitions and their relation to cross sections

Discussion of the Hessian method for fitting PDFs to data

Proposals for testing and addressing invalid assumptions in PDF analysis

Abstract

Parton distribution functions (PDFs) describe the structure of hadrons as composed of quarks and gluons. They are needed to make predictions for short-distance processes in high-energy collisions and are determined by fitting to cross section data. We review definitions of the PDFs and their relations to high-energy cross sections. We focus on the PDFs in protons, but also discuss PDFs in nuclei. We review in some detail the standard statistical treatment needed to fit the PDFs to data using the Hessian method. We discuss tests that can be used to critically examine whether the assumptions are indeed valid. We also present some ideas of what one can do in the case that the tests indicate that the assumptions fail.

Equations

$$\frac{d\sigma }{d{Q}^{2}d(P_{\mathrm{A}}\cdot q)}=\sum _a\int d\xi f_{a/\mathrm{p}}(\xi )\frac{d\hat{\sigma }}{d{Q}^{2}d(P_{\mathrm{A}}\cdot q)}.$$

$$\sigma [F]\approx \sum _a\int d\xi f_{a/A}(\xi ,{\mu }^2)\hat{\sigma }[F].$$

$$\sigma [F]\approx \sum _{a,b}\iint d{\xi }_{a}d{\xi }_b{f}_{a/A}({\xi }_a,{\mu }^2)f_{b/B}({\xi }_b,{\mu }^2)\hat{\sigma }[F].$$

$${\chi }^2(a)=\sum _{ij}(D_i-T_i(a))(D_j-T_j(a))C_{ij},$$

$$(P^{+},P^{-},{\mathbf{P}}^{\perp })=\left(P^{+},\frac{m_{\mathrm{p}}^2}{2P^{+}},0\right).$$

$$\begin{array}{rl}[b({\xi }^{\prime }& P^{+},{\mathbf{k}}_{\perp }^{\prime },s^{\prime },c^{\prime };i),b^{†}(\xi P^{+},{\mathbf{k}}_{\perp },s,c;i){]}_{+}\\ =& (2\pi {)}^{3}2\xi P^{+}\delta ({\xi }^{\prime }{P}^{+}-\xi P^{+})\delta (k_{\perp }^{\prime }-k_{\perp }){\delta }_{s^{\prime }s}{\delta }_{c^{\prime }c}.\end{array}$$

$$b(\xi P^{+},\bm{k}_{\perp},s,c;i)\big{|}{0}\big{\rangle}=0\;.$$

$$\rho ($$

$$=\frac{1}{(2\pi {)}^{3}2\xi }\sum _{s,c}{b}^{†}(\xi P^{+},{\mathbf{k}}_{\perp },s,c;i)b(\xi P^{+},{\mathbf{k}}_{\perp },s,c;i).$$

$$\int_{{\cal V}_{3}}\!d\xi\,d\bm{k}_{\perp}\ \rho(\xi P^{+},\bm{k}_{\perp};i)\big{|}{\Psi}\big{\rangle}=N({\cal V}_{3})\big{|}{\Psi}\big{\rangle}\;.$$

$$f_{i/\mathrm{p}}^{(0)}(\xi)\big{\langle}{P^{\prime}}\big{|}{P}\big{\rangle}=\int\!d\bm{k}_{\perp}\big{\langle}{P^{\prime}}\big{|}\rho(\xi P^{+},\bm{k}_{\perp};i)\big{|}{P}\big{\rangle}\;.$$

$$\begin{split}P_{\rm dy}&\psi_{i,c}(0,x^{-},\bm{x}_{\perp})={}\frac{1}{(2\pi)^{3}}\int_{0}^{\infty}\!\frac{dk^{+}}{2k^{+}}\int d\bm{k}_{\perp}\\ &\sum_{s}\,\bigg{\{}P_{\rm dy}\,u(k,s)\,e^{-ik\cdot x}b(k^{+},\bm{k}_{\perp},c,s;i)\\ &\qquad+P_{\rm dy}\,v(k,s)\,e^{+ik\cdot x}d^{\dagger}(k^{+},\bm{k}_{\perp},c,s;i)\bigg{\}}\;.\end{split}$$

$$\begin{split}f_{i/\mathrm{p}}^{(0)}&(\xi)\big{\langle}{P^{\prime}}\big{|}{P}\big{\rangle}=\frac{P^{+}}{2\pi}\int\!dy^{-}\,e^{-{\mathrm{i}}\xi P^{+}y^{-}}\int\!dx^{-}\ d\bm{x}_{\perp}\\ &\times\big{\langle}{P^{\prime}}\big{|}\bar{\psi}_{i}(0,x^{-}+y^{-},\bm{x}_{\perp})\gamma^{+}\psi_{i}(0,x^{-},\bm{x}_{\perp})\big{|}{P}\big{\rangle}\;.\end{split}$$

$$\begin{split}\big{\langle}{P^{\prime}}\big{|}\bar{\psi}_{i}(0,&x^{-}+y^{-},\bm{x}_{\perp})\gamma^{+}\psi_{i}(0,x^{-},\bm{x}_{\perp})\big{|}{P}\big{\rangle}\\ ={}&e^{{\mathrm{i}}[(P^{\prime}-P)\cdot x^{-}-(\bm{P}_{\perp}^{\prime}-\bm{P}_{\perp})\cdot\bm{x}_{\perp}]}\\ &\times\big{\langle}{P^{\prime}}\big{|}\bar{\psi}_{i}(0,y^{-},\bm{0})\gamma^{+}\psi_{i}(0,0,\bm{0})\big{|}{P}\big{\rangle}\;.\end{split}$$

$$\big{\langle}{P^{\prime}}\big{|}{P}\big{\rangle}=(2\pi)^{3}2P^{+}\delta(P^{\prime+}-P^{+})\,\delta(\bm{P}^{\prime}_{\perp}-\bm{P}_{\perp})\;.$$

$$f_{i/\mathrm{p}}^{(0)}(\xi )$$

$$\displaystyle\quad=\frac{1}{4\pi}\int\!dy^{-}\,e^{-{\mathrm{i}}\xi P^{+}y^{-}}\big{\langle}{P}\big{|}\bar{\psi}_{i}(0,y^{-},\bm{0})\gamma^{+}\psi_{i}(0)\big{|}{P}\big{\rangle}\;.$$

$$W(y^{-},0)=\mathcal{P}\exp \left(-\mathrm{i}g{\int }_0^{y^{-}}d{\bar{y}}^{-}{A}^{+}(0,{\bar{y}}^{-},0{)}_a{t}_a\right).$$

$$\begin{split}f_{i/\mathrm{p}}^{(0)}(\xi)={}&\frac{1}{4\pi}\int\!dy^{-}\,e^{-{\mathrm{i}}\xi P^{+}y^{-}}\\ &\times\big{\langle}{P}\big{|}\bar{\psi}_{i}(0,y^{-},\bm{0})\gamma^{+}W(y^{-},0)\,\psi_{i}(0)\big{|}{P}\big{\rangle}\;.\end{split}$$

$$\begin{array}{rl}{\psi }_i(0)\to & U(0){\psi }_i(0),\\ {\bar{\psi }}_i(0,y^{-},0)\to & {\bar{\psi }}_i(0,y^{-},0)U(0,y^{-},0{)}^{-1},\\ W(y-,0)\to & U(0,y^{-},0)W(y-,0)U(0{)}^{-1}.\end{array}$$

$$f_{i/\mathrm{p}}(\xi ,{\mu }^2)=$$

$$\displaystyle\times\big{\langle}{P}\big{|}\bar{\psi}_{i}(0,y^{-},\bm{0})\gamma^{+}W(y^{-},0)\,\psi_{i}(0)\big{|}{P}\big{\rangle}\;,$$

$$f_{\bar{i}/\mathrm{p}}(\xi ,{\mu }^2)=$$

$$\displaystyle\times\big{\langle}{P}\big{|}{\rm Tr}\big{[}\gamma^{+}\psi_{i}(0,y^{-}\!,\bm{0})W(y^{-},0)\bar{\psi}_{i}(0)\big{]}\big{|}{P}\big{\rangle},$$

$$f_{g/\mathrm{p}}(\xi ,{\mu }^2)=$$

$$\displaystyle\times\big{\langle}{P}\big{|}G(0,y^{-},\bm{0})^{+\nu}W(y^{-},0)\,G(0)_{\nu}^{\ +}\big{|}{P}\big{\rangle}\;,$$

$$G_{\mu \nu }^a={\partial }_{\mu }{A}_{\nu }^a-{\partial }_{\nu }{A}_{\mu }^a-g{f}^{abc}{A}_{\mu }^b{A}_{\nu }^c,$$

$$\begin{array}{rl}& \frac{d}{d\log {\mu }_{\mathrm{F}}^2}{f}_{a/\mathrm{p}}(\xi ,{\mu }_{\mathrm{F}}^2)\\ & =\sum _{\hat{a}}{\int }_{\xi }^1\frac{dz}{z}{P}_{a\hat{a}}(z,{\alpha }_{\mathrm{s}}({\mu }_{\mathrm{F}}^2))f_{\hat{a}/\mathrm{p}}(\xi /z,{\mu }_{\mathrm{F}}^2).\end{array}$$

$$P_{a\hat{a}}(z,{\alpha }_{\mathrm{s}}({\mu }_{\mathrm{F}}^2))=$$

$$+\cdots .$$

$$\sigma [F]\approx \sum _{a,b}\iint d{\xi }_{a}d{\xi }_b{f}_{a/A}({\xi }_a,{\mu }^2)f_{b/B}({\xi }_b,{\mu }^2){\hat{\sigma }}_{ab}[F].$$

$$\begin{array}{rl}{P}_{\mathrm{A}}=& \left(P_{\mathrm{A}}^{+},0,0\right),\\ P_{\mathrm{B}}=& \left(0,P_{\mathrm{B}}^{-},0\right).\end{array}$$

$$p_i=(e^{y_i}\sqrt{({\mathbf{p}}_{i,\perp }^2+m_i^2)/2},e^{-y_i}\sqrt{({\mathbf{p}}_{i,\perp }^2+m_i^2)/2},{\mathbf{p}}_{i,\perp }).$$

Full text

parton distribution functions; large hadron collider

††preprint: MS-TP-19-09, SMU-HEP-19-05

1. I Overview

2. I.1 The parton model

3. I.2 Cross sections and factorization

4. I.3 Practical issues in the theory

5. I.4 The strong coupling

6. I.5 Heavy-quark masses

7. I.6 Special kinematic regions

8. I.7 Fitting

9. II Review of theory

10. II.1 Definition of parton distribution functions

11. II.1.1 Momenta

12. II.1.2 Parton distributions in canonical field theory

13. II.1.3 Gauge invariance

14. II.1.4 Renormalization

15. II.1.5 Gluons

16. II.1.6 Evolution equation

17. II.2 Infrared safety and factorization

18. II.2.1 Kinematics

19. II.2.2 Infrared safety

20. II.2.3 Factorization

21. II.3 Treatment of heavy quarks

22. II.4 Deeply Inelastic Scattering

23. III Statistical inference in fitting the parton distributions

24. III.1 Bayes’ theorem

25. III.2 Gaussian model for the distribution of data

26. III.3 Determining theory parameters using $\chi^{2}$

27. III.4 Another definition of $\chi^{2}$

28. III.5 Dependence on the theory parameters

29. III.6 Distribution of the parameters

30. III.7 Calculation of a cross section

31. III.8 Expectation value and variance of $\chi^{2}$

32. IV Tests of performance of the fit

33. IV.1 Testing with resampled data

34. IV.2 Dependence on the number of PDF parameters

35. IV.3 Dependence on the PDF functional form

36. IV.4 Closure test

37. IV.5 Test of the nuisance parameters

38. IV.6 Test of data residuals

39. IV.7 Value of $\chi^{2}$ from an individual experiment

40. IV.8 Test of $\chi^{2}$ from individual experiments

41. IV.9 Test of consistency between experiments

42. IV.10 A more conservative way to adjust the errors

43. IV.11 Summary of measures of goodness of fit

44. IV.12 Global and dynamic tolerance

45. V Parton distributions for heavy ions

46. V.1 Universality of nuclear PDFs

47. V.2 Parameterizing the $\mathtt{A}$-dependence

48. V.3 Comparisons of nuclear PDFs

49. VI Conclusions

50. A Transformation for $\chi^{2}(D,a,\lambda)$

I Overview

As of 2018, the Large Hadron Collider (LHC) has taken a large data sample of proton-proton collisions and has been using these data to precisely measure the properties of the Higgs boson and search for physics beyond the Standard Model. The ATLAS and CMS experiments took about 150 fb*-1* each at the center-of-mass energy of $\sqrt{s}=13\;{\rm TeV}$ and more than 22 fb*-1* of data at $\sqrt{s}=8\;{\rm TeV}$. In addition, the LHCb experiment has accumulated more than 9 fb*-1* of data at various energies, including data at extreme rapidities. On top of proton-proton collisions, all experiments at the LHC, and in particular the ALICE experiment, are taking data in proton-lead and lead-lead collisions.

All predictions at the LHC are crucially dependent on knowledge of the quark and gluon content of the proton. The probability distributions of the constituents of the nucleon, collectively called partons, need to be known to make predictions in the theoretical framework of perturbative quantum chromodynamics (PQCD). The concept of a parton originates in a model by Bjorken and Feynman Bjorken (1969); Feynman (1973), where it refers to quasi-independent pointlike constituents observed inside hadrons that undergo deeply inelastic scattering. As the formalism of PQCD developed, partons were shown to be excitations of elementary quantum fields of spins 1/2 and 1 – quarks, antiquarks, gluons, and even photons – in a hadron undergoing a hard scattering process.

Understanding of the structure of the proton is being updated continuously using the wealth of LHC and other world data. Knowledge of the structure of the nucleon needed for a large class of theoretical predictions in PQCD is encoded in the collinear parton distribution functions (PDFs). The PDFs have been determined from data starting in the early 1980s Duke and Owens (1984); Eichten et al. (1984); Gluck et al. (1982). They are now determined using the method of the global QCD analysis Morfin and Tung (1991); Martin et al. (1988, 1989); Harriman et al. (1990); Owens (1991) from experimental measurements at colliders such as HERA, the Tevatron, and the LHC, and in fixed-target experiments. The PDFs are provided in a number of practically useful forms by several collaborations, including ABM Alekhin et al. (2017), HERAPDF Abramowicz et al. (2015), CT Dulat et al. (2016), CTEQ-JLab Accardi et al. (2016a), MMHT Harland-Lang et al. (2015), and NNPDF Ball et al. (2017). The modern PDF parameterizations are provided with families of “error PDF sets” Giele and Keller (1998); Giele et al. (2001); Pumplin et al. (2001) that allow the user to assess the total uncertainty on the PDFs arising from a variety of experimental and theoretical errors. Methods for statistical combination of PDF ensembles from various groups exist Gao and Nadolsky (2014); Carrazza et al. (2015a, b), and comprehensive guidelines on uses of PDFs at the LHC are published by the PDF4LHC working group Butterworth et al. (2016). As of 2020, close to seven hundred PDF ensembles from various groups are distributed as numerical tables in the standard Les Houches Accord PDF (LHAPDF) format Buckley et al. (2015); Giele et al. (2002); Whalley et al. (2005); Bourilkov et al. (2006) from a public online repository.111https://lhapdf.hepforge.org.

The parton distributions in nuclei are equally analyzed by several collaborations: EPPS Eskola et al. (2017), nCTEQ Kovařík et al. (2016), DSSZ de Florian et al. (2012), HKN Hirai et al. (2007) KA Khanpour and Atashbar Tehrani (2016) and NNPDF Abdul Khalek et al. (2019a). They are provided either as parametrizations of the nuclear PDFs themselves or the nuclear correction factor applied to a predefined reference proton PDF.

Increasingly precise requirements will be imposed on the determination of PDFs and their uncertainties during the high-luminosity (HL) phase of the LHC operation to precisely measure Higgs boson couplings and electroweak parameters de Florian et al. (2016), and to maximize the HL-LHC reach in a variety of tests of the Standard Model and new physics searches ATLAS and CMS Collaborations (2019).

The purpose of this article is to review select topics related to the theoretical definition, determination, and usage of PDFs in modern applications. We will concentrate on methodological aspects of the PDF analysis that will be of growing importance in the near-future LHC era. We primarily focus on theoretical and statistical aspects of the determination of PDFs in the nucleon and nuclei, notably, on proper theoretical definitions, statistical inference of the PDF parameterizations from the experimental data, and factorization for heavy nuclei. This work supplements the recent reviews of phenomenological applications of PDFs available in Forte and Watt (2013); Gao et al. (2018), as well as extensive comparisons Alekhin et al. (2011); Watt and Thorne (2012); Butterworth et al. (2016); Accardi et al. (2016b) of PDFs from various collaborations and QCD predictions based on these PDFs. Introductory texts on the fundamentals of QCD factorization, global PDF analysis, and collider applications of PDFs are available, e.g., in Collins (2013); Brock et al. (1995); Campbell et al. (2017).

We begin with the parton model and its relation to QCD, the field theory of the strong interactions.

I.1 The parton model

The principal aim of contemporary particle physics is to test the current theory, the Standard Model, and to look for evidence of new physics that is not included in the Standard Model. The Standard Model is a renormalizable quantum field theory with fields for leptons ($\mathrm{e},\mu,\tau$), their associated neutrinos, quarks (d,u,s,c,b,t), vector bosons ($\gamma,\mathrm{W}^{\pm},\mathrm{Z},\mathrm{g}$), and a scalar boson field, the Higgs field. The part of the theory involving the photon, $\mathrm{W}^{\pm},\mathrm{Z}$ and Higgs fields involves spontaneous symmetry breaking and is quite subtle. The part of the theory involving the gluon, g, constitutes the theory of the strong interactions, QCD. The gluon couples to quarks, but not to leptons. We cannot offer a review of the Standard Model, but we assume that the reader has some familiarity with it. See, for instance, Peskin and Schroeder (1995); Srednicki (2007); Aitchison and Hey (2012); Schwartz (2014).

One important feature of the Standard Model, and quantum field theory in general, is that the couplings that appear at the vertices of Feynman diagrams representing the theory are best considered to be dependent on a squared momentum scale $\mu^{2}$. For instance, the QCD coupling constant, $g_{\mathrm{s}}^{2}/(4\pi)\equiv\alpha_{\mathrm{s}}$, becomes $\alpha_{\mathrm{s}}(\mu^{2})$. The dependence on $\mu^{2}$ is derived from the theory, even though the value of $\alpha_{\mathrm{s}}(\mu^{2}_{0})$ at a particular scale like $\mu^{2}_{0}=M_{\mathrm{Z}}^{2}$ is a free parameter of the theory. Which value of $\mu^{2}$ is useful in addressing a particular physical problem depends on the typical momentum scale of the problem.

The electroweak coupling, $\alpha(\mu^{2})$, increases with $\mu^{2}$, but for physically relevant momentum scales, $\alpha(\mu^{2})$ is small. Thus scattering cross sections for electroweak interactions can be usefully computed as a perturbation series in the small parameter $\alpha(\mu^{2})$.

The QCD coupling constant, $\alpha_{\mathrm{s}}(\mu^{2})$, decreases with $\mu^{2}$. It is large when $\mu$ is of order 1 GeV or less, and it reduces to $\alpha_{\mathrm{s}}(\mu^{2})\approx 0.1$ at $\mu=100\ \mathrm{GeV}$. This scale dependence of $\alpha_{\mathrm{s}}(\mu^{2})$ tells us that perturbation theory should be useful for describing the parts of a physical scattering process with hadrons that involves only very large momentum scales.

However, we run into a complication. How can we use experiments involving protons to investigate the Standard Model? For instance, we know from experiment that we can make Higgs bosons in proton-proton collisions, but how can we understand the Higgs production process quantitatively? The Higgs boson production cross section depends on the internal structure of the initial-state protons. The strong coupling $\alpha_{\mathrm{s}}(\mu^{2})$ is large for $\mu$ of order the proton mass, $m_{\mathrm{p}}$, or smaller. This means that we cannot expect perturbation theory to be useful for calculating the structure of the proton as a bound state of quarks, antiquarks, and gluons.

To understand the problem and its tentative solution, consider deeply inelastic electron scattering from a proton (DIS). The DIS process is reviewed in Sec. II.4. Suppose an electron scatters from a proton of momentum $P_{\mathrm{A}}$ by exchanging a photon with momentum $q^{\alpha}$.222We denote Lorentz indices with Latin letters and use the $\{+,-,-,-\}$ sign convention for the Minkowski metric tensor. We define $Q^{2}\equiv-q^{2}$. Since the 4-momentum $q^{\alpha}$ is spacelike, we have $Q^{2}>0$. We demand that $Q^{2}$ be much larger than $1\mbox{ GeV}^{2}$. Then there is some hope for a perturbative approach that uses an expansion in powers of $\alpha_{\mathrm{s}}(Q^{2})$. The lowest-order Feynman diagram, corresponding to $t$-channel electromagnetic scattering $e+q\!\!\!\!{}^{{}^{\small(-)}}\rightarrow e+q\!\!\!\!{}^{{}^{\small(-)}}$ of an electron on a quark or antiquark, is pretty simple. But how do we relate the initial-state quark field to the proton?

Here is an approach we can follow. Let us use a “brick-wall” reference frame in which $\vec{q}$ lies entirely in the $-z$ direction, $q^{0}=0$, and in which $\vec{P}_{\mathrm{A}}$ lies entirely in the $+z$ direction. Then $|q^{z}|=Q$, and we can consider the hard quark-photon interaction to be localized in a time interval $\Delta t_{\textsc{h}}=1/Q$. In DIS, we also demand that $P_{\mathrm{A}}\cdot q$ be large, of order $Q^{2}$. Then the proton momentum is large, with $P_{\mathrm{A}}^{z}=(P_{\mathrm{A}}\cdot q)/Q$. This means that the proton is highly boosted, with a boost factor $e^{\omega}=(2P_{\mathrm{A}}\cdot q)/(m_{\mathrm{p}}Q)$. In the proton rest frame, we can suppose that the time period between successive quark-gluon interactions is of order $\Delta t_{s}^{\mathrm{rest}}\sim 1/m_{\mathrm{p}}$. In the “brick-wall” reference frame, the typical time for these soft internal interactions is $\Delta t_{\textsc{s}}\sim e^{\omega}/m_{p}=2P_{\mathrm{A}}\cdot q/(m_{\mathrm{p}}^{2}Q)$ as a consequence of relativistic time dilation. The time interval $\Delta t_{\textsc{s}}$ is much longer than the hard interaction time scale: $\Delta t_{\textsc{s}}/\Delta t_{\textsc{h}}\sim 2P_{\mathrm{A}}\cdot q/m_{\mathrm{p}}^{2}$. This argument implies that the quark-gluon interactions inside the proton are largely frozen during the time interval $\Delta t_{\textsc{h}}$ while the hard scattering interaction with the virtual photon is taking place. The struck quark is effectively free.

This gives us the parton model of Feynman and Bjorken Bjorken (1969); Feynman (1973). The word parton as used now refers to a quark, antiquark, gluon, and sometimes to a photon. Originally, partons meant just the constituents of the proton, whatever those constituents might be.

It is easy to use the parton model to describe the cross section for DIS. We assume that, inside a fast moving proton, there are partons of various types $a$. A parton can carry a momentum $k^{\alpha}$ that is a fraction $\xi$ of the momentum of the proton: $k^{\alpha}=\xi P_{\mathrm{A}}^{\alpha}$. Denote by $f_{a/\mathrm{p}}(\xi)d\xi$ the probability to find such a parton with momentum fraction between $\xi$ and $\xi+d\xi$. Let $d\hat{\sigma}/(dQ^{2}d(P_{\mathrm{A}}\cdot q))$ be the cross section for the electron to scatter from a free quark with momentum $k=\xi P_{\mathrm{A}}$, calculated in lowest-order perturbation theory. Then the cross section for this scattering from a proton should be

[TABLE]

This picture is simple and intuitive. An analogous picture covers processes such as Higgs boson production in high energy proton-proton collisions.

Unfortunately, the parton model as just described does not survive scrutiny when we include gauge field interactions at higher orders of $\alpha_{\mathrm{s}}$: one encounters contradictions as soon as one tries to calculate $d\hat{\sigma}/[dQ^{2}\,d(P_{\mathrm{A}}\cdot q)]$ beyond the leading order in QCD. The basic problem is that QCD is a quantum field theory, in which interactions among the quarks and gluons occur at all time or distance scales, including time scales much smaller than the $\Delta t_{\textsc{s}}$ scale that one would naively associate with the interior motions inside a highly boosted proton.

Nevertheless, one can turn the parton model picture into a theoretically consistent framework that includes higher-order radiative contributions. One first needs to carefully define what one means by a parton distribution function $f_{a/\mathrm{p}}(\xi)$. With a careful definition, the parton distribution functions become $f_{a/\mathrm{p}}(\xi,\mu^{2})$, with a dependence on a scale $\mu^{2}$. Then cross sections for DIS and for many processes in hadron-hadron scattering have a property known as factorization, which means that they satisfy a formula very similar to Eq. (1).

I.2 Cross sections and factorization

In this paper, we will review how the PDFs are systematically defined in the QCD theory and determined by applying statistical inference to experimental observations sensitive to the PDFs. We will also review the general formalism to estimate the uncertainty on the PDFs that results from fitting the experimental data.

It is most common to define and determine the the PDFs in the nucleon in the $\overline{\mathrm{MS}}$ factorization scheme discussed in some detail in Sec. II.1. Intuitively, these functions, $f_{a/A}(\xi,\mu^{2})$, represent the probability to find a parton of type $a$ (a gluon or a particular flavor of quark or antiquark) in a hadron of type $A$, for example a proton, as a function of the fraction $\xi$ of the momentum of the hadron that is carried by the parton. The argument $\mu^{2}$ in $f$ indicates the momentum scale at which the parton distribution function applies. The $\mu^{2}$ dependence is given by the Dokshitser-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations Dokshitzer (1977); Gribov and Lipatov (1972); Altarelli and Parisi (1977) that we describe in Sec. II.1.6. With the aid of these evolution equations, the functions $f_{a/A}(\xi,\mu^{2})$ can be determined from the functions $f_{a/A}(\xi,\mu_{0}^{2})$ at a scale $\mu_{0}^{2}$ that is typically chosen to be around $1\ \mathrm{GeV}^{2}$. The functions at scale $\mu_{0}^{2}$ cannot be calculated in perturbation theory.

The PDFs at the starting scale $\mu_{0}^{2}$ are determined from experimental data. Consider first a cross section $\sigma[F]$, defined by integrating the completely differential cross section for any number of final-state particles, multiplied by functions $F$ that describe what is measured in the final state. For instance, if we observe a single weakly interacting particle, the measurement function $F$ might be simply a product of delta functions that specify the energy and direction of momentum of the particle. We thus integrate over the momenta of particles that are not measured, giving us an inclusive cross section. The observable $\sigma[F]$ must be “infrared-safe”, as described later in Sec. II.2.2. For lepton-hadron scattering, the cross section $\sigma[F]$ is related to parton distributions by

[TABLE]

For cross sections at hadron colliders, a parton distribution function is needed for each of two colliding hadrons:

[TABLE]

We will review this formula in more detail in Sec. II.2.3.

In order to determine the PDFs at the starting scale $\mu_{0}^{2}$, one selects observables that are sensitive to different combinations of parton distributions. The parton distributions at scale $\mu_{0}^{2}$ are parameterized by a sufficiently flexible functional form. The observables are first calculated using the parton distributions, where the free parameters are given some initial values and are compared to data. The parameters are then adjusted until the theoretical predictions describe the data well.

I.3 Practical issues in the theory

Following this very simple strategy requires in reality a very detailed understanding of many facets of perturbative QCD. First, as the precision of the determination of parton distributions needs to match the experimental precision, fitted observables are calculated at the next-to-next-to-leading order (NNLO) in perturbative QCD for nucleon PDFs. The NNLO accuracy in the global fits corresponds to computing the hard cross sections by including perturbative radiative contributions suppressed by up to two powers of $\alpha_{s}$. Such predictions at NNLO for hard processes suitable for determination of PDFs are increasingly available Moch et al. (2005); Vermaseren et al. (2005); Zijlstra and van Neerven (1991); van Neerven and Zijlstra (1991); Buza and van Neerven (1997); Catani and Grazzini (2007); Catani et al. (2009); Gavin et al. (2011, 2013); Li and Petriello (2012); Gehrmann-De Ridder et al. (2016, 2018); Currie et al. (2017a, b); Campbell and Ellis (2010); Boughezal et al. (2017); Berger et al. (2016). There are also partial results at $\mathrm{N}^{3}\mathrm{L}\mathrm{O}$, such as Vermaseren et al. (2005); Moch and Rogal (2007); Moch et al. (2009); Bierenbaum et al. (2009); Ablinger et al. (2011, 2014). Theoretical predictions used for nuclear PDFs are still typically at next-to-leading-order (NLO), even though first NNLO nuclear PDF analyses exist.

Incorporating the theoretical calculations into the fit requires a careful selection of observables that are theoretically well-defined (infrared-safe) and can be calculated up to the required order. Due to the nature of higher-order calculations, the numerical evaluation can be time-consuming. This shortcoming is usually solved either by using precomputed tables of computationally slow point-by-point NNLO corrections applied to fast NLO calculations, or increasingly by using fast gridding techniques, such as the ones implemented in the fastNLO Wobisch et al. (2011), APPLGRID Carli et al. (2010), aMCFast Bertone et al. (2014b), and NNPDF FastKernel Forte et al. (2010) programs. The comparably accurate and fast DGLAP evolution of PDFs up to NNLO accuracy is implemented in a number of public codes: PEGASUS Vogt (2005), HOPPET Salam and Rojo (2009), QCDNUM Botje (2011), or APFEL Bertone et al. (2014a).

Even after implementing the measured observables and the corresponding DGLAP evolution for the PDFs at (N)NLO, one still has to address a number of issues that become important as the PDF analysis is pushed towards higher precision. On the experimental side, the NNLO PDFs are increasingly constrained by high-luminosity measurements, in which the statistical experimental errors are small, and adequate implementation of many (sometimes hundreds) of correlated systematic uncertainties is necessary. Commonly followed procedures for implementation of systematic uncertainties in the PDF fits are reviewed in Appendix A of Ref. Ball et al. (2013b). We discuss the treatment of systematic uncertainties in some detail in Sec. III.

From the side of theory, subtle radiative contributions, such as NLO electroweak or higher-twist contributions, are comparable to NNLO QCD contributions in some fitted observables. The photon constituents contribute at a fraction-of-percent level to the total momentum of the proton. The associated parton distribution for the photon can be computed very accurately using the structure functions and nucleon form factors from lepton-hadron (in)elastic scattering as the input Manohar et al. (2016, 2017). The resulting LUXqed parameterization of the photon PDF is already implemented in Bertone et al. (2017b); Nathvani et al. (2018). An alternative is to fit a phenomenological parameterization of the photon PDF at the initial scale of evolution together with the rest of the PDFs. Such phenomenological parameterizations constrained just by the global fit Ball et al. (2013a); Schmidt et al. (2016); Giuli et al. (2017) are less precise than the LUXqed form.

Even at NNLO, the residual theoretical uncertainties due to missing higher-order contributions in $\alpha_{s}$ may have an impact on PDFs. Such theory uncertainties are partly correlated in a generally unknown way across experimental data points. In addition to the traditional estimation of higher-order contributions by the variation of factorization and renormalization scales, recently, more elaborate methods for estimation of higher-order uncertainties have been explored with an eye on applications in the PDF fits, such as Olness and Soper (2010); Cacciari and Houdeau (2011); Gao (2011); Forte et al. (2014); Harland-Lang and Thorne (2019); Abdul Khalek et al. (2019c, b). See the discussion at Eq. (79).

I.4 The strong coupling

Another associated issue is the treatment of the strong coupling $\alpha_{s}(\mu^{2}_{R})$. The strong coupling depends on a scale $\mu^{2}_{R}$, called the renormalization scale. For a review, see any text on QCD, for example Collins (2013). Since $\alpha_{s}(\mu^{2}_{R})$ obeys a renormalization group equation, its value at any $\mu^{2}_{R}$ can be determined from its value at a fixed scale $\mu_{R0}^{2}$. Normally, one sets $\mu_{R0}$ to be the mass $M_{\mathrm{Z}}$ of the $Z$-boson. All QCD observables depend on $\alpha_{s}(M_{\mathrm{Z}}^{2})$, and so do the fitted PDFs. Conventionally, the world-average value of $\alpha_{s}(M_{\mathrm{Z}}^{2})$ Tanabashi et al. (2018) is derived from a combination of experimental measurements, with the tightest constraints imposed by the QCD observables that do not depend on the PDFs, notably, hadroproduction in electron-positron collisions, hadronic $\tau$-decays, and quarkonia masses.

Other useful constraints on $\alpha_{s}$ are imposed by a variety of hadron-scattering observables (lepton-hadron DIS, jet and $t\bar{t}$ production, …) that are simultaneously sensitive to PDFs, but the constraints of this class are generally weaker and more susceptible to systematic effects. As some hadronic observables of the latter class are also included in the global fit to constrain the PDFs, in principle these observables can determine both $\alpha_{s}$ and the PDFs at the same time. Consequently several treatments of $\alpha_{s}$ exist in the current PDF analyses. Most PDF groups publish some global fits that determine $\alpha_{s}(M_{\mathrm{Z}}^{2})$ and PDFs simultaneously. They typically find that the best-fit $\alpha_{s}(M_{\mathrm{Z}}^{2})$ is consistent with in the world average of $\alpha_{s}(M_{\mathrm{Z}}^{2})$, but with a considerably larger uncertainty than in the world average. ABM fits are representative of this approach Alekhin et al. (2017), as well as the dedicated studies performed in the global framework Ball et al. (2018b); Thorne et al. (2018).

As an example, the best-fit $\alpha_{s}(M_{\mathrm{Z}}^{2})=0.1166\pm 0.0027$ in CT18 at NNLO Hou et al. (2019) is consistent with the world average $\alpha_{s}(M_{\mathrm{Z}}^{2})=0.1181\pm 0.0011$ Tanabashi et al. (2018), as well as with $\alpha_{s}(M_{\mathrm{Z}}^{2})\approx 0.1176$ in Thorne et al. (2018) and a somewhat higher $\alpha_{s}(M_{\mathrm{Z}}^{2})=0.1185\pm 0.0012$ in Ball et al. (2018b). The CT18 value results as some trade-off between the DIS experiments (notably, fixed-target DIS), which collectively prefer a somewhat lower $\alpha_{s}(M_{\mathrm{Z}}^{2})\approx 0.115$, and jet+top and Drell-Yan experiments, which collectively prefer a higher $\alpha_{s}(M_{\mathrm{Z}}^{2})\approx 0.119$ Hou et al. (2019). The quoted uncertainty on $\alpha_{s}(M_{\mathrm{Z}}^{2})$ at 68% probability level varies among the PDF-fitting groups in 2020 from about 0.001 to 0.0025 depending on the adopted definitions for the uncertainty (with the CT18 uncertainty being the most conservative).

On the other hand, it is often advantageous to perform the PDF fits and determine the PDF uncertainty at a fixed world-average value of $\alpha_{s}(M_{\mathrm{Z}}^{2})$, then estimate the $\alpha_{s}$ uncertainty of the fits by using a few PDF fits with alternative $\alpha_{s}$ values. If the uncertainties obey a Gaussian probability distribution, it can be rigorously demonstrated that, to compute the total PDF+$\alpha_{s}$ uncertainty that includes all correlations, it suffices to add the resulting PDF and $\alpha_{s}$ uncertainties in quadrature Lai et al. (2010b). The empirical probability distributions in the PDF fits are indeed sufficiently close to being Gaussian, so this prescription for computing the PDF+$\alpha_{s}$ uncertainty is adopted by the majority of recent fits. For example, the PDF4LHC group recommends Butterworth et al. (2016) calculating the PDF+$\alpha_{s}$ uncertainty at the 68% confidence level by adding in quadrature the PDF uncertainty computed using 30 (100) PDF4LHC15 error sets for the world-average $\alpha_{s}(M_{\mathrm{Z}}^{2})=0.1180$, and the $\alpha_{s}$ uncertainty computed from two best-fit PDF sets for $\alpha_{s}(M_{\mathrm{Z}}^{2})=0.1165$ and $0.1195$.

I.5 Heavy-quark masses

Another issue that needs to be addressed in global fits of PDFs is the treatment of massive charm and bottom quarks. Mass effects play an important role in describing, for example, subprocesses with charm (anti)quarks in DIS. There are several approaches to treating the mass of the quark, such as the zero-mass variable-flavour-number (ZM-VFN) scheme Collins et al. (1978); Collins and Tung (1986) or the general-mass variable-flavor-number (GM-VFN) scheme Collins (1998); Aivazis et al. (1994); Krämer et al. (2000); Thorne and Roberts (1998); Thorne (2006); Forte et al. (2010). In Sec. II.3, we provide a short pedagogical introduction to the otherwise extensive topic of the treatment of masses for heavy quarks. For a more thorough review of the heavy-quark schemes we refer the reader to Butterworth et al. (2016); Accardi et al. (2016b).

I.6 Special kinematic regions

Some kinematic regions require special treatment if their respective experimental data are to be included in a global fit of PDFs.

One such region is where the typical momentum fraction $\xi$ is large, but the momentum scale of the process is not too large. Then the nonzero mass of the proton (or in general the target), normally neglected, may need to be taken into account. These target-mass corrections are discussed in detail in Schienbein et al. (2008). In the same region, deuterium nuclear corrections also have to be considered, when data in this specific kinematic region are taken on deuterium rather than proton targets Accardi et al. (2010, 2011).

The other kinematical region in need of careful treatment is the one where the momentum transfer $Q^{2}$ is low, typically below 4 $\mbox{GeV}^{2}$, at moderate or large $\xi$. Here power corrections can become important and can be taken into account, for example as in Martin et al. (2004); Alekhin et al. (2012); Thorne (2014).

Yet another such region arises in DIS at small $x$ and $Q$, roughly satisfying $Q^{2}<A_{cut}/x^{\lambda}$ with $A_{cut}\sim 0.5-1.5\mbox{ GeV}^{2}$, and $\lambda\sim 0.3$ Golec-Biernat and Wusthoff (1999); Caola et al. (2010). This is the limit where summation of small-$x$ logarithms will become necessary, and indeed, a slowdown in perturbative convergence of inclusive DIS cross sections and resulting small-$x$ PDFs is observed even at NNLO in the affected HERA region Abramowicz et al. (2015). The DIS data in this region provide valuable constraints on the small-$x$ behavior of the gluon PDF. The small-$x$ instability in DIS can be cured to a certain extent by inclusion of power-suppressed (higher-twist) contributions Harland-Lang et al. (2016) or, quite effectively in the HERA region at $Q^{2}>4\mbox{ GeV}^{2}$, by using an $x$-dependent factorization scale $\mu^{2}$ in NNLO DIS cross sections in some of the PDF sets (CT18X and CT18Z) published in Hou et al. (2019). Summation of small-$x$ logarithms, matched to NNLO, has been successfully implemented in NNPDF Ball et al. (2018a) and xFitter analyses Abdolmaleki et al. (2018). It results in even better description of the accessible small-$x$ region. Either NNLO+NLLx summation, as in Ball et al. (2018a), or the choice of a special $x$-dependent factorization scale in the fixed-order NNLO cross section, as in Hou et al. (2019), thus leads to a better description of the small-$x$ subsample of the HERA DIS data. The resulting PDFs obtained after these changes tend to have elevated gluon and strangeness components at $x<10^{-2}$ as a result of slower $Q$ dependence of the DIS cross sections at small $x$ than would be predicted at a fixed order with a standard scale $\mu=Q$. LHC predictions based on such modified PDFs, such as CT18Z, may lie outside of the nominal error bands of the PDF set with default choices, such as CT18.

I.7 Fitting

After addressing all necessary features of theory predictions such as the ones spelled out in the previous paragraphs, one compares the theory predictions to the experimental data. The process of fitting the theoretical predictions to data by adjusting the PDFs is the main focus of this review. The reason is that proper determination of PDF uncertainties will be highly important for the analysis of the high-luminosity LHC data, as the PDF uncertainty will soon dominate systematic uncertainties on the theory side in key tests of electroweak symmetry breaking, including the measurements of Higgs couplings and mass of the charged weak boson de Florian et al. (2016); ATLAS and CMS Collaborations (2019). The statistical framework of the PDF fits is fundamentally more complex than the one in the electroweak precision fits: while the parametric model of the electroweak fits is uniquely determined by the Standard Model Lagrangian, the parametric model for the parton distribution functions may change within some limits in order to optimize agreement between QCD theory and data.

Consequently, the PDF uncertainty is comprised of four categories of contributions:

Experimental uncertainties, including statistical and correlated and uncorrelated systematic uncertainties of each experimental data set; 2. 2.

Theoretical uncertainties, including the absent higher-order and power-suppressed radiative contributions, as well as uncertainties in using parton showering programs to correct the data in order to compare to fixed-order perturbative cross sections; 3. 3.

Parameterization uncertainties associated with the choice of the PDF functional form; 4. 4.

Methodological uncertainties, such as those associated with the selection of experimental data sets, fitting procedures, and goodness-of-fit criteria.

As an illustration, the left panel of Fig. 1 shows the HERAPDF2.0 parameterizations determined from the fits exclusively to DIS data. The PDF uncertainty corresponding to the PDF solutions covering 68% of the cumulative probability is comprised of the experimental, theoretical model, and parameterization components that were estimated for a select fitting methodology Abramowicz et al. (2015). Other groups may not separate all four components listed above in the total PDF uncertainty. In the right panel of Fig. 1, the CT18 NNLO PDF uncertainty bands are evaluated for 68% cumulative probability according to a two-tier goodness-of-fit criterion Lai et al. (2010a) that accounts both for the agreement with the totality of fitted data and with individual experimental data sets. The CT18 analysis includes a variety of data sets on DIS, vector boson, jet, and $t\bar{t}$ production. While this diversity of data allows one to resolve differences between PDFs of various flavors and probe a broader range of PDF parameterization forms, in practice some incompatibilities (“tensions”) between constraints on the PDFs from various experiments are introduced and need to be either eliminated or accounted for in the PDF uncertainty estimate. [The CT18 and HERAPDF2.0 PDFs are fitted to 3690 and 1130 data points, respectively. About 100 different parameterization forms have been tried in the CT18 analysis, contributing to the spread of the PDF uncertainty.] The width of the CT18 error bands thus depends on a two-level tolerance convention Lai et al. (2010a); Pumplin et al. (2002) that is adjusted so as to reflect PDF variations associated with some disagreements between experiments, parameterization and theoretical uncertainties. We notice, for example, that the CT18 error bands for some poorly constrained flavors, notably, the strangeness PDF $xs(x,\mu^{2})$ at small $x$ (green band), may be broader than the respective HERAPDF2.0 error bands at the same probability level, despite having more experimental data included in the CT18 analysis compared to HERAPDF. The wider error bands reflect, for a large part, the spread in the acceptable PDFs estimated using the CT18 flexible parameterization forms, but also some inflation of the experimental uncertainty to reflect the imperfect agreement among experiments. Sec. IV.8 shows how to examine several experiments for their agreement.

To find most likely solutions for PDFs and establish the respective uncertainties, one must answer a fundamental question: how good, actually, is each PDF fit? We explore this question and advocate using a strong set of goodness-of-fit criteria that go beyond the weak criterion based on just the value of the goodness-of-fit function $\chi^{2}$.

All PDF fitters employ some version of minimization of the goodness-of-fit function

[TABLE]

where $D_{i}$ are the data values, $T_{i}(a)$ the corresponding theory predictions, which depend on free PDF parameters $a$, and $C_{ij}$ is the covariance matrix. The goodness-of-fit function is used to assess the quality of the theoretical description of the data and to estimate the uncertainty in the determination of the fit parameters $a$. The statistical foundations that motivate the use of the goodness-of-fit functions are discussed in Sec. III. While it is most common to find the global minimum of $\chi^{2}(a)$ numerically, many insights about the global fits can be gleaned analytically in the so called Hessian approach to fitting the PDFs, first developed in Stump et al. (2001); Pumplin et al. (2001) and refined ever since. We will discuss this approach in a great detail in order to illustrate various aspects of the fits. It relies on the observation that the PDFs approximately obey the multivariate Gaussian probability distribution in the well-constrained kinematic regions, which in turn allows one to derive the key outcomes of the PDF analysis in a closed algebraic form. For example, the Hessian method is commonly used at the end of the fit to quantify the uncertainty on the resulting PDFs. There is a powerful alternative approach, discussed only briefly in this review, to find the best-fit PDFs and determine their uncertainties using stochastic (Monte-Carlo) sampling of PDFs Giele and Keller (1998); Giele et al. (2001) and PDF parameterizations by neural networks Forte et al. (2002). Although our results are demonstrated in the Hessian approximation, they also elucidate the numerical outcomes of the Monte-Carlo sampling PDF analyses such as the one used by NNPDF. They also apply to the approximate techniques for updating the published PDF ensembles with information from new data by statistical reweighting of PDF replicas in the Monte-Carlo Giele and Keller (1998); Ball et al. (2011, 2012); Sato et al. (2014) or Hessian Paukkunen and Zurita (2014); Schmidt et al. (2018); Watt and Thorne (2012) representations.

Sec. IV is devoted to the discussion of tests one can perform to determine the extent to which the fitting procedure is consistent with the statistical hypotheses used in the procedure. This leads to a discussion of the strong goodness-of-fit set of criteria.

II Review of theory

In this section, we provide a brief overview of the theory of PDFs and their relation to cross sections. We start with the definition of PDFs as matrix elements of quantum field operators. Then we discuss the factorization property of QCD, which allows us to relate certain kinds of cross sections to PDFs and perturbatively calculated quantities. Finally, we turn to the treatment of heavy quarks in these relations, although we treat this complex subject only briefly.

II.1 Definition of parton distribution functions

In this subsection, we give definitions for PDFs as matrix elements in a proton (or other hadron) of certain operators. Instead of simply stating the definitions, we motivate them from basic field theory, following the reasoning in Collins and Soper (1982). For more details, one can consult the book Collins (2013).333Our conventions follow the Particle Data Group Tanabashi et al. (2018) and the book Collins (2013). In particular, we choose the sign of the strong coupling $g$ so that the quark-gluon vertex is $-{\mathrm{i}}g\gamma^{\mu}t^{a}$. This is the opposite from the choice in Collins and Soper (1982).

II.1.1 Momenta

Consider a proton with momentum $P$ along the +$z$ direction. We define $+$ and $-$ components of vectors using $v^{\pm}=(v^{0}\pm v^{3})/\sqrt{2}$. Then $P$ has components

[TABLE]

It is helpful to think of $P^{+}$ as being very large, $9.2\ \mathrm{TeV}$ for the LHC, but the size of $P^{+}$ does not matter for the definition of PDFs.

In the exposition below, we frequently have a function of a four-vector $v$. We use two alternative notations: either $f(v)$ or $f(v^{+},v^{-},\bm{v}^{\perp})$.

We seek to define PDFs, $f_{a/\mathrm{p}}(\xi,\mu^{2})$, which can be interpreted as giving the probability density for finding a parton of flavor $a$ (a quark, antiquark, or gluon), which carries a fraction $\xi$ of $P^{+}$, in the proton with momentum $P$. This function depends on a momentum squared scale $\mu^{2}$ at which one imagines measuring the presence of the parton.

II.1.2 Parton distributions in canonical field theory

To get started with the definition of PDFs, consider an operator $b(\xi P^{+},\bm{k}_{\perp},s,c;i)$ that destroys a quark of flavor $i$ having helicity $s$, color $c$, +-momentum $\xi P^{+}$ and transverse momentum $\bm{k}_{\perp}$. This quark then carries a fraction $\xi$ of the $+$-momentum $P^{+}$ of the proton. The adjoint operator $b^{\dagger}(\xi P^{+},\bm{k}_{\perp},s,c;i)$ then creates a quark with the same quantum numbers. We normalize the creation and destruction operators to have anticommutation relations

[TABLE]

Additionally, we suppose that the vacuum state \big{|}{0}\big{\rangle} has no quarks in it, so

[TABLE]

With the quark creation and destruction operator at hand, we can construct the operator that counts the number of quarks in a region of $\xi$ and $\bm{k}_{\perp}$:

[TABLE]

The reader can verify that, if \big{|}{\Psi}\big{\rangle} is obtained by applying quark creation operators to the vacuum, then the integral of $\rho$ over a momentum-space volume ${\cal V}_{3}$ counts the number $N({\cal V}_{3})$ of quarks in ${\cal V}_{3}$:

[TABLE]

We want to define a parton density, the number of partons $f_{i/\mathrm{p}}(\xi)$ of flavor $i$ per unit $d\xi$ in a proton. We can take a matrix element of $\rho$ in a proton state to define this:

[TABLE]

Here, for simplicity, we consider the proton to be spinless, but one can substitute a spin average: \frac{1}{2}\sum_{s_{\mathrm{p}}}\big{\langle}{P^{\prime},s_{\mathrm{p}}}\big{|}\cdots\big{|}{P,s_{\mathrm{p}}}\big{\rangle}. As noted at the beginning of this section, we take the proton momentum $P$ to be along the $z$-axis, so that $\bm{P}_{\perp}=0$. However $P^{+}$ is arbitrary.

We have given $f$ a superscript $(0)$ to indicate that this is a preliminary version of the needed definition.

To make this definition more useful, we can relate the quark creation and destruction operators to the quark field operator $\psi_{i}(x)$. For this purpose, we use the version of QCD quantized on planes of equal $x^{+}=(x^{0}+x^{3})/\sqrt{2}$ instead of planes of equal time $t=x^{0}$ Kogut and Soper (1970); Bjorken et al. (1971). The fields obey canonical commutation relations on planes of equal $x^{+}$. To make this work, we use the gauge $A^{+}(x)=0$ for the gluon field. With this way of writing the theory, the two components of the four-component Dirac field projected by $P_{\rm dy}=\frac{1}{2}\,\gamma^{-}\gamma^{+}$ (such that $P_{\rm dy}^{2}=P_{\rm dy}$) are the independent dynamical fields (‘dy’) representing quarks. The dynamical part of the quark field at $x^{+}=0$ is related to quark and antiquark creation and destruction operators by

[TABLE]

Here $k\cdot x=k^{+}x^{-}-\bm{k}_{\perp}\cdot\bm{x}_{\perp}$, and $d^{\dagger}$ is an antiquark creation operator, analogous to the quark creation operator $b^{\dagger}$. The field $\psi$ carries a flavor index $i$. It also carries a color index $c$, which we normally suppress. The spinors $u$ and $v$ are the usual solutions of the free Dirac equation, normalized to $\overline{u}(k,s)\gamma^{+}u(k,s)=2k^{+}$ and $\overline{v}(k,s)\gamma^{+}v(k,s)=2k^{+}$. Then one easily finds that $P_{\rm dy}\,u(k,s)$ and $P_{\rm dy}\,v(k,s)$ depend only on the $+$ component of $k$.

When we combine Eqs. (10) and (11), we obtain, quite directly,

[TABLE]

We can eliminate the factor \big{\langle}{P^{\prime}}\big{|}{P}\big{\rangle} by using translation invariance to write

[TABLE]

Then we can perform the $x^{-}$ and $\bm{x}_{\perp}$ integrations to give delta functions that set $P^{\prime}$ to $P$. We normalize our proton state vectors to

[TABLE]

Then the delta functions from \big{\langle}{P^{\prime}}\big{|}{P}\big{\rangle} in Eq. (12) cancel. We set $P^{\prime}$ to $P$ to get

[TABLE]

We have presented this result in some detail to emphasize that the PDF for quarks is simply the proton matrix element of the number density operator for quarks as obtained in canonically quantized field theory.

II.1.3 Gauge invariance

Next, without changing $f_{i/\mathrm{p}}^{(0)}(\xi)$, we can rewrite the definition in a way that makes it gauge-invariant. The canonical field theory that our derivation has relied on makes use of the lightlike axial gauge $A^{+}(x)=0$ for the gluon field. In an arbitrary gauge, we merely insert a Wilson line factor,

[TABLE]

This is a matrix in the color indices carried by the quark fields; $t_{a}$ is the $SU(3)_{c}$ generator matrix in the $\bm{3}$ representation. The $\cal P$ indicates path ordering of the operators and matrices, with more positive $y^{-}$ values to the left. The revised definition is

[TABLE]

The factor $W$ is just 1 if we use $A^{+}(x)=0$ gauge. If we change the gauge by a unitary transformation $U(x)$, we replace

[TABLE]

Thus, when we include the operator $W$, the right-hand side of the equation is invariant under a change of gauge.

If we use a covariant (Bethe-Salpeter) wave function for the proton state, we can use Eq. (17) for perturbative calculations. The field $\psi_{i}(0)$ absorbs a quark line from the wave function. Similarly, $\bar{\psi}_{i}(0,y^{-},\bm{0})$ creates a quark line that goes into the conjugate wave function. These quark lines can emit and absorb gluons. The factor $W(y^{-},0)$ is conveniently written as $W(y^{-},\infty)$ times $W(\infty,0)$. The operators $W$ contains gluon fields that create and absorb gluons. In a simple intuitive picture, we don’t just destroy a quark at position 0, leaving its color with nowhere to go. Rather we scatter it, so that it moves to infinity along a fixed lightlike line in the minus direction, carrying its color with it. Then its color comes back to $(0,y^{-},\bm{0})$ to provide the color for the quark that we create.

II.1.4 Renormalization

The function $f_{i/\mathrm{p}}^{(0)}(\xi)$ has so far been defined using “bare” fields, a bare coupling, bare parton masses, and a bare operator product of fields in a canonical formulation of the field theory. This will not do. Even the simplest one-loop calculation reveals that the bare $\alpha_{s}$, quark masses, and $f_{i/\mathrm{p}}^{(0)}(\xi)$ contain ultraviolet (UV) divergences. Thus we need to renormalize everything. The standard way to do this is to apply $\overline{\text{MS}}$ renormalization with scale $\mu^{2}$. For this, we need to choose a number $N_{f}$ of active flavors.444 For instance, if we are following the $N_{f}=5$ convention, then neither $\alpha_{s}$ nor PDFs include contributions from top quarks. Then top-quark virtual loops can still occur within the Feynman diagrams, but they are treated using the CWZ prescription Collins et al. (1978), in which the UV divergencies that they introduce are subtracted at zero incoming momenta and do not affect scale dependence of $\alpha_{s}$ or PDFs. Depending on the circumstances, one uses different numbers of active flavors, as we will discuss in Sec. II.3. The $\overline{\text{MS}}$-renormalized entities acquire dependence on $\mu^{2}$, which can be chosen so as to improve perturbative convergence for the short-distance cross section $\widehat{\sigma}$ in Eq. (3). Physically, quark and gluon interactions at distance scales smaller than $1/\mu^{2}$ are not resolved in these objects.

This gives us our final definition for the quark distribution Collins and Soper (1982),

[TABLE]

where $W(y^{-},0)$ is given by Eq. (16). We understand now that the formulas refer to fields and couplings and field products that are renormalized with the $\overline{\text{MS}}$ prescription for all active quarks and gluons.

For antiquarks, the analogous definition is

[TABLE]

where the color generator matrices in $W(y^{-},0)$ are in the $\bar{\bm{3}}$ representation of SU(3).

We should understand that no approximations are made in Eq. (19). In particular, we do not treat the quarks as being massless. We cannot calculate $f_{i/\mathrm{p}}(\xi,\mu^{2})$ at any finite order of perturbation theory, but we could, in principle, calculate it using lattice gauge theory. In such a calculation, we would use our best estimates for the parameters in the QCD lagrangian, including the strong coupling and the quark masses.

In fact, PDFs can be calculated using lattice gauge theory Lin et al. (2018), but the accuracy of such calculations is still limited. One can obtain much better accuracy by fitting the parton distributions to data, as described in this review. However, the definition of the parton distributions is not affected by the approximations that we make in the fitting procedure. For instance the calculated cross sections used in the fit can be leading order (LO), next-to-leading order (NLO), or next-to-next-to-leading order (NNLO). The resulting fits are often referred to as LO, NLO, or NNLO. However, it is the fits that carry these designations. The functions $f_{i/\mathrm{p}}(\xi,\mu^{2})$ that we are trying to estimate are non-perturbative objects whose definitions are independent of the fitting method.

II.1.5 Gluons

In the previous subsections, we have defined the PDFs for quarks and antiquarks by beginning with the number operator for quarks or antiquarks in unrenormalized canonical field theory using null-plane quantization in $A^{+}=0$ gauge. The starting definition is then generalized to be gauge invariant and use $\overline{\text{MS}}$ renormalized operators. One can follow the same sort of logic for the gluon field. We simply state the result Collins and Soper (1982):

[TABLE]

where $G_{\mu\nu}$ is the gluon field operator,

[TABLE]

and the Wilson line operator $W(y^{-},0)$ is given by Eq. (16), now using SU(3) generator matrices $(t_{a})_{bc}=-i\ f_{abc}$ in the adjoint representation.

II.1.6 Evolution equation

Let us take a closer look at the renormalization of the PDFs. The $\overline{\text{MS}}$ renormalization of the strong coupling $\alpha_{s}\equiv g^{2}/(4\pi)$ and the fields $\psi_{i}(x)$ and $A^{\mu}(x)_{a}$ in $n=4-2\epsilon$ dimensions proceeds in the usual way by subtracting $1/\epsilon$ poles and some finite terms from two-point subgraphs, three-point subgraphs, and four-gluon subgraphs with loops containing gluons and the $N_{f}$ active quarks. Another sort of pole arises in matrix elements for PDFs from operator products like $\bar{\psi}_{i}(0,y^{-},\bm{0})\gamma^{+}\psi_{i}(0)$ in Eq. (19). Consider a graph in which a gluon is emitted from a propagator representing the quark that is destroyed by $\psi_{i}(0)$, then absorbed by a propagator representing the quark created by $\bar{\psi}_{i}(0,y^{-},\bm{0})$. This gluon line creates a loop subgraph that is UV-divergent in four dimensions.

We subtract the divergence using the $\overline{\text{MS}}$ prescription, which creates dependence of $f_{a/\mathrm{p}}(\xi,\mu^{2}_{F})$ on the factorization scale $\mu^{2}_{F}$, possibly different from the renormalization scale $\mu^{2}_{R}$. In much of our discussion, we will assume that $\mu_{R}^{2}$ and $\mu_{F}^{2}$ are the same and denote both as $\mu^{2}$.

By examining the structure of the UV divergences, one finds that the functions $f_{a/\mathrm{p}}(\xi,\mu_{\mathrm{F}}^{2})$ obey DGLAP evolution equations,

[TABLE]

The functions $f_{a/\mathrm{p}}(\xi,\mu^{2}_{\mathrm{F}})$ are nonperturbative, but, since the dependence on $\mu^{2}_{\mathrm{F}}$ arises from the UV divergences of graphs for $f_{a/\mathrm{p}}(\xi,\mu^{2}_{\mathrm{F}})$, the evolution kernels $P_{ab}(z,\alpha_{\mathrm{s}}(\mu^{2}_{\mathrm{F}}))$ are perturbatively calculable as expansions in powers of $\alpha_{\mathrm{s}}(\mu^{2}_{\mathrm{F}})$:

[TABLE]

The exact evolution kernels $P_{ab}$ have been known up to three loops (NNLO) since 2004 Moch et al. (2004); Vogt et al. (2004), with active efforts now underway on computing the four-loop terms $P_{ab}^{(4)}$, cf. Ueda (2018). It is significant that the functions $P_{a\hat{a}}^{(n)}(z)$ do not depend on quark masses. In graphs for the PDFs, there are masses in quark propagators, ({\raise 0.97221pt\hbox to5.52084pt{\hss/\hss}\hbox to0.0pt{\hssk}}+m)/(k^{2}-m^{2}+{\mathrm{i}}\epsilon). However the ultraviolet poles of these graphs are determined by the behavior of the propagators for $k\to\infty$. In this limit, the masses do not contribute. This is an advantage of using the $\overline{\text{MS}}$ scheme for renormalizing the PDFs.

II.2 Infrared safety and factorization

PDFs are used to compute cross sections in collisions of a lepton with a hadron and collisions between two hadrons. We concentrate in this subsection on hadron-hadron collisions, since these are currently the subject of investigation at the Large Hadron Collider (LHC). Lepton-hadron collisions, as in deeply inelastic scattering, are simpler.

In one sense, the use of PDFs to describe proton-proton collisions is very simple. Suppose that we are interested in the cross section $d^{2}\sigma/(dp_{\mathrm{T}}dy)$, to produce a jet with transverse momentum $p_{\mathrm{T}}$ and rapidity $y$ plus anything else in the collision of a hadron of type $A$ and a hadron of type $B$. Or, suppose that we are interested in the cross section $d\sigma/dy$, to produce an on-shell Higgs boson with rapidity $y$ plus anything else. We can consider many cases at once by saying that we are interested in a cross section $\sigma[F]$ to measure a general observable $F$ that is “infrared-safe” in the sense that will be explained in subsection II.2.2. Then the PDFs relate $\sigma[F]$ to a calculated cross section $\hat{\sigma}_{ab}[F]$ for the collision of two partons. In its briefest form, the relation is

[TABLE]

Here we sum over the possible flavors $a$ and $b$ of partons that we might find in the respective hadrons. We integrate over the momentum fractions $\xi_{a}$ and $\xi_{b}$ of these partons. Then we multiply by $\hat{\sigma}_{ab}[F]$, which plays the role of a cross section for the collision of these partons to produce the final state that we are looking for.

We say that Eq. (25) expresses “factorization.”555The word “factorization” is applied to many formulas in which a physical quantity is expressed as a convolution of a product of factors. Eq. (25) is sometimes called inclusive collinear factorization. Factorization seems simple, but it is not. First, it works only when the observable to be measured, $\sigma[F]$, has a certain property, “infrared safety.” Second, factorization is approximate, and we need to understand what is left out. Third, its validity is not self-evident, as one finds when trying to calculate $\hat{\sigma}_{ab}[F]$ beyond the leading order and encountering infinities if the calculation is not carefully formulated. Fourth, while $\hat{\sigma}_{ab}[F]$ plays the role of a cross section to produce a certain final state in the collision of two partons, actually its calculation beyond the leading order in perturbative QCD involves subtractions, as discussed in subsection II.2.2 below.

II.2.1 Kinematics

Consider a hard scattering process in the collisions of two high energy hadrons, A and B. The hadrons carry momenta $P_{\mathrm{A}}$ and $P_{\mathrm{B}}$. The hadron energies are high enough that we can simplify the equations describing the collision kinematics by treating the colliding hadrons as being massless. Then with a suitable choice of reference frame, the hadron momenta are

[TABLE]

We then imagine a parton level process in which a parton from hadron A, with flavor $a$ and momentum $\xi_{\mathrm{a}}P_{\mathrm{A}}$ collides with a parton from hadron B, with flavor $b$ and momentum $\xi_{\mathrm{b}}P_{\mathrm{B}}$. This collision produces $m$ partons with flavors $f_{i}$ and momenta $p_{i}$. Each final state parton has rapidity $y_{i}$ and transverse momentum $\bm{p}_{i,\perp}$, so that the components of its momentum are

[TABLE]

Then momentum conservation gives us

[TABLE]

II.2.2 Infrared safety

In order for the factorization to work, the observable should be infrared-safe. The basic physical idea for this was introduced in Sterman and Weinberg (1977). We will follow the development in Kunszt and Soper (1992) and define what we mean by measuring the cross section for an observable $F$ and what it means for $F$ to be infrared safe. To keep the discussion simple, we temporarily assume that all of the partons involved are light quarks and gluons, which we consider to be massless, and that the observable does not distinguish the flavors of the partons. For this simple case, we express the parton-level cross section for an observable $F$ using the definition

[TABLE]

Here we start with the cross section to produce $m$ partons with momenta $\{p_{1},\dots,p_{m}\}$. We multiply the cross section by a function $F_{m}(p_{1},\dots,p_{m})$ that specifies the measurement that we want to make on the final-state partons. These functions are taken to be symmetric under interchange of their arguments. Accordingly, we divide by the number $m!$ of permutations of the parton labels. We integrate over the momenta of the final-state partons. The transverse momentum of parton 1 and the needed momentum fractions for the incoming partons are determined by Eq. (28). Finally, we sum over the number $m$ of final-state partons.

An example may be useful. If we want to evaluate the cross section for the sum of the absolute values of the transverse momenta of the partons to be bigger than 100 GeV, then we would choose a step function $F_{m}=\theta(\sum_{i=1}^{m}|\bm{p}_{i,\perp}|>100\ \mathrm{GeV})$. More generally, the $F_{m}$ for jet cross sections are made of step functions and delta functions constructed according to the jet algorithm used.

Infrared safety is a property of the functions $F_{m}$ that relates each function $F_{m+1}(p_{1},\dots,p_{m},p_{m+1})$ to the function $F_{m}(p_{1},\dots,p_{m})$ with one fewer parton. There are two requirements needed for $F$ to be infrared safe.

First, consider the limit in which partons $m+1$ and $m$ become collinear:

[TABLE]

Here $\tilde{p}_{m}$ is a lightlike momentum and $0\leq z\leq 1$. Therefore, $p_{m}+p_{m+1}\to\tilde{p}_{m}$. We can concentrate on just partons with labels $m+1$ and $m$ because the functions $F$ are assumed to be symmetric under interchange of the parton labels. In order for $F$ to be infrared safe, we demand that

[TABLE]

in the collinear limit (LABEL:eq:collinearlimit).

Second, consider also the limit in which parton $m+1$ becomes collinear to one of the beams:

[TABLE]

or

[TABLE]

Here $\lambda\geq 0$. When $\lambda=0$, parton $m+1$ is simply becoming infinitely soft. In order for $F$ to be infrared safe, we demand that

[TABLE]

in either limit (32) or (33).

Briefly, then, when the partonic masses are negligible, infrared safety means that the result of the measurement is not sensitive to whether or not one parton splits into two almost collinear partons, and it is not sensitive to any partons that have very small momenta transverse to the beam directions.

Sometimes an observable $F$ with this property is referred to as infrared and collinear safe (IRC-safe) instead of just infrared safe (IR-safe). The meaning is the same.

Cross sections to produce jets are infrared safe as long as we use a suitable algorithm, encoded in the functions $F_{m}$, to define jets. However, the cross section to produce a jet containing a specified number of charged particles would not be infrared safe: a collinear splitting of a gluon to a quark plus an antiquark increases the number of charged particles by two.

In Eq. (29), we have cross sections to produce specified partons. This is the formula that we use to construct $\hat{\sigma}_{ab}[F]$ on the right-hand side of Eq. (25). On the left-hand side of Eq. (25), we measure the cross section $\sigma[F]$ experimentally by using the same functions $F_{m}$ applied to hadrons instead of partons. Relating a hadron cross section to a parton cross section is evidently an approximation, so that there is an error in Eq. (25) no matter how many orders of perturbation theory we use in the calculation. The physical idea is that the infrared safe measurement involves a large scale $Q$, which is incorporated into the functions $F_{m}$. Processes that involve scales $M$ with $M\ll Q$ do not substantially affect the measurement. In particular, combining partons to form hadrons involves scales $M\sim 0.3\ \mathrm{GeV}$. Thus turning partons into hadrons does not substantially affect a measurement with a much larger scale $Q$.

While we can best understand the idea of IR safety by sticking to massless partons and flavor-independent measurement functions, it is important to keep in mind that a quark does have a mass, however small or large; QCD factorization need not assume that the partons are massless. For quark masses $m_{i}$ that are small compared to the scale $Q$ of the measurement that we have in mind, we should amend the definition of infrared safety of the functions $F_{m}$ to include taking the limit $m_{i}\to 0$. Sometimes a quark mass is comparable to the scale $Q$. For instance, this typically applies to top quarks. It also applies to, say, charm quarks when we take $Q$ to be only a couple of GeV. In such cases, we need not consider limits $m_{i}\to 0$ for “heavy” quarks and can operate with parton distribution functions only for the “light” quarks. See Sec. II.3 below.

II.2.3 Factorization

With the needed preparation accomplished, we can now state how the PDFs are used to calculate the cross section for whatever observable $F$ we want – as long as $F$ is infrared safe. For this condition to apply, the observable $F$ must be sufficiently inclusive. The formula we use was stated in Eq. (25), and we restate it here in a slightly more detailed form:

[TABLE]

Our convention in Eq. (35) is to use the name $\mu_{\mathrm{F}}$, the “factorization scale,” for the scale parameter in the parton distribution functions. The physical cross section cannot depend on $\mu_{\mathrm{F}}$, so $\hat{\sigma}$ must also depend on $\mu_{\mathrm{F}}$. Note that in Eq. (35), we have not yet expanded in powers of $\alpha_{\mathrm{s}}$.

The intuitive basis for Eq. (35) is very simple. The factor $f_{a/A}(\xi_{a},\mu_{\mathrm{F}}^{2})\,d\xi_{\mathrm{a}}$ plays the role of the probability to find a parton of flavor $a$ in a hadron of flavor $A$. For the other hadron, the corresponding probability is $f_{b/B}(\xi_{b},\mu_{\mathrm{F}}^{2})\,d\xi_{\mathrm{b}}$. Then $\hat{\sigma}[F]$ plays the role of a cross section to obtain the observable $F$ from the scattering of these partons, as given in Eq. (29). Naturally, this parton level cross section depends on the parton variables $a,b,\xi_{a},\xi_{b}$, as indicated by the subscript notation. Here the differential cross sections to produce $m$ final state partons contain delta functions that relate the momentum fractions $\xi_{\mathrm{a}}$ and $\xi_{\mathrm{b}}$ to the final state parton momenta, according to Eq. (28).

A similar formula applies for lepton-hadron scattering, the process providing important constraints on the PDF parameterizations (See Sec. II.4). Then there is only one hadron in the initial state, so the formula is simpler:

[TABLE]

The cross section $\hat{\sigma}[F]$ in Eq. (35) (or Eq. (36)) has a perturbative expansion in powers of $\alpha_{\mathrm{s}}(\mu_{\mathrm{R}}^{2})$, where the renormalization scale $\mu_{\mathrm{R}}^{2}$ can be chosen independently from the factorization scale $\mu_{\mathrm{F}}^{2}$. That is,

[TABLE]

Here $B$ is the integer that tells us how many powers of $\alpha_{\mathrm{s}}$ appear in the Born level cross section: e.g. 0 for Z boson production, 2 for two jet production. Perturbative calculations can be at lowest order (LO), corresponding to one term in the expansion, next-to-lowest order (NLO) with two terms, sometimes NNLO, and, in general, $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$.

When $\hat{\sigma}$ is expanded in powers of $\alpha_{\mathrm{s}}$, as in Eq. (37), our convention is to use the name $\mu_{\mathrm{R}}$, the “renormalization scale,” for the scale parameter in $\alpha_{\mathrm{s}}$. The physical cross section cannot depend on $\mu_{\mathrm{R}}$, so the functions $\hat{\sigma}^{B+k}$ in Eq. (37) must also depend on $\mu_{\mathrm{R}}$. (Other conventions for the precise meaning of $\mu_{\mathrm{F}}$ and $\mu_{\mathrm{R}}$ are possible.)

One useful property of Eqs. (36) and (37) is that the dependence of the calculated cross section on $\mu_{\mathrm{F}}^{2}$ and $\mu_{\mathrm{R}}^{2}$ diminishes as we go to higher orders. Indeed, the cross section in nature, $\sigma[F]$, does not depend on $\mu_{\mathrm{F}}^{2}$ and $\mu_{\mathrm{R}}^{2}$. Thus if we calculate to order $\alpha_{\mathrm{s}}^{B+k}$, the derivative of the calculated cross section with respect to $\mu_{\mathrm{F}}^{2}$ and $\mu_{\mathrm{R}}^{2}$ will be of order $\alpha_{\mathrm{s}}^{B+k+1}$. Because of this property, one often uses the effect of varying $\mu_{\mathrm{F}}$ or $\mu_{\mathrm{R}}$ by a fixed factor (e.g. 2 or 1/2) to provide an estimate of the error caused by calculating only to a finite perturbative order.

Note that, in order for terms of order $\alpha_{\mathrm{s}}^{B+k}$ in $\hat{\sigma}$ to match the terms for order $\alpha_{\mathrm{s}}^{k}$ in the evolution equation for the parton distributions, we need to include at least the terms up to $\alpha_{\mathrm{s}}^{k}P^{(k)}_{a\hat{a}}(z)$ in Eqs. (LABEL:eq:pdfevolution) and (II.1.6) giving the evolution of the PDFs. Since the lowest-order term in the evolution kernel is $P^{(1)}_{a\hat{a}}(z)$, including terms up to $P^{(k)}_{a\hat{a}}(z)$ is referred to as $\mathrm{N}^{k-1}\mathrm{L}\mathrm{O}$ evolution, while we say that including terms up to $\alpha_{\mathrm{s}}^{B+k}$ in the partonic cross section gives an $\mathrm{N}^{k}\mathrm{L}\mathrm{O}$ calculation. Thus, for example, if we have an NNLO cross section calculation, we will not obtain the proper matching for the $\mu_{\mathrm{F}}$ scale dependence unless we have at least NLO evolution for the PDFs. Normally, one uses yet higher-order evolution for the PDFs because this evolution determines the accuracy with which we know the PDFs at a high scale, given experimental inputs at much lower scales.

The error terms ${\cal O}\!\left({M}/{Q}\right)$ in Eqs. (35) and (36) arise from power-suppressed contributions that are beyond the accuracy of the factorized representation Collins (2013, 1998); Collins et al. (1983, 1985); Bodwin (1985); Collins et al. (1988). No matter how many terms are included in $\hat{\sigma}$, there are contributions that are left out. These terms are suppressed by a power of $M\sim 1\ \mathrm{GeV}$ divided by a large scale parameter $Q$ that characterizes the hard scattering process to be measured. Here we choose $1\ \mathrm{GeV}$ as a nominal scale for hadronic bound state physics. This value is somewhat smaller than the scale $\mu$ at which $\alpha_{\mathrm{s}}(\mu^{2})=1/2$, somewhat larger than the proton mass, and about four times larger than the inverse of the radius of a proton.

When data are included in a PDF fit but, for these data, $M/Q$ is not small enough to be completely negligible, it may be useful to include a nonperturbative model for the power-suppressed corrections. The model can have one or more parameters that can be fit to the data. For instance, one can add a contribution to the cross section that is proportional to $\kappa^{2}/Q^{2}$ where $\kappa$ is a parameter to be fit. This possibility is especially relevant for deeply inelastic scattering (Sec. II.4), for which one can make use of a theoretical expansion known as the operator product expansion Wilson (1969) to suggest the form of the first power-suppressed contribution. See such fits, for example, in Alekhin et al. (2012, 2017).

The power-suppressed contributions in Eq. (35) arise from the approximations needed to obtain the result. For instance, if a loop momentum $l$ flows through the wave function of quarks in a proton, we have to neglect $l$ compared to the hard momenta, say the transverse momentum of an observed jet. We can illustrate the issue of power corrections and, at the same time, see something about the interplay of ultraviolet and infrared singularities by looking at a very simple model for an integral that might occur in a calculation.

Consider a model integral for a correction to the Born cross section for a process involving only one initial state hadron. In this model, the Born cross section is

[TABLE]

where $l^{2}$ represents the transverse momentum of a quark that is restricted to be of order $1\ \mathrm{GeV}^{2}$ by a hadronic wave function $|\psi(l^{2})|^{2}$. The Born level PDF is $\int_{0}^{\infty}dl^{2}|\psi(l^{2})|^{2}$. The Born hard scattering function is simply 1. At order $\alpha_{\mathrm{s}}$, we should have an $\alpha_{\mathrm{s}}$ correction to the Born PDF times the Born hard scattering function, plus the Born PDF times an $\alpha_{\mathrm{s}}$ correction to the hard scattering function.

In our model, the order $\alpha_{\mathrm{s}}$ correction to the cross section is

[TABLE]

Here $Q^{2}$ is the “hard scale,” with $l^{2}\ll Q^{2}$. We supply a factor $\left[{\mu_{\mathrm{F}}^{2}}/{k^{2}}\right]^{\epsilon}$ that mimics dimensional regularization in $n=4-2\epsilon$ dimensions. This integral is, however convergent in the infrared (IR) and in the ultraviolet (UV), so we could simply set $\epsilon\to 0$ after its computation.

For $k^{2}\gg l^{2}$, we could approximate $1/(k^{2}+l^{2})\to 1/k^{2}$. This would define a part of the integral $I$ that represents the probability distribution $|\psi(l^{2})|^{2}$ times an order $\alpha_{\mathrm{s}}^{1}$ contribution to the hard scattering function,

[TABLE]

The integral is IR divergent. We have subtracted its IR pole, proportional to $-1/\epsilon$.

For $k^{2}\ll Q^{2}$, we could approximate $Q^{2}/(k^{2}+Q^{2})\to 1$. The integral would then define an $\alpha_{s}$ correction to the PDF times the Born hard scattering function,

[TABLE]

The PDF integral is UV divergent, so we have “renormalized” it by subtracting its UV pole, proportional to $+1/\epsilon$.

We can also use

[TABLE]

The $k$-integral equals zero because it equals $1/\epsilon-1/\epsilon$.

A simple calculation using the integrands of our integrals gives

[TABLE]

Thus our integral can be decomposed into the UV integral (including an IR subtraction), the IR integral (including its UV renormalization subtraction), and a remainder that is suppressed by a factor $l^{2}/Q^{2}$, where $l^{2}$ is of order $1\ \mathrm{GeV}^{2}$. This is the power-suppressed correction.

Notice that dimensional regularization with factorization scale $\mu_{\mathrm{F}}^{2}$ followed by subtracting $1/\epsilon$ poles serves two functions. The integral $I_{\textrm{UV}}$ in Eq. (40) is a model for how one calculates the hard scattering function at NLO. This integral has an IR divergence, which is removed by subtracting its pole. This subtraction matches the subtraction necessary to renormalize UV divergence in the correction to the model parton distribution function in Eq. (41).

This model calculation is somewhat misleading because it suggests that PDFs have a useful perturbative expansion. They do not. To compute a hard cross section in real QCD, one can replace the distribution of partons in a proton by the distribution of partons in a parton, which is singular but perturbatively calculable. See Collins (2013) for details.

Not much is known about the general form of the power corrections for hadron-hadron collisions.666See, however, Qiu and Sterman (1991a, b). It is important that they are there, but, if $Q$ is of order hundreds of GeV, then the power corrections are completely negligible. However, if $Q$ is of order 5 GeV, then we ought not to claim 1% accuracy in the calculation of $\sigma[F]$, no matter how many orders of perturbation theory we use.

Eq. (35), representing inclusive collinear factorization, is the basis of every prediction for hard processes at hadron colliders like the LHC, including both Standard Model processes and processes that might produce new heavy particles. So far as we know, it is a theorem for infrared-safe QCD observables dependent on energy-momentum variables that are of the same order of magnitude. There are other formulas in QCD that go under the name of “factorization” and typically apply either to amplitudes rather than cross sections or to observables dependent on several momentum variables of disparate orders of magnitude. These include $k_{\mathrm{T}}$ factorization and soft-collinear-effective-theory (SCET) factorization. These other forms of factorization may fail, thus requiring a more complex analysis, or at least are more subject to doubt than Eq. (35). See, for example, Collins and Qiu (2007); Catani et al. (2012); Forshaw et al. (2012); Rothstein and Stewart (2016); Schwartz et al. (2018).

Early attempts to establish inclusive collinear factorization Amati et al. (1978); Ellis et al. (1979) were instructive, but incomplete. Later proofs of Eq. (35) Collins et al. (1983, 1985); Bodwin (1985); Collins et al. (1988); Collins (1998) are far from simple. They could perhaps benefit from more scrutiny than they have received. One issue is that the published proofs have considered only the Drell-Yan process, not more complex processes like jet production. A more serious issue is that there is no known general method that can deal with the boundaries between integration regions in the Feynman diagrams. On the other hand, any breakdown in the inclusive collinear factorization in Eq. (35) could lead to infinities in calculations of $\hat{\sigma}$, and no problems have been observed so far even in $\mathrm{N}^{3}\mathrm{LO}$ calculations Anastasiou et al. (2015).

II.3 Treatment of heavy quarks

In order to accurately describe data at energies from one to thousands of GeV, modern global PDF fits not only change the number $N_{f}$ of active flavors depending on the scales $\mu_{\mathrm{R},\mathrm{F}}^{2}$, but also retain relevant quark mass dependence in the hard-scattering cross sections. The most comprehensive approach to do this is to work in one of the general-mass variable-flavor-number (GM-VFN) factorization schemes Aivazis et al. (1994); Buza et al. (1998); Chuvakin et al. (2000); Kramer and Spiesberger (2004); Kniehl et al. (2005b, a); Thorne and Roberts (1998); Thorne (2006); Forte et al. (2010). The massive fixed-flavor number schemes Gluck et al. (1982, 1988) are also applicable under the right circumstances and may result in simpler predictions. Such computations are a complex subject that we cannot cover in any depth in the space available. We will illustrate some of the key ideas by mostly following Krämer et al. (2000).

Consider the perturbative calculation of an infrared-safe cross section $\sigma$ with scale $Q^{2}$ when $Q^{2}\gg m_{i}^{2}$, where $i$ denotes any of the $\mathrm{u}$, $\mathrm{d}$, $\mathrm{s}$, $\mathrm{c}$, and $\mathrm{b}$ quarks. We can greatly simplify this calculation by neglecting masses of the five quarks. But what if $Q$ is high enough, and we expect that top quarks contribute either in the final state or in the virtual loop corrections? We rarely can set $m_{\mathrm{t}}=0$, as $m_{\mathrm{t}}\approx 174\ \mathrm{GeV}$ is so large that we seldom have $Q^{2}\gg m_{\mathrm{t}}^{2}$ even at the LHC.

There is a simple answer: in the range of energies comparable to $m_{t}$, we can use the $\overline{\text{MS}}$ scheme with five active quark flavors $\mathrm{u}$, $\mathrm{d}$, $\mathrm{s}$, $\mathrm{c}$, and $\mathrm{b}$. Top quarks are included in the relevant Feynman graphs, but, in accord with the CWZ prescription, we use the zero-momentum subtraction, instead of $\overline{\text{MS}}$ subtraction, for the UV renormalization of loop subgraphs with top-quark lines. Then, terms involving $m_{t}$ appear only in the hard cross section $\widehat{\sigma}$ of Eq. (35).

One consequence of this is that the evolution equation for $\alpha_{\mathrm{s}}(\mu^{2}_{\mathrm{R}})$ uses the 5-flavor beta-function. Radiative contributions to the top quark mass also need to be renormalized, see the later part of the subsection. The 5-flavor scheme introduces nonzero PDFs $f_{i/\mathrm{p}}(\xi,\mu^{2}_{\mathrm{F}})$ for $i\in\{\mathrm{g},\mathrm{u},\bar{\mathrm{u}},\mathrm{d},\bar{\mathrm{d}},\dots,\mathrm{b},\bar{\mathrm{b}}\}$. There are no parton distributions for $i=\mathrm{t}$ or $i=\bar{\mathrm{t}}$ in this scheme. The PDFs $f_{i/\mathrm{p}}(\xi,\mu^{2}_{\mathrm{F}})$ are defined as in the previous sections and evolve using 5-flavor DGLAP kernels. Another (not obvious) result is that top-quark contributions are negligible in the limit $Q^{2}\ll m_{\mathrm{t}}^{2}$: that is, top quarks decouple when the momentum scale $Q$ of the problem is much smaller than top-quark mass.

In a general case, we can distinguish among several versions of an $N_{f}$-flavor scheme. In the zero-mass (ZM) scheme, only the Feynman graphs with massless active quarks, $q_{i}$ with $i\leq N_{f}$, are included in $\widehat{\sigma}$ in the factorized hadronic cross section (35). In the fixed-flavor number (FFN) scheme, the massless active quark contributions for $i\leq N_{f}$ are included as in the ZM scheme, and also the Feynman graphs with massive inactive (anti)quarks, $q_{i}$ with $i>N_{f}$, are included only in the hard cross section $\widehat{\sigma}$. Finally, the most complete general-mass (GM) scheme retains non-negligible quark mass terms from both active and inactive quarks in all parts of the hadronic cross section.

When discussing the region $Q^{2}\gg m_{\mathrm{b}}^{2}$, we can use the 5-flavor FFN scheme and set masses of $\mathrm{u},\mathrm{d},\mathrm{s},\mathrm{c},\mathrm{b}$ quarks to zero. Consider now the region $Q^{2}\sim m_{\mathrm{b}}^{2}$ and $Q^{2}\gg m_{\mathrm{c}}^{2}$. Here, we can use the 4-flavor FFN scheme. There are parton distributions for flavors $\mathrm{u},\mathrm{d},\mathrm{s},\mathrm{c}$, and these quarks are considered massless, while there are no parton distributions for $\mathrm{b}$ and $\bar{\mathrm{b}}$. In fact, the 4-flavor scheme is also fine for the subregion of $Q^{2}\gg m_{\mathrm{b}}^{2}$ in which $\alpha_{\mathrm{s}}\log(Q^{2}/m_{\mathrm{b}}^{2})\ll 1$, where we do not need to sum logs of $Q^{2}/m_{\mathrm{b}}^{2}$.777One can also use a 3-flavor FFN scheme. For instance, this scheme is used for fitting $\mathrm{c}$- and $\mathrm{b}$-quark production in deeply inelastic scattering in Alekhin et al. (2017). It introduces PDFs for $\mathrm{u}$, $\mathrm{d}$, and $\mathrm{s}$ quarks, which are treated as massless. The $\mathrm{c}$ and $\mathrm{b}$ quarks appear, with their masses, in $\hat{\sigma}$.

We now have two possible FFN schemes for calculating a cross section at $m_{\mathrm{b}}^{2}/Q^{2}\ll 1$ but $\alpha_{\mathrm{s}}\log(Q^{2}/m_{\mathrm{b}}^{2})\ll 1$: the 5-flavor scheme with $\alpha_{\mathrm{s}}^{(5)}(\mu^{2})$ and $f_{i/\mathrm{p}}^{(5)}(\xi,\mu^{2})$; and the 4-flavor scheme with $\alpha_{\mathrm{s}}^{(4)}(\mu^{2})$ and $f_{i/\mathrm{p}}^{(4)}(\xi,\mu^{2})$. (Here, we set $\mu_{\mathrm{R}}^{2}=\mu_{\mathrm{F}}^{2}=\mu^{2}$.) The physical predictions must be the same, order by order in $\alpha_{\mathrm{s}}$, in either scheme. This condition gives us matching relations between $\alpha_{s}$ and PDFs in the 4-flavor and 5-flavor schemes. At lowest order in $\alpha_{s}$, these relations are very simple. We should not use $\alpha_{\mathrm{s}}^{(5)}(\mu^{2})$ and $f_{i/\mathrm{p}}^{(5)}(\xi,\mu^{2})$ for calculating physical cross sections unless $\mu^{2}\gg m_{\mathrm{b}}^{2}$, but if we simply use their analytic forms for $\mu^{2}=m_{\mathrm{b}}^{2}$, we have

[TABLE]

At higher orders of $\alpha_{s}$, $f_{\mathrm{b}/\mathrm{p}}^{(5)}(\xi,m_{\mathrm{b}}^{2})\neq 0$, and the matching conditions are different and depend on whether $m_{b}$ is an $\overline{\text{MS}}$ or pole mass Ablinger et al. (2017); Bierenbaum et al. (2007); Buza et al. (1996). Then, to obtain the $f_{i/\mathrm{p}}^{(5)}(\xi,\mu^{2})$ for $\mu^{2}>m_{\mathrm{b}}^{2}$, we solve the $N_{f}=5$ evolution equation with a boundary condition (LABEL:eq:4to5flavors) at $\mu^{2}=m_{\mathrm{b}}^{2}$. One can also use a different scale, $\mu_{\mathrm{b}}^{2}\neq m_{\mathrm{b}}^{2}$, as well as alternative presriptions, for the matching Bertone et al. (2017a, 2018); Kusina et al. (2013).

We can derive analogous matching relations between the 3-flavor and 4-flavor schemes at $\mu=m_{\mathrm{c}}\approx 1.3\mbox{ GeV}$. The full range $\mu^{2}\geq m_{c}^{2}$ is then described by a sequence of the schemes with $N_{f}=3,4,$ and $5$ that together comprise a VFN scheme.

The scheme described so far is conceptually simple, but involves awkward switches between different values of $N_{f}$ at yet unspecified energy values. For instance, suppose we wish to switch between the $N_{f}=3$ FFN calculation and the $N_{f}=4$ calculation at a switching value $\mu=Q\equiv\mu_{4}$ somewhere above $m_{\mathrm{c}}$. At the corresponding value of $Q=\mu_{4}$, the calculated cross section will be discontinuous, which is detrimental in a PDF fit. If $\mu_{4}^{2}$ is too close to $m_{\mathrm{c}}^{2}$, the $N_{f}=4$ FFN cross section will miss important terms proportional to $m_{\mathrm{c}}^{2}$. If $\mu_{4}^{2}$ is too far above $m_{\mathrm{c}}^{2}$, the missing higher-order collinear logs $\alpha_{\mathrm{s}}\log(Q^{2}/m_{\mathrm{c}}^{2})$ make the $N_{f}=3$ FFN cross section unreliable in the upper range of $Q$.

To avoid that, many global fits use some version of a general-mass variable flavor number (GM-VFN) scheme, see references above. These schemes allow for $N_{f}$ switching exactly at the corresponding quark mass, $\mu_{i}=m_{i}$, and they achieve a smooth interpolation between the switching points. In such a scheme, one retains numerically nonneglible masses for any quark type in $\hat{\sigma}$. For instance, the ZM VFN and GM VFN schemes differ only in the treatment of the terms of order ${\cal O}(m_{i}^{2}/Q^{2})$ in the short-distance cross sections. They have the same mass dependence of PDFs. All that we are doing is including the essential ${\cal O}(m_{i}^{2}/Q^{2})$ terms in $\hat{\sigma}$.

The logic that we just outlined is closely followed by the simplified Aivazis-Collins-Olness-Tung (S-ACOT) scheme Aivazis et al. (1994); Krämer et al. (2000). It is proved to all $\alpha_{s}$ orders in Ref. Collins (1998) and applied in DIS up to NNLO Guzzi et al. (2012) for use in CTEQ-TEA fits. In the ACOT family of schemes, the flavor number in $\alpha_{s}$, masses, and PDFs is varied according to the CWZ prescription. Other GM-VFN schemes are perturbatively equivalent to the (S-)ACOT scheme. The reader can find comparisons between the GM-VFN approaches in Gao et al. (2018); Guzzi et al. (2012); Alekhin et al. (2010); Binoth et al. (2010).

We have discussed variable flavor number schemes with up to five active flavors. Also, at a 100 TeV collider, introduction of a top-quark PDF may be warranted, see Sec. 3.4 in Mangano et al. (2017). Groups doing global fitting often also present PDFs determined in fixed-flavor number schemes.

Heavy-quark masses. All nonzero quark masses in the calculation require renormalization. Either the heavy-quark $\overline{\text{MS}}$ masses or pole masses provided by the Review of Particle Physics Tanabashi et al. (2018) are used as input parameters when fitting the PDFs. Their values can be even extracted from PDF fits Gizhko et al. (2017); Alekhin et al. (2017); Ball et al. (2016); Gao et al. (2013). The $\overline{\text{MS}}$ mass for charm quark is better defined in pQCD and more precisely constrained by the world data. On the other hand, some pQCD calculations in the global fits use the pole mass as the input. Perturbative relations to convert the $\overline{\text{MS}}$ mass into the pole mass, or back, are known to a high order in pQCD Chetyrkin et al. (2000).

Fitted charm. The PDF for charm quarks is of special interest. Consider first the most standard treatment. Suppose, hypothetically, that we have data for a cross section at a scale a scale $Q$ around $m_{\mathrm{c}}$, and that, in this cross section, the hard interaction produces a charm quark and antiquark in the final state. The 3-flavor scheme would be useful for predicting this cross section, containing nonperturbative parameterizations of $f_{i/\mathrm{p}}^{(3)}(\xi,\mu_{0}^{2})$ for $i\in\{\mathrm{g},\mathrm{u},\bar{\mathrm{u}},\mathrm{d},\bar{\mathrm{d}},\mathrm{s},\bar{\mathrm{s}}\}$. There is no PDF for charm quarks. If we now change to a 4-flavor scheme, we have PDFs $f_{i/\mathrm{p}}^{(4)}(\xi,\mu^{2})$ for $i\in\{\mathrm{g},\mathrm{u},\bar{\mathrm{u}},\mathrm{d},\bar{\mathrm{d}},\mathrm{s},\bar{\mathrm{s}},\mathrm{c},\bar{\mathrm{c}}\}$. These would be obtained by matching to the 3-flavor PDFs at scale $\mu^{2}=m_{\mathrm{c}}^{2}$. For the charm quark, the matching gives $f_{\mathrm{c}/\mathrm{p}}^{(4)}(\xi,m_{\mathrm{c}}^{2})=0$ (at leading order). Matching gives us what may be called perturbative charm because the charm PDF for $\mu^{2}>m_{\mathrm{c}}^{2}$ is produced solely by the perturbative matching and DGLAP evolution.

Note that the leading-order $3\rightarrow 4$ matching condition $f_{\mathrm{c}/\mathrm{p}}^{(4)}(\xi,m_{\mathrm{c}}^{2})=0$ can be traced to the assumption that the power-suppressed corrections of order ${M}^{2}/{m_{c}^{2}}$ to a hypothetical cross section at $Q^{2}\sim m_{\mathrm{c}}^{2}$ would be negligible. If we think that this is perhaps not the case, we could allow $f_{\mathrm{c}/\mathrm{p}}^{(4)}(\xi,m_{\mathrm{c}}^{2})$ to be non-zero. Then, solving the DGLAP evolution equation with the boundary condition at $\mu^{2}=m_{\mathrm{c}}^{2}$ would result in a possibly non-zero value of $f_{\mathrm{c}/\mathrm{p}}^{(4)}(\xi,m_{\mathrm{c}}^{2})$, which in turn will influence the prediction of experimental results at scales $Q^{2}$ larger than $m_{\mathrm{c}}^{2}$, even $Q^{2}\gg m_{\mathrm{c}}^{2}$. If we put these experimental results into the PDF fitting program, we can fit $f_{\mathrm{c}/\mathrm{p}}^{(4)}(\xi,m_{\mathrm{c}}^{2})$. This gives us fitted charm.

The logic of fitted charm just outlined implies that, if we do consider experiments at $Q^{2}\sim m_{\mathrm{c}}^{2}$, we should allow for fitted non-perturbative charm quark contributions to the predictions for these experiments. If we use a 4-flavor scheme with $m_{\mathrm{c}}\neq 0$ as described below, these contributions can come from $f_{\mathrm{c}/\mathrm{p}}^{(4)}(\xi,\mu^{2})$ for $\mu^{2}\sim m_{\mathrm{c}}^{2}$.

Representative Feynman diagrams that contribute to the fitted component of the charm PDF can be viewed in Fig. 3 of Hou et al. (2018). These diagrams introduce contributions of order ${\cal O}\!\left({M}^{2}/{m_{c}^{2}}\right)$ to the charm PDF. The available data can accommodate or even mildly prefer contributions $f_{\mathrm{c}/\mathrm{p}}^{(4)}(\xi,m_{\mathrm{c}}^{2})$ carrying up to about 1% of the net proton’s momentum Ball et al. (2016); Jimenez-Delgado et al. (2015); Hou et al. (2018).

II.4 Deeply Inelastic Scattering

We conclude the theory overview by a brief discussion of deeply inelastic lepton scattering (DIS), $\ell(k)+A(p_{A})\to\ell^{\prime}(k^{\prime})+X$, a very important class of processes for the determination of PDFs Devenish and Cooper-Sarkar (2004). Here $\ell$ and $\ell^{\prime}$ can be either electrons, neutrinos, or muons with specified momenta $k$ and $k^{\prime}$. Hadron $A$ is a proton, nucleus, or pion with momentum $p_{A}$; and $X$ stands for unobserved particles. The interaction between the leptons and hadron proceeds by an exchange of a virtual $\gamma^{*}$, $Z$, or $W$ boson with momentum $q=k-k^{\prime}$.

Not only were the measurements in DIS historically influential in the development of QCD, while diverse DIS data from HERA and fixed targets still serve as the backbone for global fits, but also projections Abdul Khalek et al. (2018); Hobbs et al. (2019) show that DIS data will continue to provide essential constraints on the PDFs in the high-luminosity LHC era.

It is conventional to define three Lorentz-invariant variables,

[TABLE]

“Deeply inelastic” means that $Q^{2}\gg 1\ \mathrm{GeV}^{2}$, and $x_{\it bj}$ is not too small or too close to 1. The use of the variable $x_{\it bj}$ was first suggested by James Bjorken, who proposed that the cross section would have simple properties in the deeply inelastic limit Bjorken (1969).

If only the final-state lepton is observed, one often writes the spin-independent cross section as a linear combination of three “structure functions,” $F_{1}(x_{\it bj},Q^{2})$, $F_{2}(x_{\it bj},Q^{2})$, and $F_{3}(x_{\it bj},Q^{2})$. To determine the $F_{i}$’s, one needs data from two c.m. energies, $\sqrt{s}$. Otherwise, an approximate assumption is needed to extract the $F_{i}$. The structure function $F_{3}$ is nonzero only if the current violates parity.

The structure functions can be written in terms of PDFs in the form

[TABLE]

where $M\sim 1\mbox{ GeV}$. Usually, one chooses $\mu^{2}=Q^{2}$. This is really just another form of Eq. (36). It offers one nice result. In $eA$ scattering with a virtual photon exchange from the electron, at lowest order in $\alpha_{s}$, $\hat{F}_{2}$ contains a delta function that sets $\xi\to x_{\it bj}$:

[TABLE]

where $\mathcal{Q}_{i}$ is the electric charge (in units of $e$) of partons of type $i$: $Q_{\mathrm{u},\mathrm{c},\mathrm{t}}=2/3$, and $Q_{\mathrm{d},\mathrm{s},\mathrm{b}}=-1/3$. Because of the charge factors $Q_{i}^{2}$, the neutral-current DIS via the photon exchange is about four times more sensitive to up-type (anti)quark PDFs than to down-type ones. It is more difficult for DIS to constrain the $d$-quark and especially $s$-quark PDFs, so the uncertainties on these flavors tend to be higher than for $u$ and $c$ PDFs, as we already saw in Fig. 1. We caution, however, that the simple result in Eq. (47) does not hold beyond the lowest order, and that the structure functions are not to be confused with PDFs.

Note that we have used $\xi$ for the momentum fraction argument of PDFs and $x_{\it bj}$ for the kinematic variable defined in Eq. (45). In the particle physics literature, one often sees the notation $x$ used for the PDF momentum fraction argument $\xi$. In fact, we will sometimes use $x$ with this meaning later in this review. We always use $x_{\it bj}$ and not simply $x$ for the kinematic variable defined in Eq. (45).

III Statistical inference in fitting the parton distributions

In this section, we derive the key statistical results relevant for the extraction of the PDFs from the experimental data. We use the simplest possible framework, in which experimental errors can be approximated as having a Gaussian distribution, and the theoretical predictions are approximated as linear functions of the parameters used to describe the parton distribution functions. This framework is sometimes called the Hessian method Pumplin et al. (2001), since a certain matrix called the Hessian matrix plays a prominent role. The Hessian approach is motivated by the observation that many essential features of the PDF fits are captured by assuming an approximately Gaussian behavior of the underlying probability distribution of the experimental data. The PDF functional forms can be determined (“inferred”) from the experimental data by applying Bayes’ theorem Alekhin (1999) as is briefly summarized in Sec. III.1. The Hessian approximation provides a simplified solution to the problem of Bayesian inference when the PDFs are well-constrained, and the deviations from the most likely solution for PDFs are relatively small. The only non-standard feature that we add is the inclusion of a set of parameters $R_{k}$ that represent the possibility that the theory, with an ideal choice of parameters for the PDFs, is not perfect and may not fit data exactly.

III.1 Bayes’ theorem

Fitting parton distributions to data involves accounting for the statistical and systematic errors in the data. Thus we will need a statistical analysis. For this, we use a Bayesian framework in this paper. The alternative is a frequentist framework, but we find that the Bayesian framework is simple to understand and makes us more aware of assumptions that are otherwise left obscure.

We begin with Bayes’ theorem. At its base, this is a simple matter of counting. Consider a population in which each individual can have one or both of two characteristics, $T_{1}$ and $D$. For a concrete example, the population might consist of people in California. $T_{1}$ might be “has a certain genetic marker”, while $D$ might be “tests positive for this genetic marker.” Denote by $P(T_{1})$ the probability that an individual has characteristic $T_{1}$. That is, $P(T_{1})$ is the number of individuals with characteristic $T_{1}$ divided by the total number of individuals. Similarly, let $P(D)$ be the probability that an individual has characteristic $D$. Denote by $P(T_{1}|D)$ the conditional probability that an individual that is known to have characteristic $D$ has characteristic $T_{1}$. That is, $P(T_{1}|D)$ is the number of individuals with both characteristics $T_{1}$ and $D$ divided by the number of individuals with characteristic $D$. Similarly, let $P(D|T_{1})$ the conditional probability that an individual that is known to have characteristic $T_{1}$ has characteristic $D$. With these definitions, we have Bayes’ theorem:

[TABLE]

By evaluating the likelihood $P(D|T_{1})$ according to an explicit prescription presented in the next section and by having estimates for $P(D)$ and $P(T_{1})$, we could infer the “posterior probability” $P(T_{1}|D)$ by rearranging the factors in Eq. (48):

[TABLE]

Now consider another characteristic $T_{2}$, for instance “does not have the genetic marker.” Then Bayes’ theorem gives us

[TABLE]

Note that $P(D)$ cancels in this equation. If we know the quantities on the right-hand side of Eq. (50), this tells us the relative probabilities of finding characteristics $T_{1}$ and $T_{2}$ among individuals that are known to have characteristic $D$.

Eq. (50) is directly applicable and useful in cases in which there is a large number of individuals, and we know the probabilities on the right-hand side of the equation. For instance, a physician of a patient who tests positive for a rare health condition will want to use Eq. (50) to help decide on prescribing a specific treatment to address this condition.

We will, however, use Eq. (50) in a more subtle case. Suppose that there is only one individual. Suppose further that you have a subjective belief, based on your prior experience, that the probability that this individual has property $T_{i}$ is $P(T_{i})$ for $i\in\{1,2\}$. This “prior probability” could be formed based on your past observations.

Now suppose that the individual is observed to have property $D$. Assume that you know how to compute the probability $P(D|T_{i})$ for an individual to have property $D$ if the individual has property $T_{i}$. You can turn this knowledge around with the help of Bayes theorem (49) to calculate the probability $P(T_{i}|D)$ for having property $T_{i}$ on the condition that the individual was observed to have property $D$. You can also use Eq. (50) to compute the “posterior probability ratio” $P(T_{1}|D)/P(T_{2}|D)$, by multiplying the prior probability ratio $P(T_{1})/P(T_{2})$, reflecting your knowledge before the measurement, by the ratio $P(D|T_{1})/P(D|T_{2})$ of the likelihoods calculated for each property $T_{i}$. Then, you declare an updated prior probability ratio:

[TABLE]

In our physics application, $T_{1}$ and $T_{2}$ are possible theoretical models describing a physical system, or perhaps one model of the system with two choices of parameters. Even when we cannot interpret $P(T_{i})$ by counting the number of instances of the system in which $T_{i}$ holds, we often have some idea about which theory is more likely, and we express this belief as the prior probabilities $P(T_{i})$. The prior probabilities are often based on previous experiments, but are partly subjective. Thus your prior probabilities may not be the same as mine. Now we make a measurement and observe that the system has a property $D$. With these new data in mind, we update the probability ratio for the two theories according to Eq. (51).

Your new estimate is still partly subjective and may not agree with mine. However, as data accumulates, it frequently happens that the subjective nature of our probabilities ceases to really matter. Suppose, for example, that your prior probability ratio is $P(T_{1})/P(T_{2})=10$, and mine is $P(T_{1})/P(T_{2})=0.1$. Then we do not agree which theory is more likely. However, suppose that, after lots of data become available, the likelihood ratio for the data is $P(D|T_{1})/P(D|T_{2})=10^{6}$. Then your posterior probability ratio is $P(T_{1}|D)/P(T_{2}|D)=10^{7}$, while mine is $P(T_{1}|D)/P(T_{2}|D)=10^{5}$. Either ratio is very large and strongly disfavors theory $T_{2}$ in comparison to theory $T_{1}$. At this point, I just agree that you were right, and we stop discussing the matter.

Notice that our subjective probability estimates change by the same factor, equal to the likelihood ratio $P(D|T_{1})/P(D|T_{2})$. We agree on this factor because we agree on the data $D$ as well as on the calculation. For the agreement to be possible, the theories must be natural, in the sense that the theory predictions $T$ and, by extension, the probabilities depending on them, are smooth functions of theoretical parameters. The technical reason is that estimation of correlated uncertainties requires inversion of a large Hessian matrix $H_{\alpha\beta}$ constructed from derivatives $\partial T_{k}/\partial a_{\alpha}$ of model predictions $T_{k}$ with respect to parameters $a_{\alpha}$. This inversion is numerically stable only if $\partial T_{k}/\partial a_{\alpha}$ are well-behaved. Naturalness is thus required for reliable estimates of the probabilities and derived quantities, including the model-discriminating ratios $P(D|T_{1})/P(D|T_{2})$ and uncertainties on theoretical predictions.

III.2 Gaussian model for the distribution of data

We wish to fit parameters $a_{\alpha}$, $\alpha\in\{1,\dots,N_{P}\}$, to data. The data are given by values $D_{k}$, $k\in\{1,\dots,N_{D}\}$. Each $D_{k}$ represents the number of counts divided by an integrated luminosity in a certain bin of measured momenta in a certain experiment that is to be included in the fit. In the Gaussian approximation, we suppose that, after accounting for experimental errors, the data have the form

[TABLE]

Here $\langle D_{k}\rangle$ is the value that the datum $D_{k}$ would have if there were no experimental uncertainties from either counting statistics or systematic effects such as detector calibration. The value $\sigma_{k}$ is the statistical uncertainty (one standard deviation quoted in the data tables with the “$\pm$” sign) associated with $D_{k}$. The variable $\Delta_{k}$ represents fluctuations in $D_{k}$ from counting statistics or other sources. In accord with our Gaussian approximation, the $\Delta_{k}$ are normalized to be independent Gaussian random variables with mean 0 and variance 1. We will use the notation $\mathcal{N}(0,1)$ for this distribution. In a typical particle experiment, there can be experimental systematic uncertainties, for example, associated with some imprecision in measurements of luminosity or particle energies. In general, the systematic uncertainties can be decomposed into two types: random fluctuations that are fully uncorrelated point-by-point and contribute to $\sigma_{k}\Delta_{k}$ in Eq. (52); and fully correlated variations between the experimental data points, contributing to the last term in Eq. (52). Then, we would interpret $\sigma_{k}$ as the full uncorrelated uncertainty, composed from the uncorrelated statistical and systematic uncertainties added in quadrature.

Let the number of sources of correlated systematic uncertainties be $N_{\lambda}$. We represent the correlated systematic uncertainties by a correlation matrix $\sigma_{k}\beta_{kJ}$. This is a matrix with indices $k,J$. There is no sum over $k$. Then in Eq. (52) there is a sum over an index $J$ that labels the sources of correlated systematic uncertainties. The contribution to $D_{k}$ from a systematic source $J$ is written as $\sigma_{k}\beta_{kJ}\bar{\lambda}_{J}$, where $\bar{\lambda}_{J}$ is an independent random variable with an $\mathcal{N}(0,1)$ distribution.888In general, both statistical and systematic uncertainties can be asymmetric and quoted as such by the experimental groups. In practice, neglecting the asymmetry has been often acceptable in the PDF analyses.

An example may be helpful. Each $D_{k}$ is obtained by dividing a number of counts by a measured value $\mathcal{L}_{0}$ of the integrated luminosity $\mathcal{L}$ for the corresponding experiment. There is some uncertainty in the luminosity measurement, which we express by writing $\mathcal{L}=\mathcal{L}_{0}+\mathcal{L}_{0}\sigma_{\mathcal{L}}\bar{\lambda}_{J_{\mathcal{L}}}$, where $\sigma_{\mathcal{L}}$ is an estimated fractional uncertainty in the luminosity, ${J_{\mathcal{L}}}$ is the index we choose for this source of uncertainty, and $\bar{\lambda}_{J_{\mathcal{L}}}$ is the corresponding Gaussian random variable. Then in Eq. (52),

[TABLE]

We will refer to the random variables $\bar{\lambda}_{J}$ as correlated systematic error variables. They are analogous to the variables $\Delta_{k}$, whose fluctuations embody the uncorrelated errors. In one way of analyzing the data, we do not try to determine the values of the $\bar{\lambda}_{J}$. In another analysis, introduced in Sec III.4, we introduce nuisance parameters $\lambda_{J}$ that are intended to estimate the values of the $\bar{\lambda}_{J}$.

The distribution of data $D$ is given by noting that, if $f$ is a function of the $\Delta_{k}$ and $\bar{\lambda}_{J}$, its expectation value is

[TABLE]

Using Eq. (52) in Eq. (54) and observing that $\langle\Delta_{k}\rangle=0$ and $\langle\bar{\lambda}_{J}\rangle=0$, we see that the expectation value of $D_{k}$ is the value $\langle D_{k}\rangle$ that appears in Eq. (52). Then, using Eq. (52) in Eq. (54) again, and noting that $\langle\Delta_{i}\Delta_{j}\rangle=\delta_{ij}$, $\langle\bar{\lambda}_{J}\bar{\lambda}_{L}\rangle=\delta_{JL}$, and $\langle\Delta_{i}\bar{\lambda}_{J}\rangle=0$, we find that

[TABLE]

where

[TABLE]

The matrix $C^{-1}$ is the inverse of a matrix $C$ that is called the covariance matrix. When experimental systematic errors $\beta_{iJ}$ are present, $C$ is generally not a diagonal matrix.999We may note that, if we start with the matrix $C^{-1}$, instead of Eq. (52), we can diagonalize $C^{-1}$ to obtain a non-unique representation in the form of Eq. (56), which is consistent with Eq. (52). We can use this simple result to show that, if $f$ is a function of the data $D_{k}$, then its expectation value is

[TABLE]

This gives the the probability density for the data $D_{k}$, with the variables $\bar{\lambda}_{J}$ integrated out.

The relation between Eq. (55) and Eq. (57) is an important general result that we can use whenever variables $y_{i}$ are distributed in a generalized Gaussian fashion, that is, with a probability density proportional to $\exp(-\sum y_{i}y_{j}M_{ij})$ for a matrix $M$. To prove it, we use Eq. (57) to compute $\langle(D_{i}-\langle D_{i}\rangle)(D_{j}-\langle D_{j}\rangle)\rangle$. We change variables in (57) to $x_{i}=(\sqrt{C})_{ij}(D_{j}-\langle D_{j}\rangle)$.101010We can define the matrix $\sqrt{C}$ because $C$ is a real symmetric matrix with all positive eigenvalues. We use $d^{N_{D}}D=(1/\sqrt{\det C})\,d^{N_{D}}x$. Then the exponent is $-\frac{1}{2}\sum_{i}x_{i}^{2}$. We then obtain $\langle x_{k}x_{l}\rangle=\delta_{kl}$, which gives us the result in Eq. (55).

III.3 Determining theory parameters using $\chi^{2}$

We now introduce theoretical predictions $T_{k}(a)$ for the data $D_{k}$. The prediction depends on some parameters $a$, most notably the parameters of the PDFs. We have defined $\langle D_{k}\rangle$ to be the value of $D_{k}$ if the experimental errors are negligible. For a given choice of the parameters $a$, the theoretical prediction may not be correct, but if the prediction is perfect then $T_{k}(a)=\langle D_{k}\rangle$.

If we substitute $\langle D_{k}\rangle\to T_{k}(a)$ in Eq. (57), we see that the probability to obtain the experimental results $D$ if the theory represented by $T(a)$ is correct is

[TABLE]

where

[TABLE]

and the data space measure is

[TABLE]

Consider two choices, $a_{1}$ and $a_{2}$, for the parameters. Suppose that before seeing the experimental results $D$, we judge the probability that theory $T(a_{1})$ is correct to be $P(T(a_{1}))$, and we judge the probability that theory $T(a_{2})$ is correct to be $P(T(a_{2}))$. Perhaps these prior probabilities are based on previous experiments, or perhaps they are based on some sort of dynamical model of parton behavior. Whatever our prior belief was, it should be modified after we know the experimental results. Let the new probabilities based on the experimental results $D$ be $P(T(a_{1})|D)$ and $P(T(a_{2})|D)$, respectively. Then Bayes’ theorem (50) gives us

[TABLE]

The information from experiment is contained in the likelihood ratio

[TABLE]

Thus $\chi^{2}(D,a)$ is the function that we need for parameter estimation: differences in $\chi^{2}$ give us the likelihood ratio that tells us how to adjust our judgments of which parameter choices are favored, using Bayes’ theorem (50). The maximum of the likelihood $P(D|T(a))$ with respect to $a$ is achieved at the global minimum of $\chi^{2}(D,a)$.

We can use $\chi^{2}(D,a)$ as a measure of goodness of fit: the theory matches the experimental results well when the differences $D_{i}-T_{i}(a)$ are small, meaning that $\chi^{2}(D,a)$ is small. However, in our opinion, there are reasons to be critical of $\chi^{2}(D,a)$ as the only measure of goodness of fit when the number $N_{D}$ of data is large. We will return to this question in Sec. IV.

III.4 Another definition of $\chi^{2}$

Given a sample of data $D_{k}$, we can measure how well the theory matches the data by the parameter $\chi^{2}(D,a)$ defined in Eq. (59). We can rewrite Eq. (59) as

[TABLE]

with the expectation value

[TABLE]

which we derived using Eqs. (59) and (55). This is minimized when $T_{i}(a)=\langle D_{i}\rangle$.

The function $\chi^{2}(D,a)$ is derived from the differential probability $P(D|T(a))$ to obtain the data $D$ if the theory $T(a)$ is correct, cf. Eq. (58). We can define another version of $\chi^{2}$, based on the differential probability $P(D|T(a),\lambda)$ to obtain the data $D$ if the theory $T(a)$ is correct, and if the random systematic error variables $\bar{\lambda}$ take the values $\bar{\lambda}_{J}=\lambda_{J}$. Using Eqs. (54) and (52) with $\langle D_{k}\rangle$ replaced by $T_{k}(a)$, this probability is

[TABLE]

Performing the integrations, we have

[TABLE]

where

[TABLE]

Then, as in Eq. (61), the information from experiment needed to apply Bayes’ Theorem to the determination of $a$ and $\lambda$ is contained in the likelihood ratio

[TABLE]

Note that the parameters $\lambda_{J}$ are not necessarily equal to the true systematic error variables $\bar{\lambda}_{J}$. Rather, the $\lambda_{J}$ are parameters that one can fit to the data $D_{k}$. The best fit values of $\lambda_{J}$ then approximate the true $\bar{\lambda}_{J}$. The $\lambda_{J}$ are called nuisance parameters. We maximize the likelihood of obtaining the observed data with parton parameters $a$ and nuisance parameters $\lambda$ by minimizing $\chi^{2}(D,a,\lambda)$ with respect to $a$ and $\lambda$.

The function $\chi^{2}(D,a,\lambda)$ is useful if we want to use data to estimate not only the true values of the parton parameters, $\bar{a}$, but also the systematic error parameters $\bar{\lambda}$. Since we are normally not so interested in the $\bar{\lambda}$ values, the function $\chi^{2}(D,a,\lambda)$ may seem less important than the function $\chi^{2}(D,a)$. However, the function $\chi^{2}(D,a,\lambda)$ has the advantage that it does not involve the covariance matrix $C_{ij}$.

It is significant that if we fit values of $\lambda$ by minimizing $\chi^{2}(D,a,\lambda)$, we can obtain $\chi^{2}(D,a)$. With the manipulations outlined in Appendix A, we can write $\chi^{2}(D,a,\lambda)$ in an instructive form:

[TABLE]

where $C$ is the covariance matrix defined in Eq. (56), $B$ is a matrix with elements

[TABLE]

and $\lambda^{\prime}$ is a shifted version of $\lambda$,

[TABLE]

The minimum of $\chi^{2}(D,a,\lambda)$ with respect to $\lambda$ occurs when $\lambda^{\prime}=0$, corresponding to $\lambda=\lambda^{\textrm{fit}}$. Thus

[TABLE]

The PDF-fitting groups use either form of $\chi^{2}$. See, for example, the review of various conventions for $\chi^{2}$ in Appendix A of Ball et al. (2013b).

In $\chi^{2}(D,a)$, the experimental systematic errors are encoded in the matrix $C$. This form is used, e.g., by the NNPDF analyses. In $\chi^{2}(D,a,\lambda)$, we have the systematic errors expressed explicitly using parameters $\lambda$. This is the convention adapted by CTEQ analyses, starting with CTEQ6 Pumplin et al. (2002). There is then an extra term $\sum_{J}\lambda_{J}^{2}$ in $\chi^{2}(D,a,\lambda)$. We are instructed to fit the parameters $\lambda$ to the data by minimizing $\chi^{2}(D,a,\lambda)$. In the following sections, we will see how to fit the theory parameters $a$ by minimizing $\chi^{2}(D,a)$ with respect to the parameters $a$. This is then equivalent to minimizing $\chi^{2}(D,a,\lambda)$ with respect to $a$ and $\lambda$.

If we use $\chi^{2}(D,a)$, then we do not need to be concerned with the systematic error parameters $\bar{\lambda}$. With $\chi^{2}(D,a)$, we have a matrix $C_{ij}$. The fact that this matrix is not diagonal indicates that the errors are correlated. The presence of $C_{ij}$ makes the formulas a little complicated, but there are no real conceptual complications: $C_{ij}$ acts as a metric tensor on the space of the data, so that one could think of $u_{i}C_{ij}v_{j}$ as simply $u\cdot v$.

The minimum of $\chi^{2}(D,a,\lambda)$ occurs at values $\lambda^{\textrm{fit}}_{J}$ of the nuisance parameters. What is the relation between the $\lambda^{\textrm{fit}}_{J}$ and the systematic error variables $\bar{\lambda}_{J}$? The correlated error variables $\bar{\lambda}_{J}$ influence the data $D_{k}$, but the uncorrelated error variables $\Delta_{k}$ also influence the $D_{k}$ and, furthermore, we do not know the exact parton parameters $a_{\alpha}$, so one cannot expect to be able to recover the $\bar{\lambda}_{J}$ exactly from the data. However, we will see below that the $\lambda^{\textrm{fit}}_{J}$ approximate the $\bar{\lambda}_{J}$ when there are many data $D_{k}$ and the parameters $\beta_{kJ}$ that give the influence of the $\bar{\lambda}_{J}$ on the data are not too small. Specifically, we will see that the $\lambda^{\textrm{fit}}_{J}$ approximate well the $\bar{\lambda}_{J}$ when the matrix elements $B_{IJ}$ are large. This happens when the sum over the data index $k$ in Eq. (69) includes many terms, and the parameters $\beta_{kI}$ are not too small.

The analysis is simple. We begin with $\chi^{2}(D,a,\lambda)$ in Eq. (66) and substitute

[TABLE]

from Eq. (52). This gives

[TABLE]

The partial derivatives vanish at the best fit $\lambda_{I}=\lambda^{\textrm{fit}}_{I}$,

[TABLE]

so

[TABLE]

On the right-hand side of this equation, the $\bar{\lambda}_{J}$ are of order 1, the quantities $T_{k}(a)-\langle D_{k}\rangle$ should be small if we use values of $a$ fit to the data, and, since the $\Delta_{k}$ are independent random variables with $\mathcal{N}(0,1)$ distributions, $\sum_{k}\beta_{kJ}\Delta_{k}$ should have fluctuations of order a typical $\beta_{kJ}$ coefficient times the square root of the number of contributing indices $k$. Thus the quantity in braces is not large. However, the matrix elements $B^{-1}_{IJ}$ are small. Thus we expect the $\lambda^{\textrm{fit}}_{I}-\bar{\lambda}_{I}$ to be small.

If we use $\chi^{2}(D,a,\lambda)$, then the treatment of the systematic error parameters $\lambda$ is similar to the treatment of the theory parameters $a$. We obtain values $\lambda^{\textrm{fit}}\approx\bar{\lambda}$. The values $\bar{\lambda}$ are, by assumption, distributed according to the $\mathcal{N}(0,1)$ distribution. Thus the values $\lambda^{\textrm{fit}}$ should be approximately distributed with this distribution. In Sec. IV.5, we use this to test the assumptions that we have used.

We can also consider the limit in which $\sum_{k}\beta_{kI}\beta_{kJ}$ is small rather than large for some value of $I$ and all values of $J$. Then we have $B_{IJ}\approx\delta_{IJ}$. This gives us

[TABLE]

Then $\lambda^{\textrm{fit}}_{I}$ is not approximately $\bar{\lambda}_{I}$ but is instead quite small. The test in Sec. IV.5 suggests that this happens for some of the nuisance parameters.

III.5 Dependence on the theory parameters

In section III.3, we introduced theory predictions $T_{k}(a)$ for the data $D_{k}$. The $T_{k}(a)$ depend on a number $N_{P}$ of parameters $a_{\alpha}$. Now we suppose that, in the neighborhood of the global $\chi^{2}$ minimum, the functions $T_{k}(a)$ are approximately linear in the parameters. In order to keep the notation as simple as possible, we define the origin of the parameter space so that the neighborhood of the global minimum is a region near $a=0$. That is, if we were fitting a function $x^{A_{1}}(1-x)^{A_{2}}$ for parameters $\{A_{1},A_{2}\}$ and preliminary fits gave $A_{1}\approx-1.3$, $A_{2}\approx 4.5$, we would define new parameters by $A_{1}=-1.3+a_{1}$ and $A_{2}=4.5+a_{2}$. Then we would be interested in small values of $\{a_{1},a_{2}\}$. Then we assume that, for the purpose of examining the fitting procedure, a linear approximation is adequate:

[TABLE]

Here and in what follows, we use the Einstein summation convention for parameter indices $\alpha$, $\beta$, ….

Finding the best-fit PDFs does not typically rely on the approximation (77) of linear dependence of $T_{k}(a)$ on the PDF parameters $a$. However, the PDF error analysis using the Hessian method does rely on this approximation (77). We will adopt it in this review. Using Eq. (77), $\chi^{2}$ is a quadratic function of the parameters $a$. In fact, within the range of $a$ important for the fit, this quadratic dependence is a pretty reasonable approximation, as one may note, for example in Fig. 11 in Sec. IV.9. One should be careful to check whether this approximation is adequate in an actual fit, see Sec. IV.

When the statistical and systematic experimental errors that contribute to the data in Eq. (52) vanish, the data $D_{k}$ equal their expectation values $\langle D_{k}\rangle$. We suppose that there are ideal values $\bar{a}$ of the parton parameters, related to the expectation values $\langle D_{k}\rangle$ of the data by

[TABLE]

We make the definition of the ideal parameters $\bar{a}$ more precise in Eq. (97) below. In Eq. (78), we have included constants $R_{k}$ that represent imperfections in the theory, such that even when we use parameters $\bar{a}$, the theory does not match $\langle D_{k}\rangle$ exactly. Of course, one commonly assumes that the $R_{k}$ are zero, but in this review we want to at least consider the possibility that something goes wrong. For example, the imperfections represented by $R_{k}$ could arise because we omitted higher-order contributions, there is beyond-the-standard-model physics in the data but not in the theory, or the parameterization that we use for the PDFs cannot match the true PDFs exactly.

Another way to include imperfections in the theory would be to incorporate theory errors into the analysis. In Eq. (52), we can set the expectation value $\langle D_{k}\rangle$ of the data in bin $k$ to the cross section in that bin calculated exactly in the Standard Model. We call this exact cross section $\overline{T}_{k}(\bar{a})$, so that Eq. (52) becomes

[TABLE]

where $\mathcal{E}$ is the set of experimental systematic errors. Now, we do not have the exact prediction $\overline{T}_{k}(\bar{a})$ available. All we have is the cross section $T_{k}(\bar{a})$ calculated, say, at NNLO (but with true PDF and related parameters $\bar{a}$). In the linear approximation that we use, we can parameterize our ignorance in the form

[TABLE]

where $\mathcal{T}$ is a set of sources of theory errors, and the parameters $\bar{\lambda}_{J}$ are random variables that are chosen from some distribution such as $\mathcal{N}(0,1)$. Then the term $\sigma_{k}\sum_{J\in\mathcal{T}}\beta_{kJ}\bar{\lambda}_{J}$ is our estimate of theory contributions to the true cross section $\overline{T}_{k}(\bar{a})$ that are not included in our NNLO calculation $T_{k}(\bar{a})$. The parameters $\beta_{kJ}$ represent our estimate of the dependence of $\overline{T}_{k}(\bar{a})$ on the $J$-th source of theoretical uncertainty. Of course, in reality only one value, $\lambda_{J}^{\textrm{true}}$, of each $\bar{\lambda}_{J}$ will be realized in an exact calculation. The $R_{k}$ in Eq. (78) are then

[TABLE]

When we combine Eqs. (79) and (80), we obtain

[TABLE]

where now both experimental systematic errors and our estimated theory errors are included. The representations of theory errors are discussed in Abdul Khalek et al. (2019b); Olness and Soper (2010); Cacciari and Houdeau (2011); Gao (2011); Forte et al. (2014); Harland-Lang and Thorne (2019). Until recently, PDF fits typically omitted theory errors. In this review, we do not include theory errors and instead represent imperfections in the theory by the $R_{k}$ in Eq. (78).

There is a certain freedom in the definition of $\bar{a}$ and $R_{k}$ in Eq. (78). We can use this freedom to simplify the later analysis. Suppose that we say that Eq. (78) applies for ideal parameters $\bar{a}^{(0)}$ and imperfection parameters $R^{(0)}_{k}$:

[TABLE]

Let $\bar{a}^{(0)}=\bar{a}+\delta a$, where the $\delta a_{\alpha}$ are small parameters that we are free to choose. Then

[TABLE]

This gives us Eq. (78) with $R_{k}=R^{(0)}_{k}-T_{k\alpha}\delta a_{\alpha}$. What should the $\delta a_{\alpha}$ be? In the analysis below, the vector $\sum_{kj}R_{k}C_{kj}T_{j\beta}$ plays an important role. This vector equals

[TABLE]

where

[TABLE]

is the Hessian matrix, which will play a major role in the subsequent analysis. We choose

[TABLE]

so that

[TABLE]

This is a useful property, as we will see.

We define $\chi^{2}(D,a)$ by Eq. (59), so that

[TABLE]

The minimum of $\chi^{2}(D,a)$ is at parameters such that

[TABLE]

This is

[TABLE]

where

[TABLE]

Thus the fit parameters are

[TABLE]

What is the expectation value of $a^{\textrm{fit}}_{\alpha}$? To answer this question, we need a result for $\bar{a}_{\alpha}$. Using Eq. (78), we have

[TABLE]

Thus

[TABLE]

Using the definition (86) of the Hessian matrix and the property (88) of the $R_{k}$, this is

[TABLE]

Thus

[TABLE]

Comparing to Eq. (93) gives us

[TABLE]

The result is that the expectation value of $a_{\textrm{fit}}$ is $\bar{a}$:

[TABLE]

This also gives the correlations of the parameters $a_{\alpha}$ with the data $D_{i}$ and with each other. From Eqs. (98) and (55), we have

[TABLE]

Thus

[TABLE]

For $\langle(a^{\textrm{fit}}_{\alpha}-\bar{a}_{\alpha})(a^{\textrm{fit}}_{\beta}-\bar{a}_{\beta})\rangle$, Eq. (98) gives

[TABLE]

Thus

[TABLE]

Compare this to Eq. (55), $\langle(D_{i}-\langle D_{i}\rangle)(D_{j}-\langle D_{j}\rangle)\rangle=C^{-1}_{ij}$. Thus we see the importance of the Hessian matrix: its inverse is the correlation matrix for the fitted parton parameters.

III.6 Distribution of the parameters

Since the data are Gaussian distributed, and the fit parameters are linearly related to the data, the fit parameters will be Gaussian distributed. The expectation values (103) give us the distribution of the best fit parameters, analogously to what we found in Eq. (57): given a function $f(a_{\textrm{fit}}-\bar{a})$, we have

[TABLE]

It is good to be clear about where this comes from. We consider an ensemble of repetitions of the experiments. As the data fluctuate, the values of $a_{\textrm{fit}}$ fluctuate according to the distribution (III.6). However, we can turn this around. Given the data $D_{k}$, we find the corresponding best-fit parameters $a_{\textrm{fit}}$. We do not know $\bar{a}$. But if we repeat the experiments many times, the values of the difference $(a_{\textrm{fit}}-\bar{a})$ will fluctuate around zero according to Eq. (III.6). Then Eq. (III.6) gives us a measure of the error in estimating $\bar{a}$ by $a_{\textrm{fit}}$.

Eq. (III.6) applies in the $N_{P}$-dimensional space of fit parameters. It is instructive to consider how this works in a particular coordinate system. We let $\{e^{(1)},e^{(2)},\dots,e^{(N_{P})}\}$ be a set of basis vectors for the parameter space. We take these basis vectors to be orthogonal and normalized using $H$ as the metric tensor111111 We do not distinguish between upper and lower indices $\alpha,\beta,\dots$. If we did, parton parameters would have upper indices, $a^{\alpha}$, and the metric tensor $H$ would have lower indices, $H_{\alpha\beta}$. Then $e^{(n)}$ and $H^{-1}$ would have upper indices, $e_{(n)}^{\alpha}$ and $(H^{-1})^{\alpha\beta}$.:

[TABLE]

The corresponding completeness relation is

[TABLE]

(To prove this, we define $\sum_{n}e^{(n)}_{\alpha}e^{(n)}_{\beta}=A_{\alpha\beta}$ and use Eq. (105) to show that $A_{\alpha\beta}H_{\beta\gamma}e^{(n)}_{\gamma}=e^{(n)}_{\alpha}$.)

It is often useful to choose the basis vectors to be the eigenvectors of $H$: $H_{\alpha\beta}e^{(n)}_{\beta}=h_{n}\,e^{(n)}_{\alpha}$. However, any choice of basis vectors obeying Eq. (105) will do. For instance, $e^{(1)}$ could be a vector normalized to $e^{(1)}_{\alpha}H_{\alpha\beta}\,e^{(1)}_{\beta}=1$, pointing in a direction that is of particular interest. Then the other $e^{(n)}$ could be chosen to satisfy Eq. (105). We will use this construction in Sec. III.7.

Using the basis vectors $e^{(n)}$, we can expand a general vector of parameters $a$ about $a_{\textrm{fit}}$ in the form

[TABLE]

Here the argument $t$ in $a(t)$ denotes $\{t_{1},\dots,t_{N_{P}}\}$. How does $\chi^{2}$ depend on the parameters $t_{n}$? To find out, we evaluate $\chi^{2}(D,a(t))$. Using Eqs. (86), (88), (89), (94), and (98), we obtain for a general choice of $a$,

[TABLE]

Here $a_{\textrm{fit}}$ is the parameter choice that we get from fitting the data $D_{k}$, $\bar{a}$ is what we would get by averaging $a_{\textrm{fit}}$ over an imagined ensemble of experiments, and $a$ represents the parameters that we are free to vary. Then if we substitute $a(t)$ given in Eq. (107) for $a$, we get

[TABLE]

That is, varying $a$ from $a_{\textrm{fit}}$ in any direction $e^{(n)}_{\alpha}$ by $t_{n}=1$ increases $\chi^{2}$ by 1.

The distribution (III.6) of differences of $a_{\textrm{fit}}$ from $\bar{a}$ can be rewritten in the $e^{(n)}$ basis. We define

[TABLE]

Then in Eq. (III.6), we can regard $f$ as a function of the eigenvector coordinates $t_{n}$ instead of the original parameters $a^{\textrm{fit}}_{\alpha}-\bar{a}_{\alpha}$, and we can change integration variables to the $t_{n}$, giving us

[TABLE]

Each coordinate $t_{n}$ follows an $\mathcal{N}(0,1)$ distribution.

In particular, if we are interested only in the component of $\bar{a}-a_{\textrm{fit}}$ in the direction $e^{(1)}$, we can let the function $f$ in Eq. (111) depend only on $t_{1}$. Then

[TABLE]

Comparing to Eq. (III.6), we see that a “$2\sigma$” value $t_{1}$, that is $t_{1}=2$, increases $\chi^{2}$ by 4 from its best fit value.

It is of interest to understand the distribution $\rho(R^{2},N_{P})$ of $R^{2}=\sum_{n}t_{n}^{2}$,

[TABLE]

This is the $\chi^{2}$ distribution with $N_{P}$ degrees of freedom. In fitting parton distributions, $N_{P}$ is quite large, say $N_{P}=25$. The mean value of any of the $t_{n}^{2}$ is $t_{n}^{2}=1$, but the mean value of $R^{2}$ is much larger, as we will see in Sec. III.8: $\langle R^{2}\rangle=N_{P}$. Using the $t_{n}$ as coordinates, the hypersurface $\sum_{n=1}^{N_{P}}t_{n}^{2}=R^{2}$ is a sphere. In terms of the original parton parameters $\bar{a}_{\alpha}-a^{\textrm{fit}}_{\alpha}$, it is an ellipsoid. With the values of $\bar{a}_{\alpha}-a^{\textrm{fit}}_{\alpha}$ distributed according to Eq. (III.6), in order for the ellipsoid to include 95% of the points we should choose $R\approx 6.1$ for $N_{P}=25$. In contrast, if we look at just one $t_{n}$, then in order for the interval $t_{n}^{2}<R_{n}^{2}$ to include 95% of the points $t_{n}$, we should choose $R_{n}\approx 2$. But with $N_{P}=25$, the fraction $P$ of points inside an ellipsoid with $R=2$ is $P\approx 5\times 10^{-7}$. This discussion illustrates that, when we discuss the uncertainties in the determination of the parton parameters, we need to carefully distinguish whether we are discussing an uncertainty interval in one dimension or in 25 dimensions. As another consequence of the large-$N_{P}$ geometry, when the PDF probability distribution is sampled by randomly varying the PDF parameters, the overwhelming majority of such Monte-Carlo parameter replicas are likely to be bad fits with $P\sim 0$ Hou et al. (2017b). Thus, though the estimates of the first and second moments of the $N_{P}$-dimensional probability distribution from the Monte-Carlo sample converge to their true values with about 100-1000 replicas, the Monte-Carlo method tends to be highly inefficient for exploring the neighborhood $\sum_{n}t_{n}^{2}\leq R^{2}$ with $R$ of order unity. In contrast, the analytic minimization of $\chi^{2}$ by the gradient descent, as implemented in CTEQ and MMHT fits, directly finds the neighborhood $\sum_{n}t_{n}^{2}\leq R^{2}$ around the global minimum and renders the probability distribution within this neighborhood. The analytic minimization and Monte-Carlo sampling approaches thus offer complementary strengths when examining the multidimensional probability distribution of PDF parameters and associated PDF uncertainties.

III.7 Calculation of a cross section

We can now ask another question. Suppose that $\sigma$ is a cross section that is determined by the PDFs. Then $\sigma$ is a function $\sigma(a)$ of the parton parameters $a$. We would need data with no errors to determine the ideal parameters $\bar{a}$, so we never know $\sigma(\bar{a})$ exactly. However, we can fit the parameters and estimate $\sigma(\bar{a})$ by $\sigma(a_{\textrm{fit}})$.

What is the expected error resulting from using the fit parameters? To analyze this, we begin by defining the shift in the cross section,

[TABLE]

Then with a linear approximation we write

[TABLE]

with $\sigma^{\prime}_{\alpha}\equiv\partial\sigma(a)/\partial a_{\alpha}$.

Now we define a special vector in the space of parameters according to

[TABLE]

This vector is normalized to

[TABLE]

It is useful to define more vectors, $\{e(\sigma)^{(2)},\dots,e(\sigma)^{(N_{P})}\}$ such that these vectors, together with $e(\sigma)\equiv e(\sigma)^{(1)}$, form a basis for the parameter space and such that the basis vectors are orthogonal and normalized using the metric tensor $H_{\alpha\beta}$ as in Eq. (105).

With the aid of these basis vectors, we can write a general parameter vector as

[TABLE]

The corresponding change in $\chi^{2}$ is, using Eq. (III.6),

[TABLE]

Let us set $a\to\bar{a}$ in the definition Eq. (118), so that

[TABLE]

Then, according to Eq. (111), as $a_{\textrm{fit}}$ varies in an ensemble experiment sets, the expansion parameters $\{t_{1},\dots,t_{N_{P}}\}$ fluctuate following $\mathcal{N}(0,1)$ distributions.

We can use this result to analyze the fluctuations in the cross section from Eq. (115):

[TABLE]

Using Eq. (116), this becomes

[TABLE]

Because of the orthonormality condition (105), only the first term survives and we obtain

[TABLE]

Thus the fluctuations in the cross section are given entirely by the fluctuations of the parameters along the special direction $e(\sigma)^{(1)}$. There is a coefficient, $[{\sigma^{\prime}_{\gamma}H^{-1}_{\gamma\delta}\sigma^{\prime}_{\delta}}]^{1/2}$ that is larger when the cross section is a fast varying function of the parton parameters. The remaining factor, $t_{1}$ fluctuates following an $\mathcal{N}(0,1)$ distribution as the data fluctuate. That means that if we want $\Delta\sigma$ to represent, say, a two standard deviation error on $\sigma$, we set $t_{1}=2$. Furthermore, $t_{1}$ has the property that when the parameters vary from the best fit parameters according to $a=a_{\textrm{fit}}+t_{1}\,e(\sigma)^{(1)}$, the $\chi^{2}$ increases by $t_{1}^{2}$.

Formulas to compute PDF uncertainties. There is a standard practical method for calculating $\Delta\sigma$. We choose basis vectors that are not specially adapted to the cross section $\sigma(a)$. We choose the basis vectors $e^{(n)}$, $n=1,\dots,N_{P}$ to obey the orthonormality condition Eq. (105). Typically, the basis vectors $e^{(n)}$ are chosen to be eigenvectors of the Hessian matrix $H$. Commonly, there are $2N_{P}$ error fits, $a_{\textrm{fit}}\pm\tilde{t}e^{(n)}$, that come with a set of published PDFs. Here $\tilde{t}$ is defined by the published analysis. If we assume that linear approximations are adequate, then we need only $N_{P}$ error fits, $a_{\textrm{fit}}+\tilde{t}e^{(n)}$, with positive $\tilde{t}$. Error fits with two signs Nadolsky and Sullivan (2001) allow for a more complete treatment, beyond what we give here, that allows for nonlinear contributions.

One can use the basis vectors $e^{(n)}$ to evaluate the uncertainty in $\sigma(a_{\textrm{fit}})$. For each direction $n$, define

[TABLE]

Recall from Eq. (III.6) that if we set $a=a_{\textrm{fit}}+te^{(n)}$, then $\chi^{2}(D,a)$ increases by $t^{2}$ compared to $\chi^{2}(D,a_{\textrm{fit}})$. As long as we use a linear approximation, we have

[TABLE]

Now sum the squares of the $\Delta\sigma^{(n)}$:

[TABLE]

Using the completeness relation (106) for the basis vectors $e^{(n)}$, this is

[TABLE]

According to Eq. (123), we can estimate the error in $\sigma$ by

[TABLE]

where, for example, we would choose $t=2$ if we want a “$2\sigma$” error estimate. Thus we can obtain $\Delta\sigma$ by using variations in the eigenvector directions $e^{(n)}$:

[TABLE]

That is, we need to calculate $N_{P}$ error contributions $\Delta\sigma^{(n)}$ by using the parton error sets according to Eq. (124). Then adding the errors $\Delta\sigma^{(n)}$ in quadrature gives the total error $\Delta\sigma$.

There is a second standard practical method for calculating $\Delta\sigma$. This method derives from the publications of the NNPDF group Ball et al. (2010, 2013c, 2015, 2017). With the NNPDF approach, there is effectively a very large number of parameters, and the distribution of the parameters is not strictly Gaussian. The distribution of results is represented by giving a large sample of parton distribution sets. Within the linear approximations that we use in this review, one would generate a large sample of parton distributions based on parameters $a_{\alpha}=a^{\textrm{fit}}_{\alpha}+\sum_{n}t_{n}\,e^{(n)}_{\alpha}$ as in Eq. (107), with Gaussian random variables $t_{n}$. Given this sample, one calculates $\sigma(a)$ for each example and, thus, obtains a corresponding sample of $\sigma(a)$ values, from which one obtains the statistical properties of the sample such as $\langle\sigma\rangle$ and $(\Delta\sigma)^{2}$.

We have spoken of $\sigma(a)$ as being a cross section. More broadly, $\sigma(a)$ in this section could be any physical quantity that depends on the parton parameters $a$. In particular, $\sigma(a)$ could be a parton distribution function $f_{a/A}(x,\mu^{2})$ for a particular parton flavor $a$, evaluated at a particular momentum fraction $x$ and a particular scale $\mu$. Then the calculation presented above gives us an error estimate $\Delta f_{a/A}(x,\mu^{2})$.

III.8 Expectation value and variance of $\chi^{2}$

In this section, we investigate the value of $\chi^{2}$ obtained in the fit. Start with Eq. (III.6) for $\chi^{2}(D,a)$. We fit the parameters $a$ to minimize $\chi^{2}(D,a)$ for given data $D$. Using the fit parameters gives

[TABLE]

Using Eq. (98) for $a^{\textrm{fit}}_{\alpha}-\bar{a}_{\alpha}$, this is

[TABLE]

where the matrix $M$ is

[TABLE]

Note that here we have eliminated the parameters $a_{\alpha}$ entirely.

We can use this to evaluate the expectation value $\langle\chi^{2}\rangle$ of $\chi^{2}$ and its variance

[TABLE]

We use Eq. (57) for the probability distribution of the data $D$. This gives

[TABLE]

and

[TABLE]

To derive this, one can change variables in Eq. (57) to $x_{i}=\sum_{j}(\sqrt{C})_{ij}(D_{j}-\langle D_{j}\rangle)$. Then the symmetries of the integrand imply that $\langle x_{i}x_{j}\rangle\propto\delta_{ij}$ and $\langle x_{i}x_{j}x_{k}x_{l}\rangle\propto\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}$. The coefficients of proportionality are simple to evaluate, giving Eqs. (134) and (135).

Now a certain amount of algebra with the matrices leads to results containing factors

[TABLE]

The results are

[TABLE]

and

[TABLE]

These results are often used to provide an indication of whether a good fit has been found. The parameters $R_{i}$ that we have introduced represent an imperfection in the theory: if the parton distributions do not have enough available parameters, we expect $R_{i}\neq 0$. If there are enough parameters and if the rest of the theoretical model is correct, then the $R_{i}$ should vanish. In that case, $\langle\chi^{2}\rangle$ should be close to $N_{D}-N_{P}$. For example, if $N_{D}=3000$ and $N_{P}=25$, then $\chi^{2}$ should be around 2975. The square root of the variance of $\chi^{2}$ in this case is $\sqrt{5950}\approx 77$. Thus we expect to find $\chi^{2}\approx 2975\pm 77$. The fit is bad if $\chi^{2}$ is too high or too low.

IV Tests of performance of the fit

In Sec. III.3, we examined the function $\chi^{2}(D,a)$ defined in Eq. (59) that is minimized to determine PDF parameters $a$ from experimental data $D$, giving values $a_{\textrm{fit}}$. However, this fitting procedure produces correct results only if the data $D$ are reliable within their errors as given by the experiments, and if the adopted theory is actually a good description of nature for some parameter combination $a_{\mathrm{fit}}$. If $\chi^{2}(D,a_{\textrm{fit}})$ does not lie within certain limits, one can conclude that something is wrong with the fit. However, it has been realized since the inception of the global QCD analysis in late 1980’s that the value of the global $\chi^{2}(D,a_{\textrm{fit}})$ is an essential, but far from sufficient, measure of the goodness of fit. [See, e.g., Morfin and Tung (1991)].

In this section, we argue that the PDF fit should pass a number of tests in order to fulfill what one might call a strong set of goodness-of-fit criteria. Several of these tests involve looking at quantities derived from the fitting procedure that should follow a predicted distribution if the statistical assumptions on which the fit is based are valid. One can then test whether the quantities are in fact distributed as predicted. We include the distribution of the nuisance parameters, the distribution of the residuals for the fitted data, and the distribution of $\chi^{2}$ values for many subsets of the data. Another test looks at whether individual subsets of the data are statistically consistent with the global $a_{\mathrm{fit}}$ in individual directions in the space of parameters. In Sec. IV.4, we discuss the closure test, which is a powerful test of the fitting methodology and is used, in particular, by NNPDF Ball et al. (2015).

These tests, taken together, are more constraining and difficult to satisfy than the standard weak goodness-of-fit criterion based on the value of the global $\chi^{2}(D,a)$. The point of the strong set of goodness-of-fit criteria is to find places where the statistical assumptions used in the fitting procedure break down. We will find, in examples, that some of the assumptions used in the fit do break down. When this happens, the PDF error estimates that emerge from the fit cannot be realistic. One possibility is that any observed problems are the result of understated estimates of the experimental systematic errors. We propose ways to deal with this in Sections IV.9 and IV.10. The method usually used to adjust the final PDF error estimates is to apply a tolerance criterion. This is a complicated subject, which we address briefly in Sec. IV.12.

One of the tests, which we investigate in Sec. IV.8, looks at $\chi^{2}$ for subsets of the data. Data subsets are often examined visually to rule out systematic discrepancies by comparing data and theory predictions in the figures. This is a reasonable, but slow and imprecise test. It can be made more precise using the procedure in Secs. IV.7 and IV.8.

A standard goodness-of-fit criterion based on the value of the overall $\chi^{2}(D,a_{\textrm{fit}})$ is called the hypothesis-testing criterion Collins and Pumplin (2001). We have seen in Eq. (137) that, if the theory is perfect so that the $R_{j}=0$, the expectation value of this quantity is $\langle\chi^{2}(D,a_{\textrm{fit}})\rangle=N_{D}-N_{P}$. However, Eq. (138) shows that $\chi^{2}(D,a_{\textrm{fit}})$ is expected to fluctuate by about $\sqrt{2(N_{D}-N_{P})}$. Thus we surely have a bad fit at the $2\sigma$ level if $\chi^{2}(D,a_{\textrm{fit}})-(N_{D}-N_{P})$ is bigger than $2\times\sqrt{2(N_{D}-N_{P})}$, that is, about 150 for $N_{D}-N_{P}\approx 3000$.

However, if we fit the parameters $a$ without changing the functional form, a small difference in $\chi^{2}$ values of order $2^{2}=4$ is already significant, while $2\times\sqrt{2(N_{D}-N_{P})}\approx 150$ is far too large. In this restricted situation, the parameter-fitting criterion, which assigns the 95% (or 68%) probability level to the increase $\Delta\chi^{2}=4$ (or 1), adequately estimates the uncertainty on parameters of this fixed model, as long as the statistical assumptions on which the fit is based are all valid. Multiple functional forms can give good fits to the data. To the parameter-fitting uncertainty, one must add the uncertainty due to the choice of the functional form of PDFs, as is discussed in Sec. IV.3.

As outlined above, we will look at quantities derived from the fitting procedure that should follow a predicted distribution. Call the quantities $q_{j}$. When testing for a possible systematic deviation from the predicted distribution for the $q_{j}$, we find it useful to transform the observed quantity to a form $x_{j}=x(q_{j})$ such that expected distribution of the $x_{j}$ is $\mathcal{N}(0,1)$, the standard normal (Gaussian) distribution with the mean of zero and variance of one. For a quantitative estimate of the probability that the observed distribution of the $x_{j}$ was sampled from ${\cal N}(0,1)$, we can apply the standard Anderson-Darling test Anderson and Darling (1952). The test yields a “distance”$A_{\textrm{obs}}$ of the observed cumulative probability distribution of $x_{j}$ values from that for ${\cal N}(0,1)$. Then it calculates the probability $P_{\textrm{A-D}}$ that the the same number of randomly drawn $x_{j}$ values from ${\cal N}(0,1)$ will have a distance $A$ with $A>A_{\textrm{obs}}$. With this test, a) $P_{\textrm{A-D}}$ always lies between 0 and 1, b) $P_{\textrm{A-D}}$ is close to 1 (or 0) if the histogram matches the ${\cal N}(0,1)$ distribution closely (or poorly), and c) if we repeat the sampling procedure many times with data actually drawn from the ${\cal N}(0,1)$ distribution, the values of $P_{\textrm{A-D}}$ will be uniformly distributed between 0 and 1.

We begin with a preliminary question: do we have enough fitting parameters to obtain a good fit to the data?

IV.1 Testing with resampled data

The PDFs $f_{a/p}(\xi,\mu^{2})$ must use a sufficiently flexible functional form to reproduce only regular, but no random, features of the hadronic data. However, the functional form for PDFs is known only semi-qualitatively based on considerations like the positivity of cross sections, asymptotic limits at small and large $x$, and nucleon sum rules. One resorts to a phenomenological form $f_{a/p}(\xi,\mu^{2}_{0})$ for the PDFs at the initial scale $\mu_{0}^{2}$ and must decide how many parameters $a_{\alpha}$ to use. If the number $N_{P}$ of parameters is too small, the theory may not be perfect. If too many, no global minimum of $\chi^{2}$ may exist, or we may overfit the data. In this section, we will explore the dependence on the number of free parameters for a given parameterization form. One can also vary the PDF functional forms, as discussed in Sec. IV.3.

To estimate the optimal number of parameters, let us return to $\chi^{2}$, using Eq. (98) in Eq. (III.6):

[TABLE]

Suppose that we obtain a set of parameters $(a_{\alpha})_{1}$ by minimizing $\chi^{2}(D_{1},a_{1})$ for the fitted data sample $D_{1}$ with $N_{D}$ data points. Then we use the same parameters to calculate $\chi^{2}(D_{2},a_{1})$ for a control data sample $D_{2}$ that is obtained by repeating the experiment with different random fluctuations. [$N_{D}$, $R_{i}$, $C_{ij}$ are the same for $D_{1}$ and $D_{2}$.] What do we get? From Eqs. (55), (99), and (103), the $\chi^{2}$ expectation for sample $D_{2}$, but using the parameters $a_{1}$ fitted to $D_{1}$, is given by

[TABLE]

This is bigger by $+2N_{P}$ than

[TABLE]

obtained by using the data sample $D_{1}$ to evaluate the fit $a_{1}$ obtained from this data. If there are more than just a few parameters, this is a big change.

Although having “too many” parameters makes $\langle\chi^{2}(D_{2},a_{1})\rangle$ larger, we could imagine that, with “too few” parameters, we get a bad fit because the fit cannot get close to the true PDFs. In our treatment here, with “too few” parameters, $\sum_{ij}R_{i}R_{j}C_{ij}$ is large.

We can illustrate this with a toy example, showing what may happen in parameterization studies with real parton distributions, such as the one in Martin et al. (2013). Suppose we fit a test function $h_{N_{P}}(x)$ with $N_{P}$ parameters to pseudodata generated by random fluctuations around the function

[TABLE]

For the fit, we use a polynomial,

[TABLE]

If we use $h_{3}(x)$, with just three parameters, we do not get a good fit, as we see in Fig. 2. Here for typical choices $x_{j}$ of $x$, the measures $R_{j}=f(x_{j})-h_{3}(x_{j})$ of how well the “theory” matches the exact function are typically of order 0.1. However, if we use 4 parameters to form $h_{4}(x)$, we can get a good fit, with $R_{j}\sim 0.01$. If we increase $N_{P}$ beyond 4, the measures $R_{j}$ are even smaller. Then increasing $N_{P}$ when fitting data, with its fluctuations, does not help produce a better fit to the true function $f(x)$. However, increasing $N_{P}$ does make $\chi^{2}$ smaller — because we start better fitting the random fluctuations!

To illustrate what happens, let us generate “toy data” $y_{i}=f(x_{i})+0.2f(x_{i})\,r_{i}$ at coordinate values $x_{i}=0.01i-0.05$ for $i=1,\cdots,N_{D}$, where $N_{D}=100$. These are shown as scattered points in Fig. 3(a). The $r_{i}$’s are random numbers sampled from ${\cal N}(0,1)$. Then fitting the data using $h_{4}(x)$ produces the black dashed curve in Fig. 3(a). This is already a fairly good fit to $f(x)$, with $R_{j}\sim 0.1$; but we can allow ourselves even more parameters, for example by fitting $h_{13}(x)$ with 13 parameters. The 13-parameter fit is shown as the solid blue curve in Fig. 3(a). This produces a smaller value of $\chi^{2}$, but not a better fit to $f(x)$.

It appears from Fig. 3(a) that what we are doing is fitting the fluctuations in the data. To test this, we can generate a second set of data $D_{2}$ using the same $f(x)$ and a different set of $r_{i}$. We measure $\chi^{2}(D_{2},a_{1})$ of the original fit $h_{13}(x)$ to the new data sample $D_{2}$. In Fig. 3(b), for each $N_{P}$, we repeat this procedure many times and show $\chi^{2}(D_{2},a_{1})$ averaged over many such trials as a function of $N_{P}$. We see that $N_{P}=3$ is not enough: there is a substantial decrease in $\langle\chi^{2}(D_{2},a_{1})\rangle$ if we increase $N_{P}$ to $4$. However, beyond $N_{P}=4$, $\langle\chi^{2}(D_{2},a_{1})\rangle$ increases with $N_{P}$, in agreement with Eq. (140). The rise of $\chi^{2}(D_{2},a_{1})$ for $N_{P}>4$ is suggestive of overfitting the randomly fluctuating data: while increasing $N_{P}$ improves $\chi^{2}(D_{1},a_{1})$ for the fitted sample, when $N_{P}$ is too large, it increases $\chi^{2}(D_{2},a_{1})$ for the control sample, indicating that the fit adapts to random fluctuations in $D_{1}$.

Of course, the optimal number of fitted parameters will depend on the size of the fluctuations in the data. Having 20% fluctuations in the data is not representative of what one finds in data used in PDF global fits. We can repeat this same exercise using 2% fluctuations. Then the minimum of $\chi^{2}(D_{2},a_{1})$ (as in Fig. 3(b)) occurs at $N_{P}=5$ instead of $N_{P}=4$.

IV.2 Dependence on the number of PDF parameters

A strategy comparing $\chi^{2}$ values of the fitted and control samples is routinely employed to prevent overfitting of the data in the approach utilizing neural-network (NN) parton distribution functions Ball et al. (2010, 2013c, 2015, 2017). The PDF of each flavor is given by a NN of a certain configuration, which behaves essentially as a very flexible function, with its optimal number of parameters selected so as to satisfy a two-fold condition that the resulting PDFs render acceptable fits to the fitted and control samples in each run, or replica, of the global analysis. Both the fitted sample $D_{1}$ and control sample $D_{2}$ are obtained by randomly fluctuating the central values of the data according to the Gaussian distributions provided by the standard deviations of the data. The fit consists in training the NN to maximize agreement with the fitted sample. When training the NN on sample $D_{1}$, $\chi^{2}(D_{1},a_{1})$ is improved to an arbitrary accuracy by training the neural network long enough. For the control sample $D_{2}$, the $\chi^{2}(D_{2},a_{1})$ initially decreases and then grows after some number of training cycles. The training is stopped when $\chi^{2}(D_{2},a_{1})$ starts growing. The NN obtained at this point most optimally approximates both $D_{1}$ and $D_{2}$ samples without overfitting $D_{1}$.

In the traditional approach used by the groups other than NNPDF, the PDFs are parameterized by a set of fixed functional forms; if the number of free parameters $N_{P}$ is too small or too large, the fit is too poor or unstable. For example, the recent MMHT fits approximate each PDF by a 6-degree Chebyshev polynomial Martin et al. (2013); Thorne et al. (2019). CT18 error PDFs Hou et al. (2019) are parameterized with functional forms containing Bernstein polynomials and up to five free parameters.

We can illustrate this behavior with an example using real data. The CT14HERA2 parton distributions Hou et al. (2017a) are parameterized using a generic form

[TABLE]

The $x^{A_{1}}$ and $(1-x)^{A_{2}}$ prefactors capture the typical behavior in the $x\rightarrow 0$ and $x\rightarrow 1$ limits, respectively. The function $P(x;A_{3},A_{4},...)$ is constructed as a linear combination of Bernstein basis polynomials,

[TABLE]

with $b_{n,m}(x)=\left({\begin{smallmatrix}n&\\ m&\end{smallmatrix}}\right)x^{m}(1-x)^{n-m}$. Up to four Bernstein polynomials per flavor are introduced in the CT14HERA2 parameterization, with a total of $N_{P}=26$ parameters. However, we can try to choose $N_{P}\leq 26$ or $N_{P}>26$. The number $N_{P}$ can be easily varied by adding or removing Bernstein polynomials with non-zero coefficients in the functions $P(x;A_{3},A_{4},...)$ in $f_{a}(x,Q_{0})$.

We divide the CT14HERA2 data set into two equal parts, assigning each datum to the half set $D_{1}$ or half set $D_{2}$ at random. Then we fit the $N_{P}$ parameters to data set $D_{1}$, giving parameters $a_{1}$. We measure $\chi^{2}(D_{1},a_{1})$ for this fitted data using the data $D_{1}$ to which it was fitted. We also measure the $\chi^{2}(D_{2},a_{1})$ for the fitted parameters $a_{1}$ using the second half set $D_{2}$. Alternatively, we fit the $N_{P}$ parameters to data set $D_{2}$, giving parameters $a_{2}$. We measure $\chi^{2}(D_{2},a_{2})$ for this fitted data using the data $D_{2}$ to which it was fitted. Then we also measure the $\chi^{2}(D_{1},a_{2})$ for this fit using the other data set, $D_{1}$.

We show the results in Fig. 4. As we increase the number of parameters, $\chi^{2}(D_{1},a_{1})$ and $\chi^{2}(D_{2},a_{2})$ decrease. For $N_{P}<28$, $\chi^{2}(D_{2},a_{1})$ and $\chi^{2}(D_{1},a_{2})$ also decrease, although there is not much decrease beyond $N_{P}=26$. For $N_{P}>28$, the behavior is different. For the fitted samples, the $\chi^{2}$ values continue to decrease, although the numerical minimization by the MINUIT program James and Roos (1975) becomes less stable. On the other hand, $\chi^{2}(D_{2},a_{1})$ and $\chi^{2}(D_{1},a_{2})$ exhibit large fluctuations for $N_{P}>28$ and, additionally, their values increase dramatically for the largest $N_{P}$ in the figure.

We could also examine $\chi^{2}(D,a)$ for the whole data set, $D=D_{1}+D_{2}$. In this case, $\chi^{2}(D,a)$ also initially decreases as $N_{P}$ increases, but for large enough $N_{P}$, the fits become unstable: first, the $\chi^{2}$ is no longer a well-behaved function of PDF parameters; for even higher $N_{P}$, a unique best-fit PDF set itself may not be found. This behavior is an indication of overfitting. The number $N_{P}$ beyond which overfitting happens depends on the functional forms of PDFs and other factors. With the resampling exercise, we detect traces of overfitting in a half of the full data set. The presented exercise and other parameterization studies indicate that, with the parameterization forms adopted in the latest CT PDFs, about 25-30 free parameters is optimal for describing the fitted data set.

A way to understand the results in Fig. 4 is to note that the analysis in Sec. IV.1 of $\langle\chi^{2}\rangle$ was based on the assumptions (from Sec. III.5) that one needed to consider only a small region of the parameter space, and that, within this region, the theory predictions $T_{k}(a)$ are approximately linear functions of the parameters $a$. The analysis then predicts a slow change of $\chi^{2}(D_{2},a_{1})$ (and $\chi^{2}(D_{1},a_{2})$) according to Eq. (140), as we increase the number of parameters $N_{P}$. This is not what we see in Fig. 4.

To help interpret this, suppose that we use the data $D_{1}$ and determine fit parameters $a_{1}$ from $D_{1}$. We see in Fig. 4 a large increase of $\chi^{2}(D_{2},a_{1})$ and large fluctuations for $N_{P}>28$. What apparently happens is that, when there are too many parameters, the data do not strongly constrain the components of $a$ in one or more “flat” directions in the parameter space. The value of $\chi^{2}(D_{1},a)$ changes little as $a$ varies in these directions. One needs to move through a wide range of $a$ in a flat direction before $\chi^{2}(D_{1},a)$ increases significantly. In this wide range, the $T_{k}(a)$ are not approximately linear functions of $a$, and $\chi^{2}(D_{1},a)$ is not approximately quadratic. Then $\chi^{2}(D_{1},a)$ may have multiple minima at locations that correspond to quite unreasonable behavior of the PDFs. The best fit parameters $a_{1}$ can be at an unphysical location that strongly depends on the fluctuations in the data $D_{1}$. We can test how reasonable the fit $a_{1}$ is by testing it against the second data set, $D_{2}$. With too many parameters, we can get the large values of $\chi^{2}(D_{2},a_{1})$ that we see in Fig. 4. We conclude that, if one uses a family of functions for the PDFs that have a variable number of parameters, then a simple test like that shown in Fig. 4, ideally repeated for several random partitions of the data in two halves, can estimate how many parameters one should allow.

IV.3 Dependence on the PDF functional form

There is another sort of test available. The PDFs are unknown functions, but one represents them using fixed functional forms with a finite number of unknown parameters.121212The neural net approach avoids this limitation. PDF parameterization studies thus constitute an essential step of the global analyses. In the previous section, we have examined the dependence of the fit results on the number of parameters within one general functional form. Evidently, one can also change the general functional form. Thus one should examine the tests listed above and later in this section for a large class of PDF functional forms.

Universal approximation theorems Cybenko (1989); Hornik et al. (1990); Hornik (1991) demonstrate that feed-forward neural networks with a single or multiple layers can approximate any continuous function and its derivatives to arbitrary accuracy. This fundamental result motivates the Neural Network PDF approach that parameterizes parton distribution functions by neural networks trained on the fitted data.

Alternatively, modern parameterizations based on Chebyshev or Bernstein polynomials also allow one to parameterize a variety of functional behaviors Pumplin (2010); Glazov et al. (2011); Martin et al. (2013); Gao and Nadolsky (2014); Hou et al. (2019). Suppose that two choices for functional forms give good fits in the sense of reasonably meeting the array of goodness-of-fit criteria discussed in this section, but they also give two sets of PDFs that differ outside of their parameter-fitting error estimates. This would suggest that there is a theoretical systematic error associated with the choice of PDF functional form that needs to be added to the analysis. In a global analysis such as CT18, multiple parameterization forms are tried to estimate this source of PDF uncertainty. For example, the green solid lines in Fig. 5 illustrate predictions for the best-fit gluon PDFs obtained with about 100 candidate CT18 fits based on different parameterization forms. Notice that the solid lines show large deviations from each other for large $x$ and small $x$, regions in which the gluon distribution is not well constrained by data. The PDF uncertainty of the published CT18 NNLO set, shown as a blue band, is computed using a prescription that is broad enough to cover variations due to the choice of the parameterization form. For details, see Hou et al. (2019).

IV.4 Closure test

Given a PDF functional form depending on parameters $a$, PDF fitting of data determines best fit parameters $a_{\textrm{fit}}$. The fitting procedure is based on Eq. (58), which asserts that if the theory represented by parameters $\bar{a}$ is correct, then, in an ensemble of trials, the data $D_{i}$ will be distributed with a probability $P(D|T(\bar{a}))$ that involves the theory predictions $T_{i}(\bar{a})$ and the covariance matrix and the covariance matrix $C_{ij}$. Then the parameters $a^{\mathrm{fit}}$ obtained from this data will be distributed according to Eq. (III.6). We can test this. We can generate a set of pseudodata $D^{\prime}_{k}$ using Eq. (52) with the truth values $\langle D_{k}^{\prime}\rangle$ set equal to theory predictions $T_{i}(a)$ for some realistic parameters $a$, and with $\sigma_{k}$ and $\beta_{kJ}$ obtained from the experimental groups. Then we can run the fitting procedure with the generated data $D^{\prime}_{k}$, giving a new fit $a^{\prime}_{\textrm{fit}}$. We should find that the new fit $a^{\prime}_{\textrm{fit}}$ agrees with the original parameters $a$ within the errors generated by the fit. If they do not agree, we should understand why. This “closure test” is not commonly carried out by fitting groups who use the Hessian method described in this review, but it is a feature of fits by the NNPDF group Ball et al. (2015).

IV.5 Test of the nuisance parameters

We now describe a goodness-of-fit test based on the distribution of the nuisance parameters. For this test, it is useful to use the form of $\chi^{2}$ in which nuisance parameters $\lambda_{J}$ appear explicitly, $\chi^{2}(D,a,\lambda)$ as given in Sec. III.4, Eq. (66). Then we can fit values $\lambda_{\textrm{fit}}$ of $\lambda$ by minimizing $\chi^{2}(D,a,\lambda)$. The minimum value is $\chi^{2}(D,a)$, expressed in terms of the covariance matrix $C_{ij}$ in Eq. (59). We have $\chi^{2}(D,a_{\textrm{fit}},\lambda_{\textrm{fit}})=\chi^{2}(D,a_{\textrm{fit}})$. We also argued after Eq. (75) that $\lambda^{\textrm{fit}}_{J}\approx\bar{\lambda}_{J}$ in an accurate fit with enough data. Since the $\bar{\lambda}_{J}$ are independent random variables distributed according to ${\cal N}(0,1)$, the $\lambda^{\textrm{fit}}_{J}$ are expected to follow the ${\cal N}(0,1)$ distribution, too. We can test these assumptions by making a histogram of $\lambda^{\textrm{fit}}_{J}$.

We show in Fig. 6 a histogram of the best fit nuisance parameters $\lambda^{\textrm{fit}}_{J}$ for the CT14 HERA2 NNLO fit along with a dashed, red curve showing the expected Gaussian distribution. The correlated systematic errors are implemented in the CT14HERA2 fit using the matrices $\beta_{kJ}$ provided with the experimental data sets. Evidently, the observed distribution is substantially narrower than the expected distribution. The mean for the observed distribution is $-0.06$, which is very close to 0, but its standard deviation is 0.8, noticeably smaller than 1. We also show a blue, solid curve giving a Gaussian distribution with this mean and standard deviation. This curve also does not match the observed distribution well.

For a quantitative estimate of the probability that the observed distribution was sampled from ${\cal N}(0,1)$, we can apply the standard Anderson-Darling testAnderson and Darling (1952) described at the beginning of Sec. IV. For the histogram in Fig. 6, we find that $P_{\textrm{A-D}}\sim 10^{-6}$. This indicates that it is very unlikely that the $\lambda_{J}^{\textrm{fit}}$ values were generated from the expected ${\cal N}(0,1)$ distribution, as was self-evident from the figure.

There are a few more $|\lambda^{\textrm{fit}}_{J}|$ that are larger than 2 than would be expected. However, the main feature that we see in Fig. 6 is that too many $|\lambda^{\textrm{fit}}_{J}|$ are very small, indicating possibly that the corresponding $\beta_{kJ}$ values are overestimated. This could indicate that the estimates are conservative in the sense that, if one suspects that the $\beta_{kJ}$ values for a source $J$ of systematic error should be smaller but one cannot prove it with solid evidence, then the conservative approach is to leave these $\beta_{kJ}$ values unchanged.

There is another reason that the values of some of the $\lambda_{J}^{\mathrm{fit}}$ might be smaller than expected. The expectation from Eq. (75) that $\lambda_{J}^{\mathrm{fit}}\approx\bar{\lambda}_{J}$ was based on the assumption that, for each data set $E$, there are many data and few nuisance parameters. What happens if this is not the case, as suggested near the end of Sec. III.4? Suppose that we add a spurious nuisance parameter $\lambda_{L}$ with $\beta_{kL}=0$ for all data indices $k$. Then, in Eq. (66) for $\chi^{2}(D,a,\lambda)$, the spurious $\lambda_{L}$ does not contribute to the first term, which involves data, and contributes only to the “penalty” term $\sum_{J}\lambda_{J}^{2}$. Then minimizing $\chi^{2}(D,a,\lambda)$ with respect to the $\lambda_{J}$ gives $\lambda_{L}^{\mathrm{fit}}=0$. If a certain data set $E$ has many systematic errors and thus many nuisance parameters, it may be that some of linear combinations of them are spurious in this sense. Then the same linear combinations of the fitted values $\lambda_{J}^{\mathrm{fit}}$ will vanish, so that the fitted values do not spread out over the whole $\lambda_{J}$ space but only a subspace. In this case, the absolute values of these $\lambda_{J}^{\mathrm{fit}}$ will be smaller than would be given by an $\mathcal{N}(0,1)$ distribution. In the case that there are many nuisance parameters, some of which have little effect, one can approximate the matrix $C^{-1}$ in Eq. (56) so as to leave effectively fewer, but better behaving, nuisance parameters, as described in Appendix B of Hou et al. (2019).

We suggest that, when performing a fit, it is useful to check the $\lambda^{\textrm{fit}}_{J}$ distribution. Values of $\lambda^{\textrm{fit}}_{J}$ that are smaller than expected from an $\mathcal{N}(0,1)$ distribution suggest problems with the systematic errors used in the fit but are perhaps not of great concern. An excessive number of large $\lambda^{\textrm{fit}}_{J}$ values could indicate a more serious difficulty that needs to be resolved.

IV.6 Test of data residuals

Eq. (66) gives $\chi^{2}(D,a_{\textrm{fit}},\lambda_{\textrm{fit}})$ as

[TABLE]

where

[TABLE]

is called the residual for datum $k$ obtained in the fit. Using Eqs. (52), (77), and (78) with $R_{k}=0$, this is

[TABLE]

The $\Delta_{k}$, introduced in Eq. (52), follow independent $\mathcal{N}(0,1)$ distributions. The values $(a^{\textrm{fit}}_{\alpha}-\bar{a}_{\alpha})$ and $(\lambda^{\textrm{fit}}_{I}-\bar{\lambda}_{I})$ have expectation values zero, but they have fluctuations that arise from the fluctuations in the data. If there are enough data, we expect the fluctuations of $(a^{\textrm{fit}}_{\alpha}-\bar{a}_{\alpha})$ and $(\bar{\lambda}_{I}-\lambda^{\textrm{fit}}_{I})$ to be small. Then the residuals $r_{k}$ should also be approximately distributed as ${\cal N}(0,1)$.

We can test these assumptions by making a histogram of the values $r_{k}$ obtained. We show in Fig. 7 a histogram of the residuals $r_{k}$ for the CT14 HERA2 NNLO fit along with a curve showing the expected Gaussian distribution. Comparing these using the Anderson-Darling test gives $P_{\textrm{A-D}}=5.7\times 10^{-3}$. Thus we can conclude with some confidence that the observed distribution of residuals was not drawn from exactly ${\cal N}(0,1)$. However, we judge the difference between the two distributions to be not physically significant. After all, we expect the observed distribution to be only approximately an ${\cal N}(0,1)$ distribution. The mean for the observed distribution is $0.04$, which is very close to 0, and its standard deviation is 1.04, which is quite close to 1. We also show a blue, solid curve giving a Gaussian distribution with this mean and standard deviation.

The distribution of residuals is another indicator that should be checked when performing a fit. In this case, no large discrepancies are observed.

IV.7 Value of $\chi^{2}$ from an individual experiment

In this subsection, we describe a goodness-of-fit test based on a decomposition of the data into subsets. We will sometimes refer to a subset of the data as an “experiment”, although we could divide the data into subsets in different ways.

Label the subset of the data that we wish to consider by an index $E$. Let $D(E)$ refer to the data in subset $E$, that is, all data points $D_{i}$ for $i\in E$. Recall Eq. (58), giving the probability to find data $D$ if the theory $T(a)$ is correct. The analogue of this that gives the probability to find data $D(E)$ if the theory $T(a)$ is correct is

[TABLE]

where

[TABLE]

$d\mu(D(E))$ is the measure

[TABLE]

$C_{E}$ is the matrix $C_{ij}$ for $i,j\in E$ and $N_{E}$ is the number of data in the subset $E$.

If we have a good fit, then the probability $P(D(E)|T(a))$ should be not too small for each subset $E$ of the data. That is, for each subset $E$, $\chi^{2}(D(E),a)$ should not be too large.

As we already mentioned, this is a much stronger criterion than the hypothesis-testing criterion. An individual experiment $E$ may be very badly fit in a large global fit (have an unacceptably high $\chi^{2}(D(E),a)$), even while the total $\chi^{2}(D,a)$ may look reasonable. In this section we ask, how large is too large for $\chi^{2}(D(E),a)$? What should the distribution of this quantity be?

Consider structure of the covariance matrix $C_{ij}$, defined by Eq. (56). Suppose first that each experiment has independent systematic errors that are not shared among the experiments, so that each source $J$ of systematic error is associated with just one experiment $E_{J}$. Then $\beta_{iJ}=0$ unless $i\in E_{J}$. The covariance matrix is then block-diagonal: $C_{ij}=C^{-1}_{ij}=0$ unless $i\in E$ and $j\in E$ for the same experiment label $E$. Then the total $\chi^{2}(D,a)$ is a sum of contributions $\chi^{2}(D(E),a)$ from the separate experiments:

[TABLE]

It is, however, not necessary that errors for experiment $E$ are uncorrelated with the errors from other experiments. If, for some of the data sets $E$, the covariance matrix has elements $C_{ij}$ that are non-zero for $i\in E$ and $j\notin E$, then Eq. (152) will fail. However the probability to find data $D(E)$ if the theory $T(a)$ is correct is still given by Eq. (149).

We return to Eq. (150). There are $N_{E}$ data in the set $E$. Then, with the parameters $a$ fixed to ideal values, ($a=\bar{a}$, as in Sec. III.5), and flawless theory ($R_{k}=0$), we find that $T_{k}(\bar{a})=\langle D_{k}\rangle$, and the distribution of $\chi^{2}(D(E),\bar{a})$ is the standard $\chi^{2}$ distribution with $N_{E}$ degrees of freedom.

If, however, we use the best-fit parameters $a_{\textrm{fit}}$ that are constrained by experiment $E$ as well as the rest of experiments, we find that $\chi^{2}(D(E),a_{\textrm{fit}})$ is approximately equal to the standard $\chi^{2}$ distribution with $N_{E}$ degrees of freedom up to a subleading term that can be determined as follows.

From Eqs. (77), (78), and (59), we find that

[TABLE]

and

[TABLE]

What is the expectation value of this at $a=a_{\textrm{fit}}$? We can refer to Eqs. (99), (103), (55) and (101) for $a=a_{\textrm{fit}}$ to arrive at

[TABLE]

where

[TABLE]

The first term in Eq. (155) is just the number of data, $N_{E}$, in the data set. The third term is small if the theory is good ($R_{i}\approx 0$). The second term $-N_{\mathrm{P}}(E)$, requires some analysis. We interpret $N_{\mathrm{P}}(E)$ as the effective number of parameters constrained by data set $E$. To support this interpretation, we note from Eqs. (86) and (137) that if $N_{\mathrm{P}}(E)$ is summed over all data sets $E$, it gives the number $N_{\mathrm{P}}$ of parameters. Now, the number of parameters, $N_{\mathrm{P}}$ is much smaller than the total number of data, $N_{\mathrm{D}}$. This makes it plausible that $N_{\mathrm{P}}(E)$ is small compared to the number $N_{E}$ of data in data set $E$: $N_{E}\gg N_{\mathrm{P}}(E)$. Then the first term in Eq. (155) dominates.

IV.8 Test of $\chi^{2}$ from individual experiments

What can we do with $\chi^{2}(D(E),a_{\textrm{fit}})$? Its value is given by Eq. (154) with $a=a_{\textrm{fit}}$. When the best-fit parameters are close to the true ones ($a_{\textrm{fit}}\approx\bar{a}$), and theory is nearly perfect ($R_{k}\approx 0$), the $\chi^{2}(D(E),a_{\textrm{fit}})$ distribution reduces to the form

[TABLE]

that, as we already know, obeys the $\chi^{2}$ distribution with $N_{E}$ degrees of freedom. Now, $a_{\textrm{fit}}$ is not exactly $\bar{a}$, but is influenced by the data in data set $E$. Thus one might think that the number of degrees of freedom is $N_{E}-N_{\mathrm{P}}(E)$ where $N_{\mathrm{P}}(E)$, Eq. (156), is the effective number of fit parameters associated with data set $E$. However, the reasoning in the previous subsection suggests that $N_{\mathrm{P}}(E)\ll N_{E}$, so that the number of degrees of freedom is approximately just the number $N_{E}$ of data in data set $E$. This expectation can be tested by calculating $N_{\mathrm{P}}(E)$ in the case that $\chi^{2}(D(E),a_{\textrm{fit}})$ appears anomalously small.

We will now check if the distributions of the observed $\chi^{2}(D(E),a_{\textrm{fit}})$ values from the experiments $E$ in actual PDF fits are close to the ideal distributions.

When $N_{E}$ is large, say, $N_{E}\gtrsim 30$, the $\chi^{2}$ distribution with $N_{E}$ degrees of freedom approaches the Gaussian distribution with with mean $\langle\chi^{2}(D(E),a_{\textrm{fit}})\rangle\approx N_{E}$ and standard deviation $\langle(\chi^{2}(D(E),a_{\textrm{fit}})-N_{E})^{2}\rangle\approx 2N_{E}$, as we have seen in Sec. III.8. For $N_{E}\lesssim 30$, the non-Gaussian features are pronounced and $\chi^{2}$ distributions with different $N_{E}$ are not easily compared. Conveniently for our purpose, the variable

[TABLE]

fluctuates with a distribution that is quite accurately131313Other definitions Lewis (1988) of $S_{E}$ are more accurate but with the very simple form (157), the distribution function $\rho(S_{E})$ matches the Gaussian distribution $(2\pi)^{-1/2}\exp(-S_{E}^{2}/2)$ to within 0.04 for $N_{E}=5$ and to within 0.01 for $N_{E}=50$. an ${\cal N}(0,1)$ distribution Lai et al. (2010a); Fisher (1925), namely

[TABLE]

Note that the $S_{E}$ distribution is independent of $N_{E}$. The original $N_{E}$ dependence for the distribution of $\chi^{2}$ was absorbed into the definition of $S_{E}$.

To test the quality of the fit, we can plot a histogram of the $S_{E}$ values for all of the experiments (or data sets) $E$ contributing to the fit. The histogram should match the Gaussian distribution (158).

To see how this should work, we can generate $S_{E}$ values for a number of randomly generated pseudoexperiments. The number $N_{E}$ of data for each pseudoexperiment is chosen at random between 0 and 3000. For each pseudoexperiment, we generate a value of $\chi_{E}^{2}$ at random according to the standard $\chi^{2}$ distribution with $N_{E}$ degrees of freedom. Then we define $S_{E}$ for that pseudoexperiment by $S_{E}=\sqrt{2\chi^{2}_{E}}-\sqrt{2N_{E}-1}$. In the left-hand plot in Fig. 8, we show the resulting histogram of $S_{E}$ values obtained for 35 pseudoexperiments, along with the expected Gaussian distribution (158). In the right-hand plot in Fig. 8, we show the analogous histogram for 500 pseudoexperiments.

The histogram in Fig. 8 for 500 pseudoexperiments is evidently pretty close to the expected distribution (158). For the histogram for 35 pseudoexperiments, it is not so clear merely by eye. For a quantitative estimate, we can apply the standard Anderson-Darling test, described in Sec. IV.5, of the probability that the observed distribution matches ${\cal N}(0,1)$. For the left-hand histogram in Fig. 8, we find that $P_{\textrm{A-D}}=0.53$ and, for the right-hand histogram, we find that $P_{\textrm{A-D}}=0.44$. These values indicate that it is quite plausible that the $S_{E}$ values were generated from ${\cal N}(0,1)$, which, to a good approximation, they were.

Now, let us turn to the distributions of $S_{E}$ from recent NNLO global analyses shown in Fig. 9. The parameters $S_{E}$ in the histograms are computed using Eq. (157) from $\chi^{2}$ and $N_{E}$ values listed in Table 5 of Harland-Lang et al. (2015) for MMHT 2014 NNLO, Tables 2.1-2.3 and 3.1 of Ball et al. (2017) for NNPDF 3.1 and 3.0, and internal tables of $\chi^{2}$ for CT14HERA2 NNLO Hou et al. (2017a). It is obvious that variations in $S_{E}$ are broader than the standard normal distribution expected in an ideal fit to all experiments, both in the positive and negative directions. We can estimate the differences by the mean and standard deviation for each observed distribution. For CT14HERA2 and MMHT2014 fits, the means are close to zero, indicating that, while some experiments are not fit well, the other experiments are fit too well. On the other hand, for the NNPDF3.0 and NNPDF3.1 analyses, the observed mean is of order 0.8 to 1 – more experiments are not fitted well than fitted too well. For the four fits, the probability values for matching the expected ${\cal N}(0,1)$ distribution according to the Anderson-Darling test are

[TABLE]

In all four cases, it is very unlikely that the observed distribution came from the expected Gaussian distribution.

We emphasize that none of the four PDF fits described above is a good fit according to the $P_{\textrm{A-D}}$ values obtained by breaking the data into smaller data sets, even though each fit is acceptable according to its total $\chi^{2}$ value.

The expectation that the $S_{E}$ distribution should match an $\mathcal{N}(0,1)$ distribution is based in part on the assumption that the parameters $R_{k}$ representing imperfections in the theory are negligible in Eq. (155). The evident failure of the distributions in Fig. 9 to match $\mathcal{N}(0,1)$ may indicate that the theory is not precise enough to match very precise experiments. In fact, some elevated $S_{E}$ values are contributed by the most precise experiments, such as the combined HERA 1+2 DIS data Abramowicz et al. (2015) and some LHC measurements. These experiments test QCD at unprecedented (NNLO) precision and thus may reveal evidence for new dynamical mechanisms. For instance, $S_{E}\approx 5.5$ for HERA 1+2 DIS data can be reduced to $S_{E}\approx 3$ by including small-$x$ resummation in DIS or by evaluating NNLO DIS cross sections with an $x$-dependent factorization scale (see the discussion in Sec. I.6). Similarly, the description of HERA 1+2 DIS and fixed-target DIS data such as BCDMS Benvenuti et al. (1990) is improved in the NNPDF3.1 analysis as compared to NNPDF3.0 in part by introducing the “fitted charm”, an independent and possibly process-dependent nonperturbative function that has similarities to power-suppressed (“higher-twist”) terms in DIS. In Sec. II.3, we briefly reviewed the rationale for optionally including the “fitted charm” in some PDF fits and the current limitations to its theoretical understanding.

A different class of concerns arises when too many data sets have negative $S_{E}$, indicating that the the data sets are systematically described better than would be expected in a good fit in which the data fluctuations around the best-fit theory are random. When too many experiments have negative $S_{E}$, the reduced global $\chi^{2}$ may hide some problems, such as overestimated experimental uncertainties or overfitting the statistical fluctuations because of using too flexible theory. Sometimes, a very low global $\chi^{2}\ll N_{D}-N_{P}$ is taken as evidence for perfect theory that justifies aggressive estimates for PDF uncertainties based exclusively on the $\Delta\chi^{2}=1$ criterion. In fact, the resulting small PDF uncertainties would be wrong: not only it is improbable that the very low $\chi^{2}$ is caused by random fluctuations consistently with the experimental errors, the parameterization uncertainty also needs to be estimated by trying other PDF parameterization forms that render $\chi^{2}$ of up to about $N_{D}-N_{P}$.

IV.9 Test of consistency between experiments

We can carry this analysis further by asking whether different experiments imposed consistent constraints on PDF parameters $a$. To this end, we consider an observable $\sigma(a)$ that depends on the parton parameters, as in Sec. III.7. The observable could be a cross section, as suggested by the notation, or, extending the notion of “observable” a bit, it could be the value of the PDF $f_{a/p}(x,\mu^{2})$ for a particular flavor at a particular momentum fraction $x$ and scale $\mu$.

As in Sec. III.7, as long as we consider parameters $a$ that are not far from the best fit parameters $a_{\textrm{fit}}$, we can apply a linear approximation for the evaluation of $\sigma(a)$,

[TABLE]

with $\sigma^{\prime}_{\alpha}=\partial\sigma(a)/\partial a_{\alpha}$. Furthermore, if we define a special vector $e(\sigma)$ according to Eq. (116),

[TABLE]

then we found in Sec. III.7 that we can evaluate $\sigma(a)-\sigma(a_{\textrm{fit}})$ by setting

[TABLE]

in Eq. (160). Variations of $a-a_{\textrm{fit}}$ in orthogonal directions $\{e(\sigma)^{(2)},e(\sigma)^{(3)},\dots\}$ do not contribute to $\sigma(a)$. That is

[TABLE]

The result is independent of the parameters $t_{n}$. That is, $t$ directly measures $\sigma(a)$. Let us denote the corresponding cross section by $\sigma_{0}(t)$.

The parameters $a_{\textrm{fit}}$ correspond to the minimum of the global $\chi^{2}$, so that, according to Eq. (119),

[TABLE]

Furthermore, if we evaluate Eq. (119) at a general point $a_{\textrm{fit}}+t\,e(\sigma)+{\scriptstyle\sum_{n}}t_{n}e(\sigma)^{(n)}$, we get

[TABLE]

If we regard parameter points $a$ as distributed at random according to a probability density proportional to $\exp(-(\chi^{2}-\chi^{2}_{\textrm{min}})/2)$, then, to find the probability $\rho$ for $a$ to lie in a plane of constant $\sigma(a)=\sigma_{0}(t)$, we simply integrate over the other variables $t_{n}$:

[TABLE]

That is, $\chi^{2}(D,a_{\textrm{fit}}+t\,e(\sigma))$ gives both the probability for $a$ to lie at position $t$ along the line $a_{\textrm{fit}}+t\,e(\sigma)$ and the probability for $a$ to lie in the plane $\sigma(a)=\sigma_{0}(t)$ that intersects this line at position $t$.

It may be useful to note that one can find the direction of $e(\sigma)$ quite simply. Up to its normalization, $e(\sigma)$ is the vector from $a_{\textrm{fit}}$ to the point on the surface $\sigma(a)=\sigma_{0}(t)$ that minimizes $\chi^{2}$ on this surface. The standard Lagrange multiplier method Stump et al. (2001) produces this vector.

After this introduction, let us explore the role of a single experiment, $E$, in the fit. Consider $\chi^{2}(D,a)$ for $a$ that varies along the line $a=a_{\textrm{fit}}+t\,e(\sigma)$. The parameter $t$ labels distance along this line. We use one of three sets of data $D$: either all of the data, $D(\textrm{all})$, or all of the data except for the data from experiment $E$, $D({\textrm{no }E})$, or the data from experiment $E$ alone, $D(E)$. We are interested in how the function $\chi^{2}(D,a_{\textrm{fit}}+t\,e(\sigma))$ depends on $t$ when we make these different choices for what data set $D$ we use in computing $\chi^{2}$.141414The dependence of $\chi^{2}$ on the position of the parameters $a$ along the line $a_{\textrm{fit}}+t\,e(\sigma)$ near $a=a_{\textrm{fit}}$ is easily determined. The dependence on $a$ over the whole plane $\sigma(a)=\sigma_{0}$ would require knowing $\chi^{2}$ in the entire parameter space. This is not so meaningful if we use just a small subset of the data, $D(E)$.

We will ask two questions concerning the role of experiment $E$ in determining $t$.

The first question is “Does experiment $E$ exert a substantial pull on the value of $t$?” To answer this question, we ask what would happen if we omitted the data from experiment $E$ from the evaluation of $\chi^{2}$. Then the minimum value of $\chi^{2}\left(D({\textrm{no }E}),a_{\textrm{fit}}+t\,e(\sigma)\right)$ will occur at a value $t({\textrm{no }E})$ that will typically be different from the value $t({\textrm{all}})=0$ that we get using all of the data. The $1\sigma$ uncertainty in $t({\textrm{all}})$ is $\Delta t({\textrm{all}})=1$. If the difference between $t(\textrm{no }E)$ and $t(\textrm{all})=0$ is smaller than this uncertainty, then we may conclude that experiment $E$ does not exert a substantial pull on $t$. That is, for experiment $E$ to substantially pull $t$, we need

[TABLE]

where $f$ is a parameter we could pick, perhaps $f=1$.

This is illustrated in Fig. 10, which is based on the CT18 NNLO fit Hou et al. (2019). For this illustration, we choose a very conservative value of $f$, $f=0.5$.

In panel (a), we choose the gluon distribution at $x=0.01$ and $\mu=125\mbox{ GeV}$ as our observable $\sigma$. The heavy black curve is the difference of $\chi^{2}(D(\textrm{all}),a_{\textrm{fit}}+t\,e(\sigma))$ and its minimum value as a function of the parameter $t$. The corresponding values of $g(0.01,125\mbox{ GeV})$ are shown along the top of the plot. We also show curves for the differences of $\chi^{2}(D(\textrm{no }E),a_{\textrm{fit}}+t\,e(\sigma))$ and their minimum values for three choices of data sets $E$, labeled experiments B1, B2, and C. [The experiments are taken from an actual global fit.] We see that the minima of all of these curves lie in the range $-0.5<t<0.5$, indicating that none of these data sets exerts a substantial pull on $t$ in the sense of Eq. (167) with $f=0.5$.

In panel (b), we choose the gluon distribution at $x=0.3$ and $\mu=125\mbox{ GeV}$ as our observable $\sigma$. Again, the heavy black curve is constructed from $\chi^{2}(D(\textrm{all}),a_{\textrm{fit}}+t\,e(\sigma))$. We also show curves for $\chi^{2}(D(\textrm{no }E),a_{\textrm{fit}}+t\,e(\sigma))$ for two choices of data sets $E$: experiments B1 and B2. These correspond to two data sets obtained from the same experiment B for two different collision energies $\sqrt{s}$. We see that the minima of these curves lie outside the range $-0.5<t<0.5$, indicating that both of these data sets exerts a substantial pull on $t$ in the sense of Eq. (167) with $f=0.5$.

Now suppose that experiment $E$ does exert a substantial pull on $t$. Then we need to check whether the global fit solution, $t(\textrm{all})$, is consistent with what experiment $E$ says. Let us define

[TABLE]

According to experiment $E$ alone, the best fit $t(E)$ is obtained by minimizing $\Delta\chi^{2}_{E}(t)$. The $1\sigma$ uncertainty range for $t(E)$ is given by $|\Delta\chi^{2}_{E}(t)-\Delta\chi^{2}_{E}(t(E)))|=1$. Thus the result $t=0$ from the full fit is consistent with the result from experiment $E$ alone if

[TABLE]

where $n$ is a parameter that we could pick, perhaps $n=2$ for consistency within a 95% confidence interval.

This is illustrated in Fig. 11, again based on the CT18 NNLO fit and the observable $\sigma=g(0.3,(125\textrm{ GeV})^{2})$. In Fig. 10(b), we saw that two data sets, B1 and B2, make at least a marginal difference in the overall fit. Now we plot $\Delta\chi^{2}_{\textrm{Total}}(t)$ for the overall fit as a heavy black line and also $\Delta\chi^{2}_{E}(t)$ for $E=$ B1 and $E=$ B2: the data sets obtained in the same experiment B that was repeated at two collider energies $\sqrt{s}$. We also exhibit the $\Delta\chi^{2}_{E}(t)$ curves for three other data sets. We see that $\Delta\chi^{2}_{\mathrm{B1}}(t)$ is about 10 units higher at $t=0$ than it is at its minimum. Thus the consistency condition (169) with $n=2$ is violated for experiment B1. On the other hand, $\Delta\chi^{2}_{\mathrm{B2}}(t)$ is only about 3 units higher at $t=0$ than it is at its minimum. Thus the consistency condition (169) with $n=2$ is satisfied for experiment B2.

At this level of inconsistency for experiment B1, it is not credible that we are simply looking at statistical fluctuations. One simple but crude way to remove the inconsistency would be to increase the error estimates for the discrepant data set(s). To illustrate how this might work, in Fig. 12, we have refitted the PDFs by assuming increased quoted errors for the data set B1. Namely, we refitted after multiplying the quoted errors of B1 by a constant factor $\sqrt{2}$. That is, in Eq. (59), we multiply $C_{ij}$ for $i,j\in B1$ by a common factor $1/2$.

Fig. 12(a) shows differences of $\chi^{2}$ from their minimum values as functions of $t$ for the the full data set $D(\textrm{all})$ and then for $D(\textrm{no }E)$ with $E=\mathrm{B1}$ (with the rescaled errors) and B2, as in Fig. 10(b). In general, we would expect that increasing the estimated errors from certain data sets would change the position of the best fit for the observable $\sigma$ and increase the estimated error on the prediction for $\sigma$. In this case, for $g(0.3,(125\textrm{ GeV})^{2})$, neither the position of the minimum nor the estimated error changes by much.

In Fig. 12(a), we see that the minima of the two curves occur well within the region $-f<t<f$, even for $f=0.5$. Since the criterion (167) no longer indicates that these two data sets exert substantial pulls on $t$, we need not examine the criterion (169) for a discrepancy between a data set and the overall fit.

If we do examine the criterion (169), we obtain Fig. 12(b), where we show $\Delta\chi^{2}_{E}(t)$ for the the full data set $D(\textrm{all})$ and then for $D(E)$ with $E=\mathrm{B1}$ (rescaled errors) and B2. We see that $\Delta\chi^{2}_{E}(t)$ is less than four units higher at $t=0$ than it is at its minimum both for B1 and B2. These represent less than $2\,\sigma$ discrepancies, which are not nearly as alarming as the discrepancy for experiment B1 that we saw in Fig. 11.

In summary, this analysis gives us criteria for checking whether there is a problem associated with the data from experiment $E$ in determining $t$. There is a problem if experiment $E$ exerts a substantial pull on the value of $t$, Eq. (167), and if the fit based on just experiment $E$ is inconsistent with the global fit, Eq. (169). There is one set of criteria for each independent direction $e(\sigma)$ corresponding to an observable $\sigma$ and for each experiment $E$.

We explored how one can make the results from a data set $E$ more consistent with the rest of the data by simply rescaling the errors for this data set. This is a very crude method. We do not recommend using it for finding the best fit. In the following section, we explore a more subtle method.

A less precise alternative is to leave the disagreeing experiment(s) and best fit based on these experiment(s) unchanged, but increase the PDF uncertainty to reflect the incompatibility in the experimental constraints. This possibility is discussed in Sec. IV.12.

IV.10 A more conservative way to adjust the errors

Suppose that the experiment $E$ exerts a substantial pull on the fit value of the observable $\sigma$ according to the criterion (167), and $\chi^{2}$ for data set $E$ along the line $a=a_{\textrm{fit}}+t\,e(\sigma)$ is not consistent according to the criterion (169) with its value at $t=0$, the best-fit value of $t$ according to the fit to all data. Then it may be helpful to increase the experimental errors for experiment $E$. We stated that simply increasing the total error estimate for experiment $E$ by dividing $C_{ij}$ for $i,j\in E$ by a common factor is a rather crude strategy. A more focused strategy would be to add a systematic error of the form (52) such that it resolves the inconsistency of the constraints on our observable $\sigma$, but does not affect the fit in any other way. We parameterize this systematic error as

[TABLE]

where $\xi$ is a constant that we can adjust, and

[TABLE]

In the normalization factor, $H(E)_{\alpha\beta}$ is the Hessian matrix, as in Eq. (86), but only including the data from data set $E$:

[TABLE]

With this factor, $\bar{\beta}_{k}$ is independent of the normalization of the vector $e(\sigma)$. The vector $e(\sigma)$ does have a definite normalization (161), but in Eq. (170), only the direction of $e(\sigma)$ matters. A simple relation for the normalization factor $e(\sigma)_{\alpha}H(E)_{\alpha\beta}e(\sigma)_{\beta}$ is given below in Eq. (188).

Let us examine whether adding a new systematic error of this form can repair the incompatibility between data set E and the rest of the data while not much affecting the remaining fit.

From Eq. (56) with the systematic error (170) added, the covariance matrix becomes

[TABLE]

Thus

[TABLE]

Using $dC/d\xi^{2}=-C\,[dC^{-1}/d\xi^{2}]\,C$, this is

[TABLE]

It is straightforward to solve this differential equation to obtain

[TABLE]

Here $C_{ij}=C(0)_{ij}$ is the covariance matrix without the added systematic error. With the definition Eq. (171) of $\bar{\beta}_{k}$, we have

[TABLE]

so

[TABLE]

What is the effect on $\chi^{2}$ for experiment $E$ of adding this systematic error? Consider $\chi^{2}$ for experiment $E$, for parameters

[TABLE]

where the $a^{\textrm{fit}}_{\alpha}$ are the parameters from the global fit before adding the extra systematic error, and the vector $e$ could be the special vector $e(\sigma)$ for observable $\sigma$, but could also be any other vector in the many dimensional space of parameters. We have, from Eqs. (59) and (77),

[TABLE]

We can write this as

[TABLE]

where

[TABLE]

gives the contribution linear in $t$ and

[TABLE]

is the Hessian matrix including the added systematic error, but just for the data from experiment $E$. Following the notation from Eq. (172), we define

[TABLE]

Without the added systematic error, we have

[TABLE]

When we add the new systematic error, the result changes to

[TABLE]

Here we used Eq. (178) to separate from $\chi^{2}(D(E),a,0)$ the extra terms resulting from the new systematic error (proportional to $\xi^{2}$).

We will now show that the new systematic error reduces an apparent tension between data set $E$ and the rest of the data along the line associated with $\sigma$ (our observable of interest). It does not modify constraints in the other directions. In the last two terms of Eq. (IV.10), the numerators contain $[e(\sigma)_{\alpha}H(E)_{\alpha\beta}e_{\beta}]$, the inner product between the unit vector $e(\sigma)$ defining the direction associated with $\sigma$, and another (possibly orthogonal) unit vector $e$ that defines the line $a_{fit}+te$ along which we choose to scan $\chi^{2}(D(E),a_{\textrm{fit}}+te,\xi)$. When we scan along the direction $e=e(\sigma)$, the $\chi^{2}(D(E),a_{\textrm{fit}}+te,\xi)$ function changes as

[TABLE]

If we turn off the systematic error (set $\xi=0$) for a moment, from Eq. (187) we can numerically find the factor $e(\sigma)_{\alpha}H(E)_{\alpha\beta}e(\sigma)_{\beta}$ that appears here and in the definition (171):

[TABLE]

If we turn the systematic error back on by choosing $\xi\neq 0$, the second term makes $\chi^{2}$ at $t=0$ smaller as $\xi$ increases without affecting the shape of $\chi^{2}$ as a function of $t$. In the remaining terms, the coefficients of $t$ and $t^{2}$ are reduced by the same factor as $\xi$ increases. The shape of $\chi^{2}$ versus $t$ simply becomes shallower as desired, reducing the tension between data set $E$ and the rest of the data along the line $a=a_{\textrm{fit}}+te(\sigma)$.

On the other hand, in the directions that are orthogonal to the direction of variation of $\sigma$, the unit vectors $e$ obey

[TABLE]

If, for example, the parameter space is 26-dimensional, then there is a 25-dimensional vector space in which $e$ could lie and satisfy this condition. In this case, Eq. (IV.10) gives

[TABLE]

That is, $\chi^{2}$ for data set $E$ at $t=0$ is smaller when $\xi>0$, but the shape of $\chi^{2}$ as a function of $t$ is not changed at all. Thus the new systematic error is a conservative choice in that it alleviates the incompatibility problem while having a minimal effect on the rest of the fit.

Let us see how this prescription resolves the tension between experiment B1 and other experiments that we observed in Fig. 11. There we examined the parameter fit along a direction $e(\sigma)$ corresponding to $\sigma=g(0.3,(125\ \mathrm{GeV})^{2})$. We saw that the curve of

[TABLE]

for a certain experiment $E=\mathrm{B1}$ was not consistent with the choice of $t=0$ that minimizes the total $\chi^{2}$. To alleviate this problem, we add the new systematic error (170) for experiment B1 with $\xi$ defined by $1/(1+\xi^{2})\approx 4/10$. With this error added, the new $\Delta\chi^{2}_{\textrm{B1}}(t)$ for B1 now satisfies

[TABLE]

The new $\Delta\chi^{2}_{\textrm{B1}}(t)$ curve becomes flat enough so that consistency with the rest of the data along the line $a_{\textrm{fit}}+t\,e(\sigma)$ is no longer a problem.

We can now perform the global fit again, with the modified systematic error for experiment B1. Then the best fit parameters $a_{\mathrm{fit}}$ change. The direction vector $e(\sigma)$ corresponding to the observable $g(0.3,(125\ \mathrm{GeV})^{2}$ also changes. The new fit gives us a new plot analogous to Fig. 11 in which all of the $\chi^{2}$ curves have changed. The result is shown in Fig. 13. In the new fit, $\chi^{2}$ for experiment B1 is very flat, indicating that experiment B1 is now not significantly affecting the determination of $g(0.3,(125\ \mathrm{GeV})^{2})$. The best fit value of $g(0.3,(125\ \mathrm{GeV})^{2})$ has changed from 0.309 to 0.312. The estimated error on the fit value of $g(0.3,(125\ \mathrm{GeV})^{2})$, determined by the second derivative of the total $\chi^{2}$ curve with respect to $g(0.3,(125\ \mathrm{GeV})^{2})$, is about 2% larger.

IV.11 Summary of measures of goodness of fit

The statistical analysis that we have presented in Sec. III relies on several assumptions. One assumes that experimental systematic errors are adequately represented by a Gaussian distribution described by a covariance matrix given by the experiments. Typically, one also assumes that theory errors are small enough to be neglected. Additionally, one assumes that the theoretical prediction is, to a good approximation, a linear function of the PDF parameters in the region $a_{\textrm{fit}}\pm\delta a$ near the best fit parameter choice. Some of these assumptions could be wrong.

We have argued that we should not simply trust these assumptions, but should test them by using a strong set of goodness-of-fit criteria. Taken together, these tests are much more stringent than that obtained by simply noting the global $\chi^{2}$ value. If the fit passes all of these tests, we can have some confidence in the results and the errors on the results. If the fit does not pass all tests, then remediation is needed. We do not offer a fixed prescription, but we have pointed out some possibilities.

The choice that has usually been made is to leave the fit as it is but use larger error estimates on the final PDFs than those found with the parameter-fitting criterion that requires $\Delta\chi^{2}=4$ at the 95% probability level. Estimation of trustworthy PDF errors in such imperfect situation can be difficult and sometimes controversial. Before a new generation of PDFs is published, it may undergo many months of multifaceted PDF testing in order to establish the realistic estimates for PDF uncertainties. The increase over the nominal PDF errors that results from this procedure is often referred to as applying tolerance to the PDF uncertainty.

IV.12 Global and dynamic tolerance

Tolerance is relevant at the stage of determination of PDF uncertainties, after the best-fit PDF has been found. In its simplest realization, tolerance defines an allowed range for the variation

[TABLE]

of the parameters $a$ along a direction $e$, at a probability level $v$. Often, $e$ is one of the eigenvectors of $H$, $H_{\alpha\beta}e_{\beta}=he_{\alpha}$. However, any direction $e$ is a possible choice. We define the normalization of $e$ using the Hessian matrix, as in Eq. (117): $e_{\alpha}H_{\alpha\beta}e_{\beta}=1$. Then the dependence of $\chi^{2}$ on $t$ is given by Eq. (119):

[TABLE]

According to Eq. (111), if the experiments that determine $a_{\textrm{fit}}$ were repeated many times, then the component $t$ of $a_{\textrm{fit}}-\bar{a}$ in the direction $e$ would be distributed according to ${\cal N}(0,1)$. Thus, if we pick a probability $v$ and ask that

[TABLE]

with probability $v$, the limiting value $t_{\textrm{lim}}(v)$ is determined by

[TABLE]

where $p(\bar{t})$ is the Gaussian distribution ${\cal N}(0,1)$. Then $t_{\textrm{lim}}(0)\approx 0$, $t_{\textrm{lim}}(0.68)\approx 1$, $t_{\textrm{lim}}(0.8)\approx 1.3$, and $t_{\textrm{lim}}(0.95)\approx 2$.

When we believe that the above procedure misestimates the true uncertainty on $t$, we could try to find a better probability distribution $p(\bar{t})$ to use in Eq. (196). For example, we could use ${\cal N}(0,T)$ with $T^{2}>1$ as our $p(\bar{t})$,

[TABLE]

Here we use the same value $T$ for every direction vector $e$ Pumplin et al. (2001). The value $T^{2}$ in this case is referred to as the global tolerance. With $v=0.68$, the allowed variation of PDF parameters will be constrained to satisfy $-T<t<T$ along any vector direction with this prescription. With $v=0.95$, the allowed variation of PDF parameters will be constrained to satisfy $-4T<t<4T$.

The dynamic tolerance introduced by the MSTW group Martin et al. (2009) is determined by a similar consideration, by constructing $p(\bar{t})$ from $\chi^{2}$ distributions for individual experiments $E$. If $P_{N}(\chi^{2})$ is the $\chi^{2}$ distribution with $N$ degrees of freedom, we can define $\xi(N,v)$ by

[TABLE]

so that $\chi^{2}<\xi(N,v)$ with probability $v$. Note that we choose a one-sided limit here. Since, according to Eq. (157), $S=\sqrt{2\chi^{2}}-\sqrt{2N-1}$ closely obeys the ${\cal N}(0,1)$ distribution, we can relate $\xi(N,v)$ to $t_{\textrm{lim}}(v)$ to a good approximation:

[TABLE]

With this information, the “dynamic tolerance” prescription of Martin et al. (2009) assigns an allowed interval

[TABLE]

for some eigenvector direction in the following way.

We define $\chi^{2}(D(E),a)$ to be the part of $\chi^{2}$ coming from only the data in data set $E$, as in Sec. IV.9. We use $\chi^{2}(D(E),a_{\textrm{fit}}+te)$ to define limits $T_{\textrm{min}}(E)$ and $T_{\textrm{max}}(E)$ arising from data set $E$, as explained below. Then we set

[TABLE]

For every $E$, Martin et al. (2009) define the range $T_{\textrm{min}}(E)<t<T_{\textrm{max}}(E)$ by the criterion

[TABLE]

where $N_{E}$ is the number of data in data set $E$. To understand the result of applying this criterion, it is helpful to use some approximations.

First, using Eq. (199) gives

[TABLE]

Noting that $t_{\textrm{lim}}(2v-1)$ is of order 1, we see that for $N_{E}\gg 1$, this is

[TABLE]

Now we examine the left-hand side of Eq. (203). The $\chi^{2}(D(E),a_{\textrm{fit}}+te)$ is a quadratic function of $t$,

[TABLE]

with the coefficients $A_{1}(E)$ and $A_{2}(E)$ given in Eq. (185). We note that $0<A_{2}(E)<1$:

[TABLE]

The coefficient $A_{1}(E)$ could have either sign and could be large.

Inserting Eq. (205) into Eq. (204) gives

[TABLE]

The large terms, 1, here cancel exactly. This gives

[TABLE]

To understand this, we can estimate $\chi^{2}(D(E),a_{\textrm{fit}})$ by its expectation value, which, according to Eq. (137) is approximately $N_{E}$. This gives

[TABLE]

This gives upper and lower limits on $t$ for each experiment $E$. If, for example, we take $v=0.9$ then $t_{\textrm{lim}}(2v-1)=t_{\textrm{lim}}(0.8)\approx 1.3$. For simplicity, consider the case that $A_{1}(E)$ is small. If $A_{2}(E)$ is also small, then this inequality restricts $t$ only weakly. That is, $|T_{\textrm{min}}(E)|$ and $|T_{\textrm{max}}(E)|$ are large. We always have $A_{2}(E)<1$. If $A_{2}(E)$ is close to 1, then this inequality can provide a significant restriction on $t$. However, the restriction is only significant if $\sqrt{2N_{E}}$ is not too large. For data sets with many data, the restriction is always weak. Thus the most restrictive values of $|T_{\textrm{min}}(E)|$ and $|T_{\textrm{max}}(E)|$, and thus the overall values of $|T_{\textrm{min}}|$ and $|T_{\textrm{max}}|$, are likely to come from data sets in which $\sqrt{2N_{E}}$ is not too large, and $A_{2}(E)$ is not too small. For most experiments, the values of $|T_{\textrm{min}}(E)|$ or $|T_{\textrm{max}}(E)|$ tend to be substantially greater than 1.

We do not attempt to justify the definition (202) of the range for $t$ or its approximate version (209). We do note, however, that the factor $\sqrt{2N_{E}}$ in Eq. (209) is familiar: it is the standard deviation for the distribution of $\chi^{2}$, as in Eq. (138) for $N_{D}\gg N_{P}$.

We have described one approach to defining a tolerance factor. There are several other approaches. It is beyond our scope to explore these in detail.

V Parton distributions for heavy ions

The concepts discussed in this article can be applied to nuclear parton distribution functions (nPDFs), nonperturbative QCD functions that are increasingly employed to model the structure of heavy nuclei in high-energy scattering. The concept of collinear QCD factorization that is central for describing scattering of free hadrons is also relevant for the growing number of measurements in collisions of heavy nuclei. The experimental data available for constraining the nPDFs is still very limited in their span over $x$ and $\mu^{2}$. They are anticipated to grow quickly as the Large Hadron Collider and as especially the envisioned Electron-Ion Collider produce new results. We will review the key features of the nPDFs and will refer the reader to the original publications by nuclear PDF analysis groups (EPPS Eskola et al. (2017), nCTEQ Kovařík et al. (2016), DSSZ de Florian et al. (2012), HKN Hirai et al. (2007)), KA Khanpour and Atashbar Tehrani (2016) and NNPDF Abdul Khalek et al. (2019a)) for more details.

V.1 Universality of nuclear PDFs

As in the standard case of protons, the structure functions and cross-sections in collisions involving one or more nuclei are related to nPDFs via perturbative QCD factorization Collins et al. (1983, 1985); Bodwin (1985); Collins et al. (1988). Even though inclusive collinear factorization has been proven only in a few cases (e.g. deeply inelastic scattering on hadrons and lepton pair production in hadron collisions), it has formed the basis for the analysis of proton PDFs. Eqs. (35) or (36) express QCD factorization for free-nucleon QCD observables. We will now show how to extend these factorization formulas to high-energy scattering of heavy nuclei. We will denote a nucleus by $\mathtt{N}$ and its atomic number by $\mathtt{A}$.

We assume that, at sufficiently large $\sqrt{s}$, the QCD observables of interest are dominated by independent parton scatterings, in which only one parton per initial-state nucleus contributes to the hard scattering. The interactions between a parton and a nucleus can be classified in three categories Accardi et al. (2004); Qiu (2003) shown in Fig. 14. Fig. 14(c) illustrates jet quenching, the changes in a jet as it moves through a large nucleus. Fig. 14(b) illustrates two partons from the nucleus participating in the hard interaction. These effects contribute to the power suppressed corrections to the factorization formulas (35) and (36) for nuclei, but they can be important if the nucleus is big enough and the scattering is not too hard.

The additional interactions between the nucleus and the initial-state parton within the same nucleus, shown in Fig. 14(a), change the parton distributions of the nucleus and do not affect the hard scattering, thereby leaving intact the form of factorization given by Eqs. (35) and (36). The only change in the prescription replaces the free-proton PDFs by the nPDFs, which account for the additional initial-state effects and can be defined as

[TABLE]

This definition is analogous to the one in Eq. (19), but the proton matrix element of the number density operator is replaced by the nuclear one. This nPDF is defined with respect to the whole nucleus with +-momentum $P_{\mathtt{N}}^{+}$. Accordingly, the parton described by this parton distribution function carries the +-momentum $p^{+}=\xi_{\mathtt{N}}\,P_{\mathtt{N}}^{+}$. The momentum fraction $\xi_{\mathtt{N}}$ is defined as

[TABLE]

The modified factorization prescription for a cross-section in collisions of nuclei $\mathtt{N}_{1}$ and $\mathtt{N}_{2}$ can be written using the nPDFs (210) as

[TABLE]

Even in collisions of nuclei, hard processes such as production of muon pairs or sufficiently high-$p_{T}$ jets are dominated by the leading-power contributions in Eq. (212). Thus these processes can be well described using $f_{i/\mathtt{N}}(\xi_{N},\mu^{2})$, where the dependence on the scale is still governed by DGLAP equations. However, the environment of the nuclear collisions is much different from the free-nucleon collisions. For example, at small values of $\xi_{\mathtt{N}}$, the parton momenta, as viewed in the nucleus rest frame, are very small, so that the parton wave functions spread over the whole nucleus and beyond. Then “saturation” Mueller (1999); Golec-Biernat and Wusthoff (1999); Bartels et al. (2002); Hautmann and Soper (2007) or “shadowing” Armesto (2006) can substantially modify the nPDFs. The nPDFs can incorporate these and other initial-state nuclear effects, such as the “EMC” effect Aubert et al. (1983); Geesaman et al. (1995); Malace et al. (2014), and still be universal.

On the other hand, jets produced in a partonic scattering can be altered by their passage through nuclear matter, as in Fig. 14(c), unless the jet’s transverse momentum is very large. This “jet quenching” can affect jet cross sections beyond what is predicted by Eq. (212) Aad et al. (2010); Chatrchyan et al. (2011); Adam et al. (2015). Jet-quenched contributions will not factorize in the same way.

The nPDFs depend on the number of protons, $\mathtt{Z}$, and number of neutrons, $\mathtt{A}-\mathtt{Z}$. Highly nontrivial $\mathtt{A}$ dependence arises from strong interactions of partons inside the nucleus. There is also trivial $\mathtt{A}$ dependence that would be present even if the nucleons were free. Consider a very simple model, in which a nucleus with +-momentum $P_{\mathtt{N}}^{+}$ is just a collection of comoving independent protons and neutrons, in which each nucleon carries the same fraction $\xi_{p,n}=1/\mathtt{A}$ of the total momentum $P_{\mathtt{N}}^{+}$.

In this model, one could write the nPDF of the whole nucleus as

[TABLE]

where $f_{a/p}$ and $f_{a/n}$ are the parton distributions in the free proton and neutron, and $\xi_{\mathtt{A}}$ is the momentum fraction of the +-momentum of the parton with respect to the +-momentum of the nucleon. The +-momentum $p^{+}$ of parton $a$ is

[TABLE]

so that the momentum fractions $\xi_{\mathtt{A}}$ and $\xi_{\mathtt{N}}$ are related via

[TABLE]

We can use this relation to rewrite Eq. (213) in terms of momentum fraction $\xi_{\mathtt{N}}$ as

[TABLE]

In this model, $\xi_{\mathtt{N}}$ is constrained to be in the range $0\leq\xi_{\mathtt{N}}\leq 1/\mathtt{A}$, since free-nucleon PDFs vanish for $\xi_{\mathtt{A}}>1$, and each nucleon carries exactly the fraction $1/\mathtt{A}$ of the +-momentum of nucleus.

In reality, one nucleon can carry any fraction of the nucleus’ +-momentum, since the nucleons participate in Fermi motion relative to each other Bodek and Ritchie (1981); Saito and Uchiyama (1985). We still find it helpful to use the momentum fraction $\xi_{\mathtt{A}}\equiv A\xi_{\mathtt{N}}$. We now define it to be the fraction of the average +-momentum $P^{+}_{\mathtt{N}}/\mathtt{A}$ of a bound nucleon. The variable $\xi_{\mathtt{A}}$ now takes values in the interval $0\leq\xi_{\mathtt{A}}\leq\mathtt{A}$, with contributions at $1<\xi_{\mathtt{A}}<\mathtt{A}$ arising from in-nucleus motion.

The PDFs of bound nucleons in the nucleus do not coincide with the free-nucleon PDFs. However, if nuclear modifications are moderate, we can start from Eq. (213) to get a reasonable ansatz for the parameterizations of nuclear PDFs.

We define a nuclear PDF of an average nucleon in a nucleus with atomic number $\mathtt{A}$, denoted by $f_{a}^{\mathtt{A}}(\xi_{\mathtt{A}},\mu^{2})$. This nPDF has the form

[TABLE]

In Eq. (217), $f_{a/p}^{\mathtt{A}}(\xi_{\mathtt{A}},\mu^{2})$ and $f_{a/n}^{\mathtt{A}}(\xi_{\mathtt{A}},\mu^{2})$ are the PDFs in the bound proton and bound neutron. They are different from the free-nucleon PDFs $f_{a/p,n}(\xi_{\mathtt{A}},\mu^{2})$. They depend on the momentum fraction $\xi_{\mathtt{A}}$ defined above. We can relate the two types of nPDFs that we just discussed:

[TABLE]

Either the nPDFs $f_{a/\mathtt{N}}(\xi_{\mathtt{N}},\mu^{2})$ in the nucleus or the nPDFs $f_{a}^{\mathtt{A}}(\xi_{\mathtt{A}},\mu^{2})$ for an average nucleon are acceptable for use in QCD calculations. But, “trivial” $\mathtt{A}$ dependence makes it difficult to compare the nPDFs of the first kind, $f_{a/\mathtt{N}}(\xi_{\mathtt{N}},\mu^{2})$, for two different nuclei.

For example, consider the prominent feature of proton PDFs: the peaks of the up- and down-quark distributions at $\xi\approx 1/3$. Similar peaks are found in the respective nPDFs $f_{u/\mathtt{N}}(\xi_{\mathtt{N}})$ and $f_{d/\mathtt{N}}(\xi_{\mathtt{N}})$ at $\xi_{\mathtt{N}}\sim 1/(3\mathtt{A})$, i.e., the position of the peaks in these nPDFs depends on the nucleus. In addition, the respective valence-quark distributions are normalized by the sum rules in a nucleus-dependent way:

[TABLE]

In contrast, the nPDFs for an average nucleon not only take into account the trivial $A$ dependence, they also correctly incorporate the specific ratio of protons to neutrons. The nPDFs $f_{a/p}^{\mathtt{A}}(\xi_{\mathtt{A}},\mu^{2})$ of a bound proton satisfy the sum rules

[TABLE]

which are much like the sum rules for the free proton.

Experimental analyses of nuclear DIS account for the trivial $\mathtt{A}$ dependence by presenting the cross-sections or DIS structure functions not for the whole nucleus but rather per nucleon. Similarly, for collisions between two nuclei with atomic numbers $\mathtt{A}_{1}$ and $\mathtt{A}_{2}$, cross sections $\tilde{\sigma}[F]\equiv\sigma[F]/(\mathtt{A}_{1}\mathtt{A}_{2})$ per nucleon are usually quoted. The cross-section $\tilde{\sigma}[F]$ can be expressed using either type of nPDFs:

[TABLE]

One should note that even in the case of non-interacting nucleons, the cross-section per nucleon $\tilde{\sigma}[F]$ would differ from the cross-section in proton-proton collisions due to the different flavor decomposition.

To summarize, the trivial $\mathtt{A}$ dependence reflecting the sheer number of the nucleons can be captured by using the ansatz (217) for the nPDF $f_{a}^{\mathtt{A}}(\xi_{\mathtt{A}},\mu^{2})$ per average bound nucleon. On the right-hand side of Eq. (217), we introduced the PDFs $f_{a/p}^{\mathtt{A}}$ and $f_{a/n}^{\mathtt{A}}$ for bound protons and neutrons that acquire non-trivial $\mathtt{A}$ dependence from a combination of nuclear effects. Their parameterization at the input scale $\mu_{0}$ is discussed in the next section.

V.2 Parameterizing the $\mathtt{A}$-dependence

In principle one can extract the nPDFs $f_{a/p}^{\mathtt{A}}(\xi_{\mathtt{A}},\mu^{2})$ of a bound proton from experimental data for each nucleus separately, without constructing a comprehensive model for initial-state nuclear effects. The current nuclear scattering data, however, are insufficient to determine the complete set of nPDFs for any single nucleus. The dependence of nuclear effects on $\xi$ (also denoted as $x$), $\mathtt{A}$ and $\mathtt{Z}$ is assumed to be unknown from the first principles. Thus, it must be determined a global fit to experimental data. To assemble all scattering data taken on various nuclei within a common global analysis, a number of simplifying assumptions need to be made.

First, given that the nuclear modifications in the bound-proton PDFs $f_{a/p}^{\mathtt{A}}$ are expected to be small, it makes sense to use the free-proton PDFs $f_{a/p}$ as the baseline for the parameterization of $f_{a/p}^{\mathtt{A}}$.

Second, to use the available data, one makes an assumption that the bulk of the nuclear corrections depends only on $\mathtt{A}$, the total number of nucleons of either isospin.

Third, the current data are not sufficient to constrain the nPDFs for momentum fractions $\xi_{\mathtt{A}}>1$, so all nPDF analyses assume that $0<\xi_{\mathtt{A}}<1$.151515In the following text we will denote the momentum fraction $\xi_{\mathtt{A}}$ restricted to the interval $\xi_{\mathtt{A}}\in(0,1)$ for simplicity by $x$. However, future fits may be able to include the region $1<\xi_{\mathtt{A}}$.

Fourth, we need to decide how to introduce the $\mathtt{A}$ dependence in $f_{a/p}^{\mathtt{A}}$. In practice, one of two approaches is taken.

The first approach introduces nuclear correction factors $R_{a}(x,\mathtt{A})$ at the input scale $\mu^{2}_{0}$:

[TABLE]

In Eq. (222) $f_{a/p}(x,\mu_{0}^{2})$, the corresponding PDF for a free proton, is held fixed during any nPDF analysis. The PDF $f_{a/n}^{\mathtt{A}}(x,\mu_{0}^{2})$ of a bound neutron is related to $f_{a/p}^{\mathtt{A}}(x,\mu_{0}^{2})$ by charge symmetry. All free parameters associated with the nuclear modification are contained in $R_{a}$. For example, the EPPS16 analysis Eskola et al. (2017) uses the following piecewise expression:

[TABLE]

where $\alpha=10x_{a}$, and all parameters $a_{k}$, $b_{k}$ and $c_{k}$ implicitly depend on the atomic number $\mathtt{A}$ and the PDF flavor $a$. A similar approach that employs a nuclear correction factor is followed by HKN07 Hirai et al. (2007) and DSSZ de Florian et al. (2012). Each analysis uses a different proton baseline, cf. Hirai et al. (2007), de Florian et al. (2012), and Eskola et al. (2017).

The second approach presented in Kovařík et al. (2016) does not operate with the nuclear correction factors $R_{a}$. It rather parameterizes the whole nPDF $f_{a/p}^{A}(x,\mu_{0}^{2})$ with a flexible functional form used for the free-proton PDF $f_{a/p}(x,\mu_{0}^{2})$, but with $\mathtt{A}$-dependent free parameters. As an example, in the nCTEQ15 analysis, the explicit parameterization at the input scale is

[TABLE]

and, similarly to the underlying CTEQ6 parameterization Pumplin et al. (2002), the parton combination $\bar{d}/\bar{u}$ is given by a different form:

[TABLE]

All free parameters $c_{k}$ depend on the atomic number,

[TABLE]

The coefficient $c_{k,0}=c_{k}(\mathtt{A}=1)$ is the underlying proton coefficient. It is held constant during the nCTEQ analysis. In the case of the nCTEQ15 analysis, the underlying proton coefficients $c_{k,0}$ are set to the coefficients from the proton analysis performed in Owens et al. (2007).

Over the years, there have been other approaches to including nuclear effects in analyses of parton distribution functions. For example, the first next-to-leading order analysis of nuclear parton distribution functions performed in de Florian and Sassot (2004) used a convolution approach to relate the nuclear PDF of an average nucleon to the one of a free nucleon as

[TABLE]

More recently the analysis of proton PDFs using the neural networks was also extended to the analysis of nuclear PDFs Abdul Khalek et al. (2019a). In that case the dependence on the number of nucleons is introduced as an additional parameter to the neural network. We will not go into details of any of these approaches as they have not (yet) been used in a global analysis of all nuclear scattering data.

V.3 Comparisons of nuclear PDFs

The nPDFs $f_{a/p}^{\mathtt{A}}(x,\mathtt{A},\mu_{0}^{2})$ of the bound proton are determined from experimental data sets taken on many different nuclei. Most of the data are still coming from deeply inelastic scattering and are provided in the form of nuclear correction factors

[TABLE]

The more recent data from neutrino DIS are provided as double differential cross-sections $d^{2}\sigma/(dx\,dQ^{2})$. The collider data from Fermilab, RHIC and the LHC are also provided as differential cross-sections (per nucleon). The coverage of the relevant nuclear world data is nowhere close to that of the data available for free-nucleon PDFs. Many features of nPDFs are still poorly known, especially outside of the interval $0.01\lesssim\xi_{\mathtt{A}}\lesssim 0.5$. Most notably, no data constrain the nuclear gluon PDF at low momentum fractions.

The comparison of different nuclear PDF ensembles is a little trickier than comparing free-proton PDFs. The deficit of precise data introduces strong sensitivity to the prior and methodological assumptions, such as the kinematic cuts, nPDF parameterization form, or the choice of the baseline free-proton PDFs.

The methods introduced in Sec. IV can illustrate the differences between the various nPDF analyses. First, in Fig. 15 we show the distributions of $S_{E}$, defined in Sec. IV.8, from four recent NLO global nPDF analyses. As in the case of the proton analyses shown in Fig. 9, the distributions of $S_{E}$ for the nPDF analyses are broader than the standard normal distribution ${\cal N}(0,1)$ expected from an ideal fit. Looking at the means and standard deviations of the distributions of $S_{E}$ shown in Fig. 15, we see that, except for the HKN07 analysis, all means are negative, indicating that more experiments were fitted too well. This indicates possibly overestimated uncertainties in multiple nuclear experiments. The prior assumptions made in the HKN07 analysis do not allow for a good description of many Drell-Yan total cross-section measurements by E772 and E866 experiments at Fermilab. Consequently, the $S_{E}$ distribution for HKN07 has its mean shifted to the right, and it is wider. Some caution is needed when comparing the $S_{E}$ distributions between the analyses in detail. For example, one entry with high $S_{E}$ in the EPPS16 analysis is the double-differential neutrino DIS cross-section from the CHORUS collaboration. This experiment is not included in the nCTEQ15 and HKN07 analyses. In the DSSZ, it is included only in the form of the structure functions $F_{2}$.

We can quantify the observation that the $S_{E}$ distributions are far from the ideal ${\cal N}(0,1)$ distribution using the Anderson-Darling test. The probability values that the distributions for the four nPDF analyses were drawn from ${\cal N}(0,1)$ are

[TABLE]

With the possible exception of the DSSZ distribution, the Anderson-Darling test confirms that it is very unlikely that the distributions in question come from the expected Gaussian distribution. This is reminiscent of what we found in the proton case in Eq. (159); however, in three cases out of four, the nuclear data are fitted too well, rather than too poorly.

The momentum fraction dependence of nPDFs is often examined by plotting scale-dependent nuclear correction factors,

[TABLE]

where $f_{i/p}(x,\mu^{2})$ is the baseline free-proton PDF. In Fig. 16, we turn to a comparison of the EPPS16, nCTEQ15 and DSSZ nuclear PDFs presented as these nuclear-correction factors. We show $R_{i}(x,Q^{2},\mathtt{A})$ at $Q=10$ GeV for lead ($\mathtt{A}=208$), for which the nuclear effects are the largest. Broadly speaking, we can conclude that all three nPDF families are consistent with each other within the indicated uncertainties. Upon a closer inspection, we see that the central values of $R_{i}$ differ substantially among the three nPDF sets, for a large part due to the strong dependence on the abovementioned methodological assumptions, and most prominently due to the choice of the parameterization form. Furthermore, even though it is conventional to compare the ratios $R_{i}$ rather than nPDFs themselves, this quantity artificially introduces a dependence on the proton baseline. Much of the dependence on the baseline is absent when one compares the bound-proton PDFs $f_{a/p}^{\mathtt{A}}(x,\mathtt{A},\mu^{2})$ directly.

The other notable difference among the results in Fig. 16 is their strikingly different uncertainties. One source of the differences are the various definitions of the uncertainties. All nPDF analyses employ some version of the global tolerance criterion that is based solely on the global $\chi^{2}$, cf. Sec. IV.12. The DSSZ analysis uses the simplest version of the tolerance: their uncertainties correspond to varying the underlying parameters along the eigenvector directions (see Eq. (193)) by $t=T=\sqrt{30}$. Both nCTEQ15 and EPPS16 analyses first examine a version of the dynamical tolerance, as described in Sec. IV.12, to estimate proper global tolerances for their final nPDF uncertainties. They determine the limits $T^{i}_{\textrm{min}}$ and $T^{i}_{\textrm{max}}$ according to Eq. (201) using the probability $v=0.90$ for each eigenvector direction $e_{i}$. Then, a global tolerance is constructed by averaging the changes in $\chi^{2}$ over all eigenvector directions as

[TABLE]

The nCTEQ15 analysis has 16 free parameters ($N_{P}=16$) and generates the error PDFs in a standard manner for the global tolerance of $T^{2}=35$.

The EPPS16 analysis was the first to include the LHC data from proton-lead collisions. It uses 20 free parameters, their prescription given by Eq. (231) yields the global tolerance of $T^{2}=52$.

If all nuclear PDF analyses were to use the same nuclear data in a specific range of momentum fractions, and all analyses had a flexible parameterization form, the uncertainties would be very similar. At present, the compared nPDF analyses do not fit the same data. Furthermore, as the four nPDF analyses rely on the traditional minimization of global $\chi^{2}$, introducing more free nPDF parameters that can be constrained by the nuclear data would lead to unstable global fits. In Secs. IV.1 and IV.2, we showed how one can find the optimal number $N_{P}$ of free parameters needed to obtain a stable fit to a given set of hadronic data. For the current nPDF analyses, the optimal number of free parameters appears to be no more than 15-20. Adding new data, for example the LHC data that are included in the EPPS16 analysis, allows one to expand the constraints to a wider range of momentum fractions or new parton flavor combinations. With more LHC data expected in the near future, it will be possible to open up additional free parameters in the initial nPDF parameterizations, leading to a more realistic estimate of uncertainties on nuclear PDFs.

VI Conclusions

We have reviewed certain aspects of the fitting of collinear parton distribution functions (PDFs) to data. This is a very large field. We have concentrated on just a few areas that could be of interest for the readers who use the PDFs or are interested in the rich subject of the global QCD analysis.

First, we have described the basic definition of what parton distribution functions are, and how they relate to the description of data. We have also provided definitions and a brief description for parton distributions in nuclei instead of in just protons and neutrons.

Second, we have described the basic statistical treatment needed to fit the PDFs using what is often called the Hessian method. Our description is simplified compared to what is actually used in current PDF fits. Most importantly, we have assumed that, in a small enough neighborhood of the best fit, the theory predictions $T_{k}(a)$ are approximately linear functions of the parameters $a$. This is not exactly the case, but this approximation is reasonably good when the PDF uncertainties are small and allows us to derive results in a closed form. Working within this framework, we have explored the statistical reasoning behind the fitting procedure and have derived analytic expressions for the key results of a PDF fit, such as expectation values and uncertainties.

We have then provided a battery of tests to critically examine whether the statistical assumptions are consistent with certain statistical measures that result from the fit. With these tests, one can identify specific features or data subsets in the multidimensional QCD fits that may indicate discrepancies between experimental measurements and theoretical predictions. Without insisting on a concrete recipe, we present some ideas of what one can do in the case of inconsistency.

Acknowledgments

We thank the Kavli Institute for Theoretical Physics at University of California in Santa Barbara, supported by National Science Foundation grant PHY11-25915, for hospitality during the start of the work on this review. We also thank the Aspen Center for Physics, supported by National Science Foundation grant PHY-1607611, for its hospitality during a later stage of our work. We are grateful to CTEQ colleagues for insightful discussions. Work at UO was supported by the United States Department of Energy under grant DE-SC0011640. Work at SMU was supported by the U.S. Department of Energy under grant DE-SC0010129.

Appendix A Transformation for $\chi^{2}(D,a,\lambda)$

In this appendix, we relate the form (66) for $\chi^{2}(D,a,\lambda)$ to the form Eq. (68), in which it is apparent that the minimum of $\chi^{2}(D,a,\lambda)$ with respect to the variables $\lambda$ is $\chi^{2}(D,a)$. We begin with $\chi^{2}(D,a,\lambda)$ as given in Eq. (66),

[TABLE]

where $B_{IJ}$ was defined in Eq. (69),

[TABLE]

and where

[TABLE]

Completing the square in the variables $\lambda$ gives

[TABLE]

Define shifted variables $\lambda$,

[TABLE]

and the matrix

[TABLE]

This gives

[TABLE]

The matrix $\widetilde{C}_{ij}$ is, in fact, the covariance matrix $C_{ij}$. To prove this, use the definition (56) of $C^{-1}_{ij}$, calculate $\sum_{j}\widetilde{C}_{ij}C^{-1}_{jk}$, and simplify the product using $\sum_{j}\beta_{jJ}\beta_{jL}=B_{JL}-\delta_{JL}$. The calculation gives

[TABLE]

so that $\widetilde{C}_{ij}=C_{ij}$.

We arrive at the form of $\chi^{2}(D,a,\lambda)$ given in Eq. (68):

[TABLE]

It is clear that minimizing $\chi^{2}(D,a,\lambda)$ with respect to $\lambda$, which is equivalent to setting $\lambda_{I,J}^{\prime}=0$, leaves only the first term in Eq. (240), which is $\chi^{2}(D,a)$ according to Eq. (59).