9. Quantum Chromodynamics

Revised August 2023 by J. Huston (Michigan State U.), K. Rab-
bertz (KIT) and G. Zanderighi (MPP, Garching; München Tech.
U.).

<table>
  <tbody>
    <tr>
      <td>9.1</td>
      <td>Basics</td>
      <td>153</td>
    </tr>
    <tr>
      <td>9.1.1</td>
      <td>Running coupling</td>
      <td>153</td>
    </tr>
    <tr>
      <td>9.1.2</td>
      <td>Quark masses</td>
      <td>154</td>
    </tr>
    <tr>
      <td>9.2</td>
      <td>Structure of QCD predictions</td>
      <td>154</td>
    </tr>
    <tr>
      <td>9.2.1</td>
      <td>Fully inclusive cross sections</td>
      <td>154</td>
    </tr>
    <tr>
      <td>9.2.2</td>
      <td>Processes with initial-state hadrons</td>
      <td>155</td>
    </tr>
    <tr>
      <td>9.2.3</td>
      <td>Cross sections with phase-space restrictions</td>
      <td>157</td>
    </tr>
    <tr>
      <td>9.2.4</td>
      <td>Accuracy of predictions</td>
      <td>161</td>
    </tr>
    <tr>
      <td>9.3</td>
      <td>Experimental studies of QCD</td>
      <td>161</td>
    </tr>
    <tr>
      <td>9.3.1</td>
      <td>Hadronic final-state observables</td>
      <td>162</td>
    </tr>
    <tr>
      <td>9.3.2</td>
      <td>QCD measurements at colliders</td>
      <td>163</td>
    </tr>
    <tr>
      <td>9.4</td>
      <td>Determinations of the strong coupling constant</td>
      <td>167</td>
    </tr>
    <tr>
      <td>9.4.1</td>
      <td>Hadronic τ decays and low Q<sup>2</sup> continuum:</td>
      <td>168</td>
    </tr>
    <tr>
      <td>9.4.2</td>
      <td>Heavy quarkonia decays:</td>
      <td>168</td>
    </tr>
    <tr>
      <td>9.4.3</td>
      <td>PDF fits:</td>
      <td>169</td>
    </tr>
    <tr>
      <td>9.4.4</td>
      <td>Hadronic final states of e<sup>+</sup>e<sup>−</sup> annihilations:</td>
      <td>169</td>
    </tr>
    <tr>
      <td>9.4.5</td>
      <td>Observables from hadron-induced collisions:</td>
      <td>169</td>
    </tr>
    <tr>
      <td>9.4.6</td>
      <td>Electroweak precision fit:</td>
      <td>171</td>
    </tr>
    <tr>
      <td>9.4.7</td>
      <td>Lattice QCD:</td>
      <td>171</td>
    </tr>
    <tr>
      <td>9.4.8</td>
      <td>Determination of the world average value of α<sub>s</sub>(m<sub>Z</sub><sup>2</sup>):</td>
      <td>172</td>
    </tr>
  </tbody>
</table>

9.1 Basics

Quantum Chromodynamics (QCD), the gauge field theory that describes the strong interactions of colored quarks and gluons, is the SU(3) component of the SU(3)×SU(2)×U(1) Standard Model of Particle Physics. The Lagrangian of QCD is given by

$$
\mathcal{L} = \sum_q \bar{\psi}_{q,a} (i\gamma^\mu \partial_\mu \delta_{ab} - g_s \gamma^\mu t_{ab}^C A_\mu^C - m_q \delta_{ab}) \psi_{q,b} - \frac{1}{4} F_{\mu\nu}^A F^{A\mu\nu}, \quad (9.1)
$$

where repeated indices are summed over. The $\gamma^\mu$ are the Dirac
$\gamma$-matrices. The $\psi_{q,a}$ are quark-field spinors for a quark of flavor
$q$ and mass $m_q$, with a color-index $a$ that runs from $a = 1$ to
$N_c = 3$, i.e. quarks come in three “colors.” Quarks are said to be
in the fundamental representation of the SU(3) color group.

The $A_\mu^C$ correspond to the gluon fields, with $C$ running from 1 to $N_c^2-1=8$, i.e. there are eight kinds of gluon. Gluons transform under the adjoint representation of the SU(3) color group. The $t_{ab}^C$ correspond to eight 3 × 3 matrices and are the generators of the SU(3) group (cf. the section on “SU(3) isoscalar factors and representation matrices” in this Review, with $t_{ab}^C \equiv \lambda_{ab}^C/2$). They encode the fact that a gluon’s interaction with a quark rotates the quark’s color in SU(3) space. The quantity $g_s$ (or $\alpha_s = g_s^2/(4\pi)$) is the QCD coupling constant. Besides quark masses, which have electroweak origin, it is the only fundamental parameter of QCD. Finally, the field tensor $F_{\mu\nu}^A$ is given by

$$
F_{\mu\nu}^{A} = \partial_{\mu} A_{\nu}^{A} - \partial_{\nu} A_{\mu}^{A} - g_s f_{ABC} A_{\mu}^{B} A_{\nu}^{C}, \\
[t^A, t^B] = i f_{ABC} t^C, \tag{9.2}
$$

where the $f_{ABC}$ are the structure constants of the SU(3) group.

Neither quarks nor gluons are observed as free particles. Hadrons are color-singlet (i.e. color-neutral) combinations of quarks, anti-quarks, and gluons.

Ab-initio predictive methods for QCD include lattice gauge the-
ory and perturbative expansions in the coupling. The Feynman
rules of QCD involve a quark-antiquark-gluon (q ¯ qg) vertex, a 3-
gluon vertex (both proportional to g$_{s}$ ), and a 4-gluon vertex (pro-
portional to g 2 $_{s}$). A full set of Feynman rules is to be found for
example in Refs. [1,2].

Adopting a standard notation where repeated indices are
summed over, useful color-algebra relations include: t<sub>ab</sub><sup>A</sup>t<sub>bc</sub><sup>A</sup> =
C<sub>F</sub>δ<sub>ac</sub>, where C<sub>F</sub> ≡ (N<sub>c</sub><sup>2</sup> − 1)/(2N<sub>c</sub>) = 4/3 is the color fac-
tor (“Casimir”) associated with gluon emission from a quark;

$f_{ACD}f_{BCD} = C_A\delta_{AB}$, where $C_A \equiv N_c = 3$ is the color factor associated with gluon emission from a gluon; $t_{ab}^A t_{ab}^B = T_R\delta_{AB}$, where $T_R = 1/2$ is the color factor for a gluon to split to a $q\bar{q}$ pair.

There is freedom for an additional *CP*-violating term to be present in the QCD Lagrangian, $\theta \frac{\alpha_s}{8\pi} F^A\mu\nu \tilde{F}^A_{\mu\nu}$, where $\theta$ is an additional free parameter, and $\tilde{F}_{\mu\nu}^A$ is the dual of the gluon field tensor, $\frac{1}{2}\epsilon_{\mu\nu\sigma\rho}F^{A\sigma\rho}$, with $\epsilon_{\mu\nu\sigma\rho}$ being the fully antisymmetric Levi-Civita symbol. Experimental limits on the electric dipole moment (EDM) of ultracold neutrons [3,4] and atomic mercury [5] constrain the QCD vacuum angle to satisfy $|\theta| \lesssim 10^{-10}$. For an overview of current and future experiments measuring EDMs see [6]. Further discussion is to be found in Ref. [7] and in the Axions section in the Listings of this Review.

This section will concentrate mainly on perturbative aspects of
QCD as they relate to collider physics. Related textbooks and
lecture notes include Refs. [1, 2, 8–11]. Aspects specific to Monte
Carlo event generators are reviewed in the dedicated section 43.
Lattice QCD is also reviewed in a section of its own, Sec. 17,
with further discussion of perturbative and non-perturbative as-
pects to be found in the sections on “Quark Masses”, “The
CKM quark-mixing matrix”, “Structure functions”, “Fragmenta-
tion Functions”, “High Energy Soft QCD and Diffraction”, “Pas-
sage of Particles Through Matter” and “Heavy-Quark and Soft-
Collinear Effective Theory” in this Review.

9.1.1 Running coupling

In the framework of perturbative QCD (pQCD), predictions for
observables are expressed in terms of the renormalized coupling
$\alpha_s(\mu_R^2)$, a function of an (unphysical) renormalization scale $\mu_R$.
When one takes $\mu_R$ close to the scale of the momentum transfer $Q$
in a given process, then $\alpha_s(\mu_R^2) \simeq Q^2$) is indicative of the effective
strength of the strong interaction in that process.

The coupling satisfies the following renormalization group equa-
tion (RGE):

$$
\mu_R^2 \frac{d\alpha_s}{d\mu_R^2} = \beta(\alpha_s) = -(b_0\alpha_s^2 + b_1\alpha_s^3 + b_2\alpha_s^4 + \dots), \quad (9.3)
$$

where $b_0 = (11C_A - 4n_f T_R)/(12\pi) = (33 - 2n_f)/(12\pi)$ is referred
to as the 1-loop $\beta$-function coefficient, the 2-loop coefficient is $b_1 =$
$(17C_A^2 - n_f T_R(10C_A + 6C_F))/(24\pi^2) = (153 - 19n_f)/(24\pi^2)$, and
the 3-loop coefficient is $b_2 = (2857 - \frac{5033}{9}n_f + \frac{325}{27}n_f^2)/(128\pi^3)$ for
the SU(3) values of $C_A$ and $C_F$. Here $n_f$ is the number of quark
flavors. The 4-loop coefficient, $b_3$, is to be found in Refs. [12, 13],
while the 5-loop coefficient, $b_4$, is in Refs. [14–18]. The coefficients
$b_2$ and $b_3$ (and beyond) are renormalization-scheme-dependent
and given here in the modified minimal subtraction scheme (MS)
[19], by far the most widely used scheme in QCD and the one
adopted in the following.

The minus sign in Eq. (9.3) is the origin of “asymptotic free-
dom” [20, 21], i.e. the fact that the strong coupling becomes
weak for processes involving large momentum transfers (“hard
processes”). For momentum transfers in the 0.1–1 TeV range,
$\alpha_s \sim 0.1$, while the theory is strongly interacting for scales around
and below 1 GeV.

The $\beta$-function coefficients, the $b_i$, are given for the coupling of an *effective theory* in which $n_f$ of the quark flavors are considered light ($m_q \ll \mu_R$), and in which the remaining heavier quark flavors decouple from the theory. One may relate the coupling for the theory with $n_f + 1$ light flavors to that with $n_f$ flavors through an equation of the form

$$
\alpha_s^{(n_f+1)}(\mu_R^2) = \alpha_s^{(n_f)}(\mu_R^2) \times \\
\qquad \times \left( 1 + \sum_{n=1}^{\infty} \sum_{\ell=0}^{n} c_{n\ell} [\alpha_s^{(n_f)}(\mu_R^2)]^n \ln^\ell \frac{\mu_R^2}{m_h^2} \right), \quad (9.4)
$$

where $m_h$ is the mass of the $(n_f+1)^{th}$ flavor, and the first few
$c_{n\ell}$ coefficients are $c_{11} = \frac{1}{6\pi}$, $c_{10} = 0$, $c_{22} = c_{11}^2$, $c_{21} = \frac{11}{24\pi^2}$,
and $c_{20} = -\frac{11}{72\pi^2}$ when $m_h$ is the $\overline{\text{MS}}$ mass at scale $m_h$, while
$c_{20} = \frac{7}{24\pi^2}$ when $m_h$ is the pole mass (mass definitions are discussed below in Sec. (9.1.2) and in the review on "Quark Masses"). Terms up to $c_{4\ell}$ with $0 \le \ell \le 4$ are to be found in Refs. [22, 23]. Numerically, when one chooses $\mu_R = m_h$, the matching is a modest effect, owing to the zero value for the $c_{10}$ coefficient. Relations between $n_f$ and $(n_f+2)$ flavors where the two heavy flavors are close in mass are given to three loops in Ref. [24].

Working in an energy range where the number of flavors is taken constant, a simple exact analytic solution exists for Eq. (9.3) only if one neglects all but the $b_0$ term, giving $\alpha_s(\mu_R^2) = (b_0 \ln(\mu_R^2/\Lambda^2))^{-1}$. Here $\Lambda$ is a constant of integration, which corresponds to the scale where the perturbatively-defined coupling would diverge. Its value is indicative of the energy range where non-perturbative dynamics dominates. A convenient approximate analytic solution to the RGE that includes terms up to $b_4$ is given by solving iteratively Eq. (9.3)

$$ \alpha_s(\mu_R^2) \approx \frac{1}{b_0 t} \left( 1 - \frac{b_1 \ell}{b_0^2 t} + \frac{b_1^2 (\ell^2 - \ell - 1) + b_0 b_2}{b_0^4 t^2} + \right. \\ \left. \qquad + \frac{b_1^3 (-2\ell^3 + 5\ell^2 + 4\ell - 1) - 6b_0 b_2 b_1 \ell + b_0^2 b_3}{2b_0^6 t^3} + \right. \\ \left. \qquad + \frac{18b_0 b_2 b_1^2 (2\ell^2 - \ell - 1) + b_1^4 (6\ell^4 - 26\ell^3 - 9\ell^2 + 24\ell + 7)}{6b_0^8 t^4} \right. \\ \left. \qquad + \frac{-b_0^2 b_3 b_1 (12\ell + 1) + 2b_0^2 (5b_2 + b_0 b_4)}{6b_0^8 t^4} \right), \quad (9.5) $$

with $t \equiv \ln \frac{\mu_R^2}{\Lambda^2}$ and $\ell = \ln t$, again parameterized in terms of a constant $\Lambda$. Note that Eq. (9.5) is one of several possible approximate 4-loop solutions for $\alpha_s(\mu_R^2)$, and that a value for $\Lambda$ only defines $\alpha_s(\mu_R^2)$ once one knows which particular approximation is being used. An alternative to the use of formulas such as Eq. (9.5) is to solve the RGE exactly, numerically (including the discontinuities, Eq. (9.4), at flavor thresholds). In such cases the quantity $\Lambda$ does not directly arise (though it can be defined, cf. Eqs. (1-3) of Ref. [25]). For these reasons, in determinations of the coupling, it has become standard practice to quote the value of $\alpha_s$ at a given scale (typically the mass of the Z boson, $m_Z$) rather than to quote a value for $\Lambda$.

A discussion of determinations of the coupling and a graph illustrating its scale dependence (“running”) are to be found in Section 9.4.

### 9.1.2 Quark masses

Free quarks have never been observed, which is understood as a result of a long-distance, confining property of the strong QCD force: up, down, strange, charm, and bottom quarks all hadronize, i.e. become part of a meson or baryon, on a timescale $\sim 1/\Lambda$; the top quark instead decays before it has time to hadronize. This means that the question of what one means by the quark mass is a complex one, which requires one to adopt a specific prescription. A perturbatively defined prescription is the pole mass, $m_q$, which corresponds to the position of the divergence of the quark propagator. This is close to one’s physical picture of mass. However, when relating it to observable quantities, it suffers from a badly behaved perturbative series which makes it ambiguous to an amount related to $\Lambda$ (see e.g. Ref. [26–28]). An alternative is the $\overline{\text{MS}}$ mass, $\overline{m}_q(\mu_R^2)$, which depends on the renormalization scale $\mu_R$.

Results for the masses of heavier quarks are often quoted either as the pole mass or as the $\overline{\text{MS}}$ mass evaluated at a scale equal to the mass, $\overline{m}_q(\overline{m}_q^2)$; light quark masses are often quoted in the $\overline{\text{MS}}$ scheme at a scale $\mu_R \sim 2$ GeV. The pole and $\overline{\text{MS}}$ masses are related by a series that starts as $m_q = \overline{m}_q(\overline{m}_q^2)(1 + \frac{4\alpha_s(\overline{m}_q^2)}{3\pi}) + \mathcal{O}(\alpha_s^2)$), while the scale-dependence of $\overline{\text{MS}}$ masses is given at lowest order by

$$ \mu_R^2 \frac{d\overline{m}_q(\mu_R^2)}{d\mu_R^2} = \left[ -\frac{\alpha_s(\mu_R^2)}{\pi} + \mathcal{O}(\alpha_s^2) \right] \overline{m}_q(\mu_R^2). \quad (9.6) $$

A more detailed discussion is to be found in a dedicated section of the Review, “Quark Masses”, with detailed formulas also in Ref. [29] and references therein.

In perturbative QCD calculations of scattering processes, it is customary to employ an approximation where the masses of all quarks, whose magnitudes are much smaller than the momentum transfer involved in the process, are neglected or set to zero.

## 9.2 Structure of QCD predictions

### 9.2.1 Fully inclusive cross sections

The simplest observables in perturbative QCD are those that do not depend on the initial-state hadrons and are fully inclusive, disregarding specific details of the final state. An example of such an observable is the total cross section for the process $e^+e^- \to$ hadrons at a center-of-mass energy $Q$. It can be expressed as follows:

$$ \frac{\sigma(e^+e^- \to \text{hadrons}, Q)}{\sigma(e^+e^- \to \mu^+\mu^-, Q)} \equiv R(Q) = R_{\text{EW}}(Q)(1 + \delta_{\text{QCD}}(Q)), \quad (9.7) $$

where $R_{\text{EW}}(Q)$ is the purely electroweak prediction for the ratio and $\delta_{\text{QCD}}(Q)$ is the correction due to QCD effects.

For the sake of simplicity, we can focus on energies where $Q \ll m_Z$, where the primary contribution to the process comes from photon exchange (while disregarding electroweak and finite-quark-mass corrections) denoted as $R_{\text{EW}} = N_c \sum_q e_q^2$, with $e_q$ representing the electric charges of the quarks. The QCD correction reads

$$ \delta_{\text{QCD}}(Q) = \sum_{n=1}^{\infty} c_n \cdot \left( \frac{\alpha_s(Q^2)}{\pi} \right)^n + \mathcal{O}\left( \frac{\Lambda^4}{Q^4} \right). \quad (9.8) $$

The first four terms in the $\alpha_s$ series expansion are [30]

$$ c_1 = 1, \qquad c_2 = 1.9857 - 0.1152n_f, \quad (9.9a) $$

$$ c_3 = -6.63694 - 1.20013n_f - 0.00518n_f^2 - 1.240\eta, \quad (9.9b) $$

$$ c_4 = -156.61 + 18.775n_f - 0.7974n_f^2 + 0.0215n_f^3 - (17.828 - 0.575n_f)\eta, \quad (9.9c) $$

with $\eta = (\sum e_q)^2 / (N_c \sum e_q^2)$ and $N_c = 3$. For corresponding expressions including also Z exchange and finite-quark-mass effects, see Refs. [31–33].

A related series holds also for the QCD corrections to the hadronic decay width of the $\tau$ lepton, which essentially involves an integral of $R(Q)$ over the allowed range of invariant masses of the hadronic part of the $\tau$ decay (see e.g. Ref. [34]). The series expansions for QCD corrections to Higgs-boson hadronic (partial) decay widths in the limit of heavy top quark and massless light flavors at $N^4\text{LO}$ are given in Ref. [35].

One characteristic feature of Eqs. (9.8) and (9.9) is that the coefficients of $\alpha_s^n$ increase order by order: calculations in perturbative QCD tend to converge more slowly than would be expected based just on the size of $\alpha_s$. The situation is significantly worse near thresholds or in the presence of tight kinematic cuts. Another feature is the existence of an extra “power-correction” term $\mathcal{O}(\Lambda^4/Q^4)$ in Eq. (9.8), which accounts for contributions that are fundamentally non-perturbative. All high-energy QCD predictions involve power corrections $(\Lambda/Q)^p$, although typically the suppression of these corrections with $Q$ is smaller than given in Eq. (9.8), where $p = 4$. The exact power $p$ depends on the observable and, for many processes and observables, it is possible to introduce an operator product expansion and associate power-suppressed terms with specific higher-dimension (non-perturbative) operators [36].

**Scale dependence.** In Eq. (9.8) the renormalization scale for $\alpha_s$ has been chosen equal to $Q$. The result can also be expressed in terms of the coupling at an arbitrary renormalization scale $\mu_R$,

$$ \delta_{\text{QCD}}(Q) = \sum_{n=1}^{\infty} \bar{c}_n \left( \frac{\mu_R^2}{Q^2} \right) \cdot \left( \frac{\alpha_s(\mu_R^2)}{\pi} \right)^n + \mathcal{O}\left( \frac{\Lambda^4}{Q^4} \right), \quad (9.10) $$
where $\bar{c}_1(\mu_R^2/Q^2) \equiv c_1$, $\bar{c}_2(\mu_R^2/Q^2) = c_2 + \pi b_0 c_1 \ln(\mu_R^2/Q^2)$, $\bar{c}_3(\mu_R^2/Q^2) = c_3 + (2b_0 c_2 \pi + b_1 c_1 \pi^2) \times \ln(\mu_R^2/Q^2) + b_0^2 c_1 \pi^2 \ln^2(\mu_R^2/Q^2)$, etc. Given an infinite number of terms in the $\alpha_s$ expansion, the $\mu_R$ dependence of the $\bar{c}_n(\mu_R^2/Q^2)$ coefficients will exactly cancel that of $\alpha_s(\mu_R^2)$, and the final result will be independent of the choice of $\mu_R$: physical observables do not depend on unphysical scales.<fnref n="1" />

With just terms up to some finite $n = N$, a residual $\mu_R$ dependence will remain, which implies an uncertainty on the prediction of $R(Q)$ due to the arbitrariness of the scale choice. This uncertainty will be $O(\alpha_s^{N+1})$, i.e. of the same order as the neglected higher-order terms. For this reason it is customary to use the scale dependence of QCD predictions as an estimate of the uncertainties due to neglected terms. One usually takes a central value for $\mu_R \sim Q$, in order to avoid the poor convergence of the perturbative series that results from the large $\ln^{n-1}(\mu_R^2/Q^2)$ terms in the $\bar{c}_n$ coefficients when $\mu_R \ll Q$ or $\mu_R \gg Q$. Uncertainties are then commonly determined by varying $\mu_R$ by a factor of two up and down around the central scale choice. This is not a rigorous prescription for determining the uncertainty, but is motivated by the requirement that there should not be large logarithms introduced into a calculation by large ratios of scales. A more detailed discussion on the accuracy of theoretical predictions and on ways to estimate the theoretical uncertainties can be found in Section 9.2.4.

## 9.2.2 Processes with initial-state hadrons

**Deep-Inelastic Scattering.** To illustrate the key features of QCD cross sections in processes with initial-state hadrons, let us consider deep-inelastic scattering (DIS), $e+p \to e+X$, where an electron $e$ with four-momentum $k$ emits a highly off-shell photon (momentum $q$) that interacts with the proton (momentum $p$). For photon virtualities $Q^2 \equiv -q^2$ far above the squared proton mass (but far below the squared $Z$ mass), the differential cross section in terms of the kinematic variables $Q^2$, $x = Q^2/(2p \cdot q)$ and $y = (q \cdot p)/(k \cdot p)$ is

$$ \frac{d^2\sigma}{dxdQ^2} = \frac{4\pi\alpha^2}{2xQ^4} \left[ (1+(1-y)^2)F_2(x, Q^2) - y^2F_L(x, Q^2) \right], \qquad (9.11) $$

where $\alpha$ is the electromagnetic coupling and $F_2(x, Q^2)$ and $F_L(x, Q^2)$ are proton structure functions, which encode the interaction between the photon and the proton. In the presence of parity-violating interactions (e.g. $\nu p$ scattering) an additional $F_3$ structure function is present. For an extended review, including equations for the full electroweak and polarized cases, see Sec. 18 of this Review.

Structure functions are not calculable in perturbative QCD, nor is any other cross section that involves QCD interactions and initial-state hadrons. To zeroth order in $\alpha_s$, the structure functions are given directly in terms of non-perturbative parton (quark or gluon) distribution functions (PDFs),

$$ F_2(x, Q^2) = x \sum_q e_q^2 f_{q/p}(x), \quad F_L(x, Q^2) = 0, \qquad (9.12) $$

where $f_{q/p}(x)$ is the non-perturbative PDF for quarks of type $q$ inside the proton, i.e. the number density of quarks of type $q$ inside a fast-moving proton that carry a fraction $x$ of its longitudinal momentum.

PDFs are primarily determined from data in global fits. PDF sets are available at different orders in perturbation theory (LO, NLO and NNLO) based on the order at which cross sections used in the fits are calculated. The evolution equations for the PDFs from one scale to another matches the perturbative accuracy of the cross sections. In modern global PDF fits, data are included from

DIS, Drell-Yan (DY), jets and $t\bar{t}$ processes, and more LHC collider data, with the global PDF fits using 3000-4000 data points. There is a large change in the PDFs from LO to NLO, with a much smaller change from NLO to NNLO. LO PDFs can be unreliable for collider predictions, especially at low and high $x$. The uncertainties for the resulting PDFs are determined primarily from the experimental uncertainties of the data that serve as input to the global PDF fits. There is also a component due to the limitations of the parameterizations used, as discussed e.g. in Ref. [37], although this may be reduced with the use of more flexible forms. The PDF uncertainties can either be determined through a Hessian approach or through the use of Monte Carlo replicas. It is now relatively straightforward to convert results from one approach to the other. There has been a recent combination (PDF4LHC21) of the CT18, MSHT20 and NNPDF3.1 PDFs, each produced from a global PDF fit, that serves to provide a more comprehensive view of both the central PDFs and of their uncertainties [38]. Recently, theoretical uncertainties related to missing higher orders have been included in global PDF determinations but so far only at NLO [39–42]. By tying renormalization scales for related processes, and factorization scales for all processes, together in the global fit, an estimate of the higher order uncertainties can be added to the PDF fit. An alternative approach [43] is to attempt to deal with cross sections rather than PDFs themselves.

Another approach is to use as much available information as possible from N³LO to examine the impact of such information on the resultant PDFs [44]. The result, of course, is sensitive to the higher order information that is still missing, for example from an incomplete knowledge of the splitting functions and the N³LO matrix elements.

Since the original work of Ref. [45], one of the intriguing questions is whether lattice calculations can provide valuable insights for extracting PDFs. Despite initial debates surrounding the underlying methodologies [46, 47], it is now understood that lattice data in Euclidean space can be connected to light-cone PDFs through factorization theorems, once lattice observables have been renormalized and extrapolated to the continuum limit [48]. Initial determinations based on lattice data have begun to emerge [48–61], and the comparisons with determinations based on conventional data show great promise. However, the uncertainties in lattice-based determinations are still quite large, therefore, for practical applications, PDFs are currently derived from global fits of experimental data (cf. Sec. 18 of this Review and also Refs. [62, 63]). It is likely that lattice simulations for light-cone PDFs will have an impact in the short-medium term if they focus on flavor components (strange and charm) and on the kinematic regions that are poorly covered by the experimental data [64]. Similarly, lattice determinations of transverse momentum dependent PDF (TMDs) and generalized parton distributions (GPDs) will be more relevant for phenomenology since the experimental data for these quantities are much less abundant and precise.

The result in Eq. (9.12), with PDFs $f_{q/p}(x)$ that are independent of the scale $Q$, corresponds to the “quark-parton model” picture, in which the photon interacts with point-like free quarks, or equivalently, one has incoherent elastic scattering between the electron and individual constituents of the proton. As a consequence, in this picture also $F_2$ and $F_L$ are independent of $Q$ [65]. When including higher orders in pQCD, one has

$$ F_2(x, Q^2) = x \sum_{n=0}^{\infty} \frac{\alpha_s^n (\mu_R^2)}{(2\pi)^n} \times \\ \times \sum_{i=q,g} \int_x^1 \frac{dz}{z} C_{2,i}^{(n)}(z, Q^2, \mu_R^2, \mu_F^2) f_{i/p}\left(\frac{x}{z}, \mu_F^2\right) + \mathcal{O}\left(\frac{\Lambda^2}{Q^2}\right). \qquad (9.13) $$

<footnote n="1">¹ With respect to pQCD there is an important caveat to this statement: at sufficiently high orders, perturbative series generally suffer from “renormalon” divergences $\alpha_s^n n!$ (reviewed in Ref. [26]). This phenomenon is not usually visible with the limited number of perturbative terms available today. However, it is closely connected with non-perturbative contributions and sets a limit on the possible precision of perturbative predictions. The cancellation of scale dependence will also ultimately be affected by this renormalon-induced breakdown of perturbation theory.</footnote>

<footnote>Just as in Eq. (9.10), we have a series in powers of $\alpha_s(\mu_R^2)$, each term involving a coefficient $C_{2,i}^{(n)}$ that can be calculated using Feynman diagrams. At variance with the parton model, the PDFs in pQCD depend on an additional scale, the factorization scale $\mu_F$, whose significance will be discussed in the following. Another important difference is the additional integral over $z$.</footnote>
The parton that comes from the proton can undergo a splitting before it enters the hard scattering. As a result, the $C_{2,i}^{(n)}$ coefficients are functions that depend on the ratio, $z$, of the parton's momentum before and after radiation, and one must integrate over that ratio. For the electromagnetic component of DIS with light quarks and gluons, the zeroth order coefficient functions are $C_{2,q}^{(0)} = e_q^2 \delta(1-z)$ and $C_{2,g}^{(0)} = 0$. Corrections are known up to $O(\alpha_s^3)$ (next-to-next-to-next-to-leading order, N³LO) for both electromagnetic [66] and weak currents [67, 68]. For heavy-quark production they are known to $O(\alpha_s^2)$ [69, 70] (next-to-leading order, NLO, insofar as the series starts at $O(\alpha_s)$). The current status of the theoretical description of unpolarized and polarized DIS is summarized in Ref. [71]. For precise comparisons of LHC cross sections with theoretical predictions, the photon PDF of the proton is also needed. It has been computed precisely in Ref. [72, 73] and has now been implemented in most global PDF fits. More recently, the PDF for leptons has also been computed [74]. The photon PDF has an elastic and an inelastic component, where the photon is generated from the intact proton or from the splitting of quarks inside the proton, respectively, while the lepton PDF is generated by photons splitting to collinear leptons.

The majority of the emissions that modify a parton's momentum are collinear (parallel) to that parton, and do not depend on whether the parton will interact with a photon. It is natural to view these emissions as modifying the proton's structure rather than being part of the coefficient function for the parton's interaction with the photon. Technically, one uses a procedure known as collinear factorization to give a well-defined meaning to this distinction, most commonly through the MS factorization scheme, defined in the context of dimensional regularization. The MS factorization scheme involves an arbitrary choice of factorization scale, $\mu_F$, whose meaning can be understood roughly as follows: emissions with transverse momenta above $\mu_F$ are included in the $C_{2,q}^{(n)}(z, Q^2, \mu_R^2, \mu_F^2)$; emissions with transverse momenta below $\mu_F$ are accounted for within the PDFs, $f_{i/p}(x, \mu_F^2)$. While collinear factorization is generally believed to be valid for suitable (sufficiently inclusive) observables in processes with hard scales, Ref. [75], which reviews the factorization proofs in detail, is cautious in the statements it makes about their exhaustivity, notably for the hadron-collider processes which we shall discuss below. Interesting considerations on the current status of our understanding of factorization can be found in Refs. [76, 77]. While the transverse degrees of freedom have been integrated over for collinear PDFs, it is also possible to produce transverse-momentum dependent (or unintegrated) PDFs where these degrees of freedom remain. For a recent comprehensive review of transverse-momentum-dependent parton distribution functions and fragmentation functions see Ref. [78].

For collinear PDFs, the resulting dependence on $\mu_F$ is described by the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations [79], which to leading order (LO) read<fnref n="2" />

$$ \mu_F^2 \frac{\partial f_{i/p}(x, \mu_F^2)}{\partial \mu_F^2} = \sum_j \frac{\alpha_s(\mu_F^2)}{2\pi} \int_x^1 \frac{dz}{z} P_{i \leftarrow j}^{(1)}(z) f_{j/p}\left(\frac{x}{z}, \mu_F^2\right), \quad (9.14) $$

with, for example, $P_{q \leftarrow g}^{(1)}(z) = T_R(z^2 + (1-z)^2)$. The other LO splitting functions are listed in Sec. 18 of this Review, while results up to NLO, $\alpha_s^2$, and NNLO, $\alpha_s^3$, are given in Refs. [80] and [81] respectively. At N³LO accuracy, only partial results are currently available in Refs. [82–86]. Splitting functions for PDFs in the helicity dependent case are given in Ref. [87].

Beyond LO, the coefficient functions are also $\mu_F$ dependent, for example $C_{2,i}^{(1)}(x, Q^2, \mu_R^2, \mu_F^2) = C_{2,i}^{(1)}(x, Q^2, \mu_R^2, Q^2) - \ln(\frac{\mu_F^2}{Q^2}) \sum_j \int_x^1 \frac{dz}{z} \times C_{2,j}^{(0)}(\frac{x}{z}) P_{j \leftarrow i}^{(1)}(z)$. Progress in higher-order QED and mixed QED-QCD corrections to the splitting functions can be found in refs. [88, 89].

As with the renormalization scale, the choice of factorization scale is arbitrary, but if one has an infinite number of terms in the perturbative series, the $\mu_F$-dependencies of the coefficient functions and PDFs will compensate each other fully. Given only N terms of the series, a residual $O(\alpha_s^{N+1})$ uncertainty is associated with the ambiguity in the choice of $\mu_F$. As with $\mu_R$, varying $\mu_F$ provides an input in estimating uncertainties on predictions. In inclusive DIS predictions, the default choice for the scales is usually $\mu_R = \mu_F = Q$.

As is the case for the running coupling, in DGLAP evolution one can introduce flavor thresholds near the heavy quark masses: below a given heavy quark's mass, that quark is not considered to be part of the proton's structure, while above it is considered to be part of the proton's structure and evolves with massless DGLAP splitting kernels. With appropriate parton distribution matching terms at threshold, such a variable flavor number scheme (VFNS), when used with massless coefficient functions, gives the full heavy-quark contributions at high $Q^2$ scales. For scales near the threshold, it is instead necessary to appropriately adapt the standard massive coefficient functions to account for the heavy-quark contribution already included in the PDFs [90–92].

At sufficiently small $x$ and $Q^2$ in inclusive DIS, resummation of small $x$ logarithms may be necessary [93, 94]. This may in fact have been observed in Refs. [95] based on HERA data [96], in a kinematic region where useful information for PDFs for collider predictions is present.

In situations, in which the center-of-mass energy $\sqrt{s}$ is much larger than all other momentum-transfer scales in the problem (e.g. $Q$ in DIS, $m_b$ for $b\bar{b}$ production in $pp$ collisions, etc.), each power of $\alpha_s$ beyond LO can be accompanied by a power of $\ln(s/Q^2)$ (or $\ln(s/m_b^2)$, etc.). This is variously referred to as the high-energy, small-$x$ or Balitsky-Fadin-Kuraev-Lipatov (BFKL) limit [94, 97, 98]. Currently it is possible to account for the dominant and first sub-dominant [99, 100] power of such logarithms at each order of $\alpha_s$, and also to estimate further sub-dominant contributions that are numerically large (see Refs. [101–104] and references therein). Progress towards NNLO is discussed in Ref. [105].

Physically, the summation of all orders in $\alpha_s$ can be understood as leading to a growth with $s$ of the gluon density in the proton. At sufficiently high energies this implies non-linear effects (commonly referred to as parton saturation), whose treatment has been the subject of intense study (see for example Refs. [106–108] and references thereto).

## Hadron-hadron collisions.

The extension to processes with two initial-state hadrons can be illustrated with the example of the total (inclusive) cross section for W boson production in collisions of hadrons $h_1$ and $h_2$, which can be written as

$$ \sigma(h_1 h_2 \rightarrow W + X) = \sum_{n=0}^{\infty} \alpha_s^n (\mu_R^2) \times \sum_{i,j} \int dx_1 dx_2 f_{i/h_1}(x_1, \mu_F^2) f_{j/h_2}(x_2, \mu_F^2) \times \hat{\sigma}_{ij}^{(n)} \rightarrow W + X \left(x_1 x_2 s, \mu_R^2, \mu_F^2\right) + \mathcal{O}\left(\frac{\Lambda^2}{m_W^4}\right), \quad (9.15) $$

where $s$ is the squared center-of-mass energy of the collision. At LO, $n=0$, the hard (partonic) cross section $\hat{\sigma}_{ij \rightarrow W+X}^{(0)}(x_1 x_2 s, \mu_R^2, \mu_F^2)$ is simply proportional to $\delta(x_1 x_2 s - m_W^2)$, in the narrow W-boson width approximation (see Sec. 51 of this Review for detailed expressions for this and other hard scattering cross sections). It is non-zero only for choices of $i, j$ that can directly give a $W$, such as $i=u, j=d$. At higher orders, $n \ge 1$, new partonic channels contribute, such as $gq$, and $x_1 x_2 s \ge m_W^2$ in the narrow W-boson width approximation.

Equation (9.15) involves a collinear factorization between the hard cross section and the PDFs, just like Eq. (9.13). As long as the same factorization scheme is used in DIS and $pp$ or $p\bar{p}$ (usually the MS scheme), then PDFs extracted in DIS can be directly used in $pp$ and $p\bar{p}$ predictions [75, 109] (with the anti-quark distributions

<footnote n="2">²LO is generally taken to mean the lowest order at which a quantity is non-zero.</footnote>
in an anti-proton being the same as the quark distributions in a proton).

A number of fully inclusive cross-sections have been computed up to N³LO, i.e. including corrections up to relative order $\alpha_s^3$, in recent years. These include inclusive Higgs production through gluon fusion in the large $m_t$ limit, calculated at N³LO in Refs. [110, 111]. Higgs boson pair production via gluon fusion in the same approximation was computed at N³LO in Ref. [112]. Calculations at this order, differential in the Higgs rapidity [113, 114] and fully differential [115] have been presented recently. Bottom-induced Higgs boson production has also been computed at N³LO [116]. Vector-boson fusion single- [117] and double-pair [118] production is also known to N³LO [117] in the factorized approximation. Neutral [119, 120] and charged [121] Drell-Yan processes, as well as associated Higgs production [116] have also been computed at N³LO. Differential rapidity distributions are available, but not yet with fiducial cuts. A number of public codes are now available which allow for phenomenological studies at N³LO [116, 117, 122–124]. The uncertainty band derived by varying the renormalization and factorization scales in the perturbative predictions is not always contained within the corresponding uncertainty band of the previous order, as might be expected from a well-behaved process. This has been ascribed to cancellations between channels at NNLO that lead to underestimates of the uncertainty at that order. This effect does not seem to be present in other processes, such as Higgs boson production through gluon-gluon fusion. Note that there is currently a mis-match in N³LO calculations in that PDFs at the same order are not yet available. Preliminary approximated N³LO PDFs [44] suggest that N³LO effects can be large, especially as far as heavy-quark distributions and low scales are concerned, as well as the gluon distribution in a region important for Higgs boson production through gluon-gluon fusion. The fully inclusive hard cross sections for several other processes are known to NNLO,<fnref n="3" /> for instance, Higgs-boson production in association with a vector boson [125], Higgs-pair production in the large $m_t$ approximation [126], and with full $m_t$ dependence [127], top-antitop production [128, 129], bottom-anti-bottom production [130] and vector-boson pair production [131–133]. For a comprehensive overview on other recent NNLO 2 → 2 calculations see Ref. [134, 135]. Other NNLO calculations with fiducial cuts will be discussed in Sec. 9.2.3.

### Photoproduction.

$\gamma p$ (and $\gamma \gamma$) collisions are similar to $p p$ collisions, with the subtlety that the photon can behave in two ways: there is “direct” photoproduction, in which the photon behaves as a point-like particle and takes part directly in the hard collision, with hard sub-processes such as $\gamma g \to q \bar{q}$; there is also resolved photoproduction, in which the photon behaves like a hadron, with non-perturbative partonic substructure and a corresponding PDF for its quark and gluon content, $f_{i/\gamma}(x, Q^2)$. While useful to understand the general structure of $\gamma p$ collisions, the distinction between direct and resolved photoproduction is not well defined beyond leading order, as discussed for example in Ref. [136].

## 9.2.3 Cross sections with phase-space restrictions

QCD final states always consist of hadrons, while perturbative QCD calculations deal with partons. Physically, an energetic parton fragments (“showers”) into many further partons, which then, on later timescales, undergo a transition to hadrons (“hadronization”). Fixed-order perturbation theory captures only a small part of these dynamics. This does not matter for the fully inclusive cross sections discussed above: the showering and hadronization stages are approximately unitary, i.e. they do not substantially change the overall probability of hard scattering, because they occur long after it has taken place (they introduce at most a correction proportional to a power of the ratio of timescales involved, i.e. a power of $\Lambda/Q$, where $Q$ is the hard scattering scale).

Less inclusive measurements, in contrast, may be affected by the extra dynamics. For those sensitive just to the main directions

of energy flow (jet rates, event shapes, cf. Sec. 9.3.1) fixed-order perturbation theory is often still adequate, because showering and hadronization do not substantially change the overall energy flow. This means that one can make a prediction using just a small number of partons, which should correspond well to a measurement of the same observable carried out on hadrons. For observables that instead depend on distributions of individual hadrons (which, e.g., are the inputs to detector simulations), it is mandatory to account for showering and hadronization. The range of predictive techniques available for QCD final states reflects this diversity of needs of different measurements.

While illustrating the different methods, we shall for simplicity mainly use expressions that hold for $e^+e^-$ scattering. The extension to cases with initial-state partons will be mostly straightforward (space constraints unfortunately prevent us from addressing diffraction and exclusive hadron-production processes; extensive discussion is to be found in Sec. 20 of this Review and in Refs. [137, 138]).

### 9.2.3.1 Soft and collinear limits

Before examining specific predictive methods, it is useful to be aware of a general property of QCD matrix elements in the soft and collinear limits. Consider a squared tree-level matrix element $|M_n^2(p_1, ..., p_n)|$ for the process $e^+e^- \to n$ partons with momenta $p_1, ..., p_n$, and a corresponding phase-space integration measure $d\Phi_n$. If particle $n$ is a gluon, which becomes collinear (parallel) to another particle $i$ and additionally its momentum tends to zero (is “soft”), the matrix element simplifies as follows,

$$ \lim_{\theta_{in} \to 0, E_n \to 0} d\Phi_n |M_n^2(p_1, \dots, p_n)| \\ = d\Phi_{n-1} |M_{n-1}^2(p_1, \dots, p_{n-1})| \frac{\alpha_s C_i}{\pi} \frac{d\theta_{in}^2}{\theta_{in}^2} \frac{dE_n}{E_n}, \quad (9.16) $$

where $C_i = C_F (C_A)$ if $i$ is a quark (gluon). This formula has non-integrable divergences both for the inter-parton angle $\theta_{in} \to 0$ and for the gluon energy $E_n \to 0$, which are mirrored also in the structure of divergences in loop diagrams. These divergences are important for at least two reasons: firstly, they govern the typical structure of events (inducing many emissions either with low energy or at small angle with respect to hard partons); secondly, they will determine which observables can be calculated within fixed-order perturbative QCD.

### 9.2.3.2 Fixed-order predictions

Let us consider an observable $O$ that is a function $O_n(p_1, ..., p_n)$ of the four-momenta of the $n$ final-state particles in an event (either partons or hadrons). In what follows, we shall consider the cross section for events weighted with the value of the observable, $\sigma_O$. As examples, if $O_n \equiv 1$ for all $n$, then $\sigma_O$ is just the total cross section; if $O_n \equiv \hat{\tau}(p_1, ..., p_n)$ where $\hat{\tau}$ is the value of the thrust for that event (see Eq. (9.22) in Sec. 9.3.1.2), then the average value of the thrust is $\langle \tau \rangle = \sigma_O/\sigma_{tot}$; if $O_n \equiv \delta(\tau - \hat{\tau}(p_1, ..., p_n))$ then one gets the differential cross section as a function of the thrust, $\sigma_O \equiv d\sigma/d\tau$.

In the expressions below, we shall omit write the non-perturbative power correction term, which for most common observables is proportional to a single power of $\Lambda/Q$.

**Leading Order.** If the observable $O$ is non-zero only for events with at least $n$ final-state particles, then the LO QCD prediction for the weighted cross section in $e^+e^-$ annihilation is

$$ \sigma_O^{LO} = \alpha_s^{n-2} (\mu_R^2) \int d\Phi_n |M_n^2(p_1, \dots, p_n)| O_n(p_1, \dots, p_n), \quad (9.17) $$

where the squared tree-level matrix element, $|M_n^2(p_1, ..., p_n)|$, including relevant symmetry factors, has been summed over all subprocesses (e.g. $e^+e^- \to q \bar{q} q \bar{q}$, $e^+e^- \to q \bar{q} g$) and has had all factors of $\alpha_s$ extracted in front. In processes other than $e^+e^-$ collisions, the center-of-mass energy of the LO process is generally not fixed, and so the powers of the coupling are often brought inside the integrals, with the scale $\mu_R$ chosen event by event, as a function of the event kinematics.

<footnote n="3">³Processes with jets or photons in the final state have divergent cross sections unless one places a cut on the jet or (dressed) photon momentum. Accordingly, they are discussed below in Section 9.2.3.2.</footnote>
Other than in the simplest cases (see the review on Cross Sections in this Review), the matrix elements in Eq. (9.17) are usually calculated automatically with programs such as CompHEP [139], MADGRAPH [140], ALPGEN [141], COMIX/SHERPA [142], and HELAC/PHEGAS [143]. Some of these (CompHEP, MADGRAPH) use formulas obtained from direct evaluations of Feynman diagrams. Others (ALPGEN, HELAC/PHEGAS and COMIX/SHERPA) use methods designed to be particularly efficient at high multiplicities, such as Berends-Giele recursion [144], which builds up amplitudes for complex processes from simpler ones (see also Refs. [145–147] for other tree-level calculational methods).

The phase-space integration is usually carried out by Monte Carlo sampling, because of the high dimensionality of the integration and in order to deal with the possibly involved kinematic cuts that are used in the corresponding experimental measurements. Perturbatively calculable observables should be insensitive to the emission of soft and collinear radiation. Because of the divergences in the matrix element, Eq. (9.16), at lowest order the integral converges only if the observable vanishes for kinematic configurations in which one of the $n$ particles is arbitrarily soft or it is collinear to another particle. As an example, the cross section for producing any configuration of $n$ partons will lead to an infinite integral, whereas a finite result will be obtained for the cross section for producing $n$ deposits of energy (or jets, see Sec. 9.3.1.1), each above some energy threshold and well separated from each other in angle.

At a practical level, LO calculations can be carried out for $2 \to n$ processes with $n \lesssim 6-10$. The exact upper limit depends on the process, the method used to evaluate the matrix elements (recursive methods are more efficient), and the extent to which the phase-space integration can be optimized to work around the large variations in the values of the matrix elements.

**NLO.** Given an observable that is non-zero starting from $n$ final-state particles, its prediction at NLO involves supplementing the LO result, Eq. (9.17), with the $2 \to (n+1)$-particle squared tree-level matrix element ($|M_{n+1}^2|$), and the interference of a $2 \to n$ tree-level and $2 \to n$ 1-loop amplitude (2Re($M_n M_{n,1-\text{loop}}^*$)),

$$ \begin{aligned} \sigma_{\mathcal{O}}^{\text{NLO}} = \sigma_{\mathcal{O}}^{\text{LO}} &+ \alpha_s^{n-1}(\mu_R^2) \int d\Phi_{n+1} |M_{n+1}^2(p_1, \dots, p_{n+1})| \\ & \times \mathcal{O}_{n+1}(p_1, \dots, p_{n+1}) + \alpha_s^{n-1}(\mu_R^2) \\ & \times \int d\Phi_n 2\text{Re}[M_n(p_1, \dots, p_n) M_{n,1-\text{loop}}^*(p_1, \dots, p_n)] \\ & \times \mathcal{O}_n(p_1, \dots, p_n). \end{aligned} \quad (9.18) $$

Relative to LO calculations, two important issues appear in the NLO calculations. Firstly, the extra complexity of loop-calculations relative to tree-level calculations means that automated calculations started to appear only about fifteen years ago (see below). Secondly, loop amplitudes are infinite in four dimensions, while tree-level amplitudes are finite, but their integrals over the phase space are infinite, due to the divergences of Eq. (9.16). These two sources of infinities have the same soft and collinear origins and cancel after the integration only if the observable $\mathcal{O}$ satisfies the property of infrared and collinear safety, which means that the observable is non-sensitive to soft emissions or to collinear splittings, i.e.

$$ \begin{aligned} \mathcal{O}_{n+1}(p_1, \dots, p_s, \dots, p_n) & \to \mathcal{O}_n(p_1, \dots, p_{s-1}, p_{s+1}, \dots, p_n) \\ & \text{if } p_s \to 0 \\ \mathcal{O}_{n+1}(p_1, \dots, p_a, p_b, \dots, p_n) & \to \mathcal{O}_n(p_1, \dots, p_a + p_b, \dots, p_n) \\ & \text{if } p_a \parallel p_b. \end{aligned} \quad (9.19) $$

Examples of infrared-safe quantities include event-shape distributions and jet cross sections (with appropriate jet algorithms, see below). Unsafe quantities include the distribution of the momentum of the hardest QCD particle (which is not conserved under collinear splitting), observables that require the complete absence of radiation in some region of phase space (e.g. rapidity gaps or 100% isolation cuts, which are affected by soft emissions), or the particle multiplicity (affected by both soft and collinear emissions). The non-cancellation of divergences at NLO due to infrared or collinear unsafety compromises the usefulness not only of the NLO calculation, but also that of a LO calculation, since LO is only an acceptable approximation if one can prove that higher-order terms are smaller. Infrared and collinear unsafety usually also implies large non-perturbative effects.

As with LO calculations, the phase-space integrals in Eq. (9.18) are usually carried out by Monte Carlo integration, so as to facilitate the study of arbitrary observables. Various methods exist to obtain numerically efficient cancellation among the different infinities. These include notably dipole [148], FKS [149] and antenna [150] subtraction.

Thanks to new ideas like the OPP method [151], generalized [152] and D-dimensional [153] unitarity, on-shell methods [154], and on-the-fly reduction algorithms [155], recent years have seen a breakthrough in the calculation of one-loop matrix elements (for reviews on unitarity based method see Ref. [156, 157]). Thanks to these innovative methods, automated tools for NLO calculations have been developed and a number of programs are available publicly: MADGRAPH5_aMC@NLO [140] and HELAC-NLO [158] provide full frameworks for NLO calculations; GoSAM [159], NJET [160], OPENLOOPS [161] and RECOLA [162] calculate just the 1-loop part and are typically interfaced with an external tool for a combination with the appropriate tree-level amplitudes. Other tools such as NLOJET++ [163], MCFM [164], VBFNLO [165], the PHOX family [166] or BLACKHAT [167] implement analytic calculations and provide full frameworks to compute NLO cross sections for selected classes of processes. Recently, a lot of attention has also been paid to the calculation of NLO electroweak corrections. Electroweak corrections are especially important for transverse momenta significantly above the W and Z masses, because they are enhanced by two powers of $\ln p_t/m_W$ for each power of the electroweak coupling, and close to Sudakov peaks$^4$, where most of the data lie and the best experimental precision can be achieved. In some cases the above programs can be used to calculate also NLO electroweak or beyond-standard-model corrections [168–174].

Given the progress in QCD and EW fixed-order computations, the largest unknown from fixed-order corrections is often given by the mixed QCD-electroweak corrections of $\mathcal{O}(\alpha_s\alpha)$. These mixed two-loop corrections are often available in an approximate form [175–188]. The first complete computation of the mixed QCD-EW corrections to the neutral-current Drell-Yan process appeared recently [189–191]. For a review on EW corrections to collider processes see Ref. [192].

**NNLO.** NNLO is considerably more complicated than NLO as it involves a further order in $\alpha_s$, consisting of: the squared $(n+2)$-parton tree-level amplitude, the interference of the $(n+1)$-parton tree-level and 1-loop amplitudes, the interference of the $n$-parton tree-level and 2-loop amplitudes, and the squared $n$-parton 1-loop amplitude.

Each of these elements involves large numbers of soft and collinear divergences, satisfying relations analogous to Eq. (9.16) which now involve multiple collinear or soft particles and higher loop orders (see e.g. Refs. [193–195]). Arranging for the cancellation of the divergences after numerical Monte Carlo integration has been one of the significant challenges of NNLO calculations, as has been the determination of the relevant 2-loop amplitudes. For the cancellations of divergences a wide range of methods has been developed. Some of them [196–205] retain the approach, inherent in NLO methods, of directly combining the separate loop and tree-level amplitudes. Others combine a suitably chosen, partially inclusive $2 \to n$ NNLO calculation with a fully differential $2 \to n+1$ NLO calculation [206–212]. The $q_T$-subtraction method was extended to deal with mixed QCD-QED corrections at NNLO [213]. For an overview of NNLO subtraction methods see Ref [214].

Quite a number of processes have been calculated differentially at NNLO so far. The state of the art for $e^+e^-$ collisions is $e^+e^- \to$

<footnote>$^4$For the definition of Sudakov form factors or Sudakov peaks see e.g. Refs. [1, 2, 11].</footnote>
3 jets [215–218].

For DIS, dijet production is known at NNLO [219, 220] and the description jet production has been recently pushed even to N³LO using the Projection-to-Born method [221, 222]. For hadron colliders, all $2 \to 1$ processes are known, specifically vector boson [223, 224] and Higgs boson production in the large $m_t$ limit [206, 225]. The finite top-mass corrections at this order have also been computed [226]. This calculation eliminates one important source of theoretical uncertainty to inclusive Higgs production at the LHC. Substantial progress has been made in the years for hadron-collider 2 → 2 processes, with calculations having been performed for nearly all relevant processes: ZZ [132] WW [131] and WZ [227], γγ [228, 229], Zγ [230] and Wγ [231] (many of these color singlet processes are available also in MCFM [232, 233] or MATRIX [133]), inclusive photon [234, 235], γ+jet [235, 236], W+jet [208, 237], Z+jet [236, 238, 239] H+jet [240–243], WH [244] and ZH [245], s-channel [246] and t-channel single-top [247–250], $t\bar{t}$ production [251, 252], dijet production [253, 254], W production in association with a c-jet [255, 256] and HH [126] (in large-top-mass approximation, see also the exact (two-loop) NLO result [127]). Recently, also the NNLO corrections to identified B-hadron hadro-production have been computed [257]. The frontier of NNLO calculation has now reached the complexity of $2 \to 3$ processes. The first $2 \to 3$ process known at NNLO has been Higgs production through vector-boson fusion, using an approximation in which the two underlying DIS-like $q \to qV$ scatterings are factorized, the so-called structure function approximation [211, 258]. Corrections beyond the structure function approximation are expected to be small, on the order of a percent or less [259]. More recently, first genuine $2 \to 3$ LHC processes have been described at NNLO accuracy, including three photon [260, 261], two photons and one jet [262], two jets and one photon [263], three-jets [264] and $Wb\bar{b}$ production [265]. The calculation of three-jet production might be relevant for future extractions of $\alpha_s(m_Z^2)$ from three-jet observables at the LHC [266].

Cross sections at the LHC are most often measured with fiducial cuts, for example on the transverse momenta and rapidities of the measured objects, restricting the measurement to regions where the objects have a good efficiency to be detected and are well-reconstructed. Ideally, the theoretical predictions should also be constructed at the same fiducial level; the other possibility is to extrapolate the experimental fiducial results to the full phase space. Such an extrapolation, however, requires the extrapolation (typically using a parton shower Monte Carlo) to be accurate over the full phase space.

Comparisons of fixed order predictions to fiducial measurements can sometimes result in the presence of large logarithms which degrade the accuracy of the prediction. Such is the case, for example, in the calculation of the Higgs rapidity distribution, in the diphoton final state. The imposition of transverse momenta and rapidity cuts on the two photons leads to an uncertainty notably greater at N³LO then at NNLO, due to the presence of these logs. One possible solution is to change the form of the kinematic cut on the photons, to a product of the two photon transverse momenta [267], reducing the impact of the logs; another solution is to perform a resummation of the logs, restoring the expected precision [268].

The Les Houches precision wishlist compiles predictions needed to fully exploit the data that will be taken at the High Luminosity LHC [269]. Most of the needed calculations require accuracy of at least NNLO QCD and NLO EW, and many require the prediction of $2 \to 3$ processes, such as $W/Z^+ \ge 2$ jets, $H^+ \ge 2$ jets, and $t\bar{t}H$ to NNLO. The latter process has been computed at NNLO in Ref. [270] in the approximation where the Higgs boson in the double-virtual contribution is soft. This calculation considerably reduces the theory uncertainty on this process.

As discussed in this section, calculations at NLO can now be relatively easily generated by non-experts using the programs mentioned above. By now there are also a number of publicly available tools to compute a range of processes at NNLO accuracy in QCD, such as EERAD3 [216] (for $e^+e^-$ collisions), [206], DYNNLO [224], DYTURBO [271], FEHIP [225], FEWZ [272], HTURBO [273], MATRIX [133], MCFM [233] or SUSHI [124, 274]. However, in some

cases NNLO calculations can be too complex and CPU-intensive to allow such an approach. In these cases, the relevant matrix element information for a specific observable can be stored by means of interpolation grids developed originally for cross sections at NLO [275–277] and recently extended to include electroweak corrections [278]. The application to NNLO has been demonstrated for $t\bar{t}$ and for DIS jet cross sections in Refs. [279, 280]. Each such interpolation grid corresponds to one fixed observable with specific selection criteria and binning, but is flexible with respect to the renormalization or factorization scale, the PDF, or the $\alpha_s$ evolution chosen for the cross-section computation. An even more flexible method, at the cost of requiring large amounts of storage space, saves huge numbers of partonic events in the form of ROOT n-tuples, which allow predictions to be generated on-the-fly for many observables of a particular process [281, 282]. Recent developments of both techniques for pp collisions at NNLO are described in Ref. [269].

### 9.2.3.3 Resummation

Many experimental measurements place tight constraints on emissions in the final state. For example, in $e^+e^-$ events, that (one minus) the thrust should be less than some value $\tau \ll 1$, or, in $pp \to Z$, events that the Z-boson transverse momentum or the transverse momentum of the accompanying jet should be much smaller than the Z-boson mass. A further example is the production of heavy particles or jets near threshold (so that little energy is left over for real emissions) in DIS and pp collisions.

In such cases, the constraint vetoes a significant part of the integral over the soft and collinear divergence of Eq. (9.16). As a result, there is only a partial cancellation between real emission terms (subject to the constraint) and loop (virtual) contributions (not subject to the constraint), causing each order of $\alpha_s$ to be accompanied by a large coefficient $\sim L^2$, where e.g. $L = \ln \tau$ or $L = \ln(m_Z/p_t^2)$. One ends up with a perturbative series, whose terms go as $\sim (\alpha_s L^2)^n$. It is not uncommon that $\alpha_s L^2 \gg 1$, so that the perturbative series converges very poorly if at all.<fnref n="5" /> In such cases one may carry out a “resummation”, which accounts for the dominant logarithmically enhanced terms to all orders in $\alpha_s$, by making use of known properties of matrix elements for multiple soft and collinear emissions, and of the all-orders properties of the divergent parts of virtual corrections, following original works such as Refs. [283–292] and also through soft-collinear effective theory [293, 294] (cf. also the section on “Heavy-Quark and Soft-Collinear Effective Theory” in this Review, as well as Ref. [295]).

For cases with double logarithmic enhancements (two powers of logarithm per power of $\alpha_s$), there are two classification schemes for resummation accuracy. Writing the cross section including the constraint as $\sigma(L)$ and the unconstrained (total) cross section as $\sigma_{tot}$, the series expansion takes the form

$$ \sigma(L) \approx \sigma_{tot} \sum_{n=0}^{\infty} \sum_{k=0}^{2n} R_{nk} \alpha_s^n (\mu_R^2) L^k, \quad L \gg 1, \quad (9.20) $$

and leading log (LL) resummation means that one accounts for all terms with $k = 2n$, next-to-leading-log (NLL) includes additionally all terms with $k = 2n - 1$, etc. Often $\sigma(L)$ (or its Fourier or Mellin transform) exponentiates,<fnref n="6" />

$$ \sigma(L) \approx \sigma_{tot} \exp \left[ \sum_{n=1}^{\infty} \sum_{k=0}^{n+1} G_{nk} \alpha_s^n (\mu_R^2) L^k \right], \quad L \gg 1, \quad (9.21) $$

where one notes the different upper limit on $k$ ($\le n+1$) compared to Eq. (9.20). This is a more powerful form of resummation: the

<footnote n="5">⁵To be precise one should be aware of two causes of the divergence of perturbative series. That which interests us here is associated with the presence of a new large parameter (e.g. ratio of scales). It is distinct from the “renormalon” induced factorial divergences of perturbation theory discussed above.</footnote>

<footnote n="6">⁶Whether or not this happens depends on the quantity being resummed. A classic example involves two-jet rate in $e^+e^-$ collisions as a function of a jet-resolution parameter $y_{cut}$. The logarithms of $1/y_{cut}$ exponentiate for the $k_t$ (Durham) jet algorithm [296], but not [297] for the JADE algorithm [298] (both are discussed below in Sec. 9.3.1.1).</footnote>
$G_{12}$ term alone reproduces the full LL series in Eq. (9.20). With the form Eq. (9.21) one still uses the nomenclature LL, but this now means that all terms with $k=n+1$ are included, and NLL implies all terms with $k=n$, etc.

For a large number of observables, NLL resummations are available in the sense of Eq. (9.21) (see Refs. [299–301] and references therein). NNLL has been achieved for the DY and Higgs-boson $p_t$ distributions [302–305] (also available in the CuTe [306], HRes [307] and ResBos [308] families of programs and also differentially in vector-boson decay products [309]) and related variables [310], for the $p_t$ of vector-boson pairs [311, 312], for the back-to-back energy-energy correlation in $e^+e^-$ [313], the jet broadening in $e^+e^-$ collisions [314], the jet-veto survival probability in Higgs and Z boson production in $pp$ collisions [315–317]$^7$ an event-shape type observable known as the beam thrust [318], hadron-collider jet masses in specific limits [319] (see also Ref. [320]), the production of top anti-top pairs near threshold [321–323] (and references therein), and high-$p_t$ $W$ and $Z$ production [324]. Automation of NNLL jet-veto resummations for different processes has been achieved in Ref. [325] (cf. also the NLL automation in Ref. [326]), while automation for a certain class of $e^+e^-$ observables has been achieved in Ref. [327]. $N^3LL$ resummations are available for the thrust variable, C-parameter and heavy-jet mass in $e^+e^-$ annihilations [328–330], for $p_t$ distribution of the Higgs boson [331, 332] and weak gauge bosons [333–335] and for Higgs- and vector-boson production near threshold [336]. In order to make better contact with experimental measurements, recent years have seen an increasing interest in resummations in exclusive phase-space regions and joint resummations [337–346]. Finally, there has also been considerable progress in resummed calculations for jet substructure, whose observables involve more complicated definitions than is the case for standard resummations, see e.g. Refs. [347–351] and references therein (see also Sec. 9.3.1.3).

### 9.2.3.4 Fragmentation functions

Since the parton-hadron transition is non-perturbative, it is not possible to perturbatively calculate quantities such as the energy-spectra of specific hadrons in high-energy collisions. However, one can factorize perturbative and non-perturbative contributions via the concept of fragmentation functions. These are the final-state analogue of the parton distribution functions which are used for initial-state hadrons. Like parton distribution functions, they depend on a (fragmentation) factorization scale and satisfy a DGLAP evolution equation.

It should be added that if one ignores the non-perturbative difficulties and just calculates the energy and angular spectrum of partons in perturbative QCD with some low cutoff scale $\sim \Lambda$ (using resummation to sum large logarithms of $\sqrt{s}/\Lambda$), then this reproduces many features of the corresponding hadron spectra [352]. This is often taken to suggest that hadronization is “local”, in this sense it mainly involves partons that are close both in position and in momentum.

Section 19 of this Review provides further information and references on these topics, including also the question of heavy-quark fragmentation.

### 9.2.3.5 Parton-shower Monte Carlo event generators

Parton-shower Monte Carlo (MC) event generators like PYTHIA [353–355], HERWIG [356–358], and SHERPA [359]$^8$ provide fully exclusive simulations of QCD events at the level of measurable particles, the so-called “particle level” or “hadron level”. Here, “measurable” refers to color-neutral particles with mean lifetimes $\tau$ long enough to be associated with tracks or decay vertices in particle detectors. Usually, this requires mean decay lengths $c\tau$ of around 10 mm. As such MC event generators are a crucial tool for all applications that involve simulating the response of detectors

to QCD events. Here we give only a brief outline of how they work and refer the reader to Sec. 43 and Ref. [361] for a full overview.

In general, we expect parton-shower matched predictions to differ from the underlying fixed-order results in regions where (1) there is a large sensitivity to jet shapes (for instance small R jets), (2) there is a restriction in phase space such that soft gluon resummation effects become important, (3) the observable contains multiple disparate scales, (4) there are perturbative instabilities at fixed order, e.g. related to kinematical cuts, and (5) the observable is sensitive to higher multiplicity states than those described by the fixed-order calculation (for an explicit study of some of these effects see e.g. [362]).

The MC generation of an event involves several stages. It starts with the random generation of the kinematics and partonic channels of whatever hard scattering process the user has requested at some high scale $Q_0$ (for complex processes, this may be carried out by an external program). This is followed by a parton shower, usually based on the successive random generation of gluon emissions (and $g \to q\bar{q}$ splittings). Leading contributions to the shower have emissions that are ordered according to some ordering variable. Common choices of scale for the ordering of emissions are virtuality, transverse momentum or angle. Each emission is generated at a scale lower than the previous emission, following a (soft and collinear resummed) perturbative QCD distribution, which depends on the momenta of all previous emissions. Parton showering stops at a scale of order 1 GeV, at which point a hadronization model is used to convert the resulting partons into hadrons. One widely-used model involves stretching a color “string” across quarks and gluons, and breaking it up into hadrons [363, 364]. Another breaks each gluon into a $q\bar{q}$ pair and then groups quarks and anti-quarks into colorless “clusters”, which then give the hadrons [356]. As both models are tuned primarily to LEP data, the cluster and string models provide similar results for most observables sensitive to hadronization [362]. For $pp$ and $\gamma p$ processes, modeling is also needed to treat the collision between the two hadron remnants, which generates an underlying event (UE), usually implemented via additional $2 \to 2$ scatterings (“multiple parton interactions”, MPI) at a scale of a few GeV, following Ref. [365]. The parameter values for the MPI models must be determined from fits to minimum-bias and/or underlying-event observables of LHC collision data. As the different MC event generators are adapted to essentially the same measurements, ideally the respective MPI implementations should lead to similar predictions for each program. One complication, however, is some non-universality of the underlying event among different physics processes.

A deficiency of the soft and collinear approximations that underlie parton showers is that they may fail to reproduce the full pattern of hard wide-angle emissions, important, for example, in many new physics searches. It is therefore common to use LO multi-parton matrix elements to generate hard high-multiplicity partonic configurations as additional starting points for the showering, supplemented with some prescription (CKKW [366], MLM [367]) for consistently merging samples with different initial multiplicities. Monte Carlo generators, as described above, compute cross sections for the requested hard process that are correct at LO.

A wide variety of processes are available in MC implementations that are correct to NLO, using the MC@NLO [368] or POWHEG [369] prescriptions, through the MADGRAPH5_aMC@NLO [140], POWHEG-BOX [370] and SHERPA [142] programs. Techniques have also been developed to combine NLO plus parton shower accuracy for different multiplicities of final-state jets [371]. While NLO+PS accurate predictions can now be implemented in the above tools rather easily and are available for a range of processes, given the advances in NNLO calculations it is natural to seek for NNLO+PS accurate Monte Carlo as well. The two main approaches today are the MiNNLO [372] and Geneva [373] ones. NNLO plus shower accuracy was achieved first for Drell-Yan and Higgs production [374, 375] and is now available for several $2 \to 2$ color singlet processes [376–390], as well as for heavy-quark pair-production [391–393].

It is important to understand/verify the accuracy of the par-

<footnote>$^7$A veto on the jet phase space can be severe, for example by requiring exactly zero jets above a given transverse momentum cut accompanying a Higgs boson, or relatively mild, for example by placing a transverse momentum cut of 30 GeV on the measurement of the production of a Higgs boson with one or more jets. In general, inclusive cross sections are preferable, as uncertainties on both the theoretical and experimental sides are smaller.</footnote>

<footnote>$^8$The program ARIADNE [360] has also been widely used for simulating $e^+e^-$ and DIS collisions.</footnote>
ton shower predictions in the Monte Carlo programs. A general framework for assessing the logarithmic accuracy of parton-shower algorithms has been formulated, based on their ability to reproduce the singularity structure of multi-parton matrix elements, and their ability to reproduce logarithmic resummation results [394]. The dominant contributions relevant for the extension of parton showers to higher logarithmic accuracy have been computed [395–399] and included in some algorithms.

There exist ways to improve on current parton-shower algorithms [400, 401] and to demonstrate the parton shower accuracy through a comparison to analytic next-to-leading logarithmic calculations for a range of observables [401]. Considerable progress has been achieved in recent years [400–412]. References [401, 404, 405, 408, 409] introduced the first shower implementations towards NLL, but these implementations did not include spin-correlation effects and lacked full color dependence at LL level. Ref. [402] incorporated full color effects at the NLL level, while spin correlations at NLL were included in Refs. [406, 407], resulting in an NLL accurate parton shower for proton-proton collisions [410, 411]. The modifications necessary for processes such as Deep Inelastic Scattering (DIS) and Vector Boson Fusion (VBF) are addressed in [412]. The next steps involve incorporating non-perturbative effects, accounting for quark masses, and tuning NLL parton shower on experimental data.

### 9.2.4 Accuracy of predictions

Estimating the accuracy of perturbative QCD predictions is not an exact science. It is often said that LO calculations are accurate to within a factor of two. This is based on experience with NLO corrections in the cases where these are available. In processes involving new partonic scattering channels at NLO and/or large ratios of scales (such as jet observables in processes with vector bosons, or the production of high-$p_t$ jets containing B-hadrons), the ratio of the NLO to LO predictions, commonly called the “K factor”, can be substantially larger than two. NLO corrections tend to be large for processes for which there is a great deal of color annihilation in the interaction. In addition, NLO corrections tend to decrease as more final state legs are added.

For calculations beyond LO, one approach to estimate the perturbative uncertainty is to base it on the last known perturbative correction; this may lead to misleading results if new sub-processes are present at the next-higher order. A more widely used method is to estimate it from the change in the prediction when varying the renormalization and factorization scales around a central value $Q$ that is taken close to the physical scale of the process. A conventional range of variation is $Q/2 < \mu_R, \mu_F < 2Q$, varying the two scales independently with the restriction $\frac{1}{2}\mu_R < \mu_F < 2\mu_R$ [413]. This constraint limits the risk of misleadingly small uncertainties due to fortuitous cancellations between the $\mu_F$ and $\mu_R$ dependence when both are varied together, while avoiding the appearance of large logarithms of $\mu_R^2/\mu_F^2$ when both are varied completely independently. Where possible, it can be instructive to examine the two-dimensional scale distributions ($\mu_R$ vs. $\mu_F$) to obtain a better understanding of the interplay between $\mu_R$ and $\mu_F$. This procedure should not be assumed to always estimate the full uncertainty from missing higher orders, but it does indicate the size of one important known source of higher-order ambiguity.<fnref n="9" />

Most $2 \to 2$ processes at the LHC are now known to NNLO. This typically results in a large reduction in the uncertainty from that obtained at NLO. However, care must be taken in the estimate of the uncertainty in processes containing jets, as accidental cancellations of the scale uncertainties may result an artificial reduction of the scale uncertainty; for some jet $R$ values, the uncertainty may even be zero. There are several possibilities for providing a more realistic estimate of the scale uncertainty for such processes [362], but none have been widely adopted at the LHC.

Calculations that involve resummations usually have an additional source of uncertainty associated with the choice of argument of the logarithms being resummed, e.g. $\ln(2\frac{p_t^Z}{m_Z})$ as opposed

to $\ln(\frac{1}{2}\frac{p_t^Z}{m_Z})$, as well as a prescription to switch off resummation effects when the logarithm is not large. In addition to varying renormalization and factorization scales, it is common practice to vary the argument of the logarithm by a suitable factor in either direction with respect to the default argument.

The accuracy of QCD predictions is limited also by non-perturbative (or hadronization) corrections, which typically scale as a power of $\Lambda/Q$.<fnref n="10" /> For measurements that are directly sensitive to the structure of the hadronic final state, the corrections are usually linear in $\Lambda/Q$. The non-perturbative corrections are further enhanced in processes with a significant underlying event (i.e. in $pp$ and $p\bar{p}$ collisions) and the impact of non-perturbative corrections is larger in cases where the perturbative cross sections fall steeply as a function of $p_t$ or some other kinematic variable, for example in inclusive jet $p_t$ or dijet mass spectra. Under high-luminosity running conditions, such as 13 TeV at the LHC, there can be on the order of 50 minimum-bias interactions occurring at each beam-beam crossing. This additional energy needs to be subtracted, and is typically removed by means of a rapidity-dependent transverse energy density determined on an event-by-event basis [419]. This subtraction, of necessity, also removes the underlying event, which must be added back in by means of MC event generator modelling if one wants to restore the measured event to the particle level.

Non-perturbative corrections are commonly estimated from the difference between Monte Carlo events at the “parton level” and at particle level. Parton level refers to the stage of the parton shower, where the evolution is stopped at an energy scale of typical hadron masses of a few GeV. An issue to be aware of is that “parton level” is not a uniquely defined concept. For example, in a MC event generator such a procedure depends on an arbitrary and tunable internal cutoff scale that separates the parton showering from the hadronization. In contrast, no such cutoff scale exists in an NLO or NNLO partonic calculation. The uncertainties in these corrections are often estimated by comparing different tunes of the various MC event generators. It should be noted that such estimates are not guaranteed to fully cover the true uncertainties.

Alternative methods exist for estimating hadronization corrections, that attempt to analytically deduce non-perturbative effects in one observable based on measurements of other observables (see the reviews [26, 420]). While they directly address the problem of different possible definitions of parton level, it should also be said that they are far less flexible than Monte Carlo programs and not always able to provide equally good descriptions of the data.

One of the main issues is whether the fixed-order partonic final state of a NLO or NNLO prediction can match the parton shower in its ability to describe the experimental jet shape (minus any underlying event). Calculations at NNLO provide a better match to the parton shower predictions than do NLO ones, as might be expected from the additional gluon available to describe the jet shape. The hadronization predictions appear to work for both orders, but at an unknown accuracy.

## 9.3 Experimental studies of QCD

Since we are not able to directly measure partons (quarks or gluons), but only hadrons and their decay products, a central issue for every experimental study of perturbative QCD is establishing a correspondence between observables obtained at the parton and the hadron level. The only theoretically sound correspondence is achieved by means of *infrared and collinear safe* quantities, which allow one to obtain finite predictions at any order of perturbative QCD.

As stated above, the simplest case of infrared- and collinear-safe observables are inclusive cross sections. More generally, when measuring fully inclusive observables, the final state is not analyzed at all regarding its (topological, kinematical) structure or its composition. Basically the relevant information consists in the rate of a process ending up in a partonic or hadronic final state. In $e^+e^-$ annihilation, widely used examples are the ratios of partial widths or branching ratios for the electroweak decay of particles into hadrons or leptons, such as $Z$ or $\tau$ decays, (cf. Sec. 9.2.1).

<footnote n="9">⁹Various studies have been carried out on how to estimate uncertainties from missing higher orders that go beyond scale variations [414–418].</footnote>

<footnote n="10">¹⁰In some circumstances, the scale in the denominator could be a smaller kinematic or physical scale that depends on the observable.</footnote>
Such ratios are often favored over absolute cross sections or par-
tial widths because of large cancellations of experimental and the-
oretical systematic uncertainties. The strong suppression of non-
perturbative effects, $O(\Lambda^4/Q^4)$, is one of the attractive features
of such observables, however, at the same time, the sensitivity to
radiative QCD corrections is small, which for example leads to a
larger statistical uncertainty when using them for the determina-
tion of the strong coupling constant. In the case of $\tau$ decays not
only the hadronic branching ratio is of interest, but also moments
of the spectral functions of hadronic $\tau$ decays, which sample dif-
ferent parts of the decay spectrum and thus provide additional
information. Other examples of fully inclusive observables are
structure functions (and related sum rules) in DIS. These are ex-
tensively discussed in Sec. 18 of this Review.
On the other hand, often the structure or composition of the
final state are analyzed, and cross sections differential in one or
more variables characterizing this structure are of interest. Ex-
amples are jet rates, jet substructure, event shapes or transverse
momentum distributions of jets or vector bosons in hadron colli-
sions. The case of fragmentation functions, i.e. the measurement
of hadron production as a function of the hadron momentum rel-
ative to some hard scattering scale, is discussed in Sec. 19 of this
Review.

On the other hand, often the structure or composition of the
final state are analyzed, and cross sections differential in one or
more variables characterizing this structure are of interest. Ex-
amples are jet rates, jet substructure, event shapes or transverse
momentum distributions of jets or vector bosons in hadron colli-
sions. The case of fragmentation functions, i.e. the measurement
of hadron production as a function of the hadron momentum rel-
ative to some hard scattering scale, is discussed in Sec. 19 of this
*Review*.

It is worth mentioning that, besides the correspondence between
the parton and hadron level, also a correspondence between the
hadron level and the actually measured quantities in the detector
has to be established. The simplest examples are corrections for
finite experimental acceptance and efficiencies. Whereas accep-
tance corrections essentially are of theoretical nature, since they
involve extrapolations from the measurable (partial) to the full
phase space, other corrections such as for efficiency, resolution and
response are of experimental nature. For example, measurements
of differential cross sections such as jet rates require corrections in
order to relate, e.g., the energy deposits in a calorimeter to the jets
at the hadron level. Typically detector simulations and/or data-
driven methods are used in order to obtain these corrections. Care
should be taken here in order to have a clear separation between
the parton-to-hadron level and hadron-to-detector level correc-
tions. Finally, for the sake of an easy comparison to the results of
other experiments and/or theoretical calculations, it is suggested
to provide, whenever possible, measurements corrected for detec-
tor effects and/or all necessary information related to the detector
response (e.g., the detector response matrix). Any fiducial phase
space for measurements should be defined as close as possible to
the detector-level selection in order to minimize model-dependent
extrapolations. A versatile repository for storing such information
is the use of Rivet routines [421].

9.3.1 Hadronic final-state observables

9.3.1.1 Jets

In hard interactions, final-state partons and hadrons appear
predominantly in collimated bunches, which are generically called
jets. To a first approximation, a jet can be thought of as a hard
parton that has undergone soft and collinear showering and then
hadronization. Jets are used both for testing our understanding
and predictions of high-energy QCD processes, and also for iden-
tifying the hard partonic structure of decays of massive particles
such as top quarks and W, Z and Higgs bosons.

In order to map observed hadrons onto a set of jets, one uses a
*jet definition*. The mapping involves explicit choices: for example
when a gluon is radiated from a quark, for what range of kinemat-
ics should the gluon be part of the quark jet, or instead form a
separate jet? Good jet definitions are infrared and collinear safe,
simple to use in theoretical and experimental contexts, applicable
to any type of inputs (parton or hadron momenta, charged par-
ticle tracks, and/or energy deposits in the detectors) and lead to
jets that are not too sensitive to non-perturbative effects.

An extensive treatment of the topic of jet definitions is given
in Ref. [422] (for $e^+e^-$ collisions) and Refs. [423–425]. Here we
briefly review the two main classes: cone algorithms, extensively
used at hadron colliders before the LHC, and sequential recombi-
nation algorithms, more widespread in $e^+e^-$ and ep colliders and
at the LHC.

Very generically, most (iterative) cone algorithms start with

some seed particle *i*, sum the momenta of all particles *j* within a
cone of opening-angle *R*, typically defined in terms of differences
in rapidity and azimuthal angle. They then take the direction of
this sum as a new seed and repeat until the direction of the cone
is stable, and call the contents of the resulting stable cone a jet
if its transverse momentum is above some threshold $p_{t,\text{min}}$. The
parameters *R* and $p_{t,\text{min}}$ should be chosen according to the needs
of a given analysis.

There are many variants of the cone algorithm, and they differ
in the set of seeds they use and the manner in which they ensure
a one-to-one mapping of particles to jets, given that two stable
cones may share particles (“overlap”). The use of seed particles is
a problem w.r.t. infrared and collinear safety. Seeded algorithms
are generally not compatible with higher-order (or sometimes even
leading-order) QCD calculations, especially in multi-jet contexts,
as well as potentially subject to large non-perturbative correc-
tions and instabilities. Seeded algorithms (JetCLU, MidPoint,
and various other experiment-specific iterative cone algorithms)
are therefore to be deprecated. Such algorithms are not used at
the LHC, but were at the Fermilab Tevatron, where data still
provide useful information, for example for global PDF fits.<fnref n="11" /> A
modern alternative is to use a seedless variant, SIScone [427].

Sequential recombination algorithms at hadron colliders (and in
DIS) are characterized by a distance $d_{ij} = \min(k_{t,i}^{2p}, k_{t,j}^{2p}) \Delta_{ij}^2 / R^2$
between all pairs of particles *i*, *j*, where $\Delta_{ij}$ is their separation
in the rapidity-azimuthal plane, $k_{t,i}$ is the transverse momen-
tum w.r.t. the incoming beams, and *R* is a free parameter. At
the LHC, *R* is typically in the range from 0.4 to 0.8, although
analyses can also use jet sizes up to 1.0-1.2. They also involve
a “beam” distance $d_{iB} = k_{t,i}^{2p}$. One identifies the smallest of
all the $d_{ij}$ and $d_{iB}$, and if it is a $d_{ij}$, then *i* and *j* are merged
into a new pseudo-particle (with some prescription, a recombina-
tion scheme, for the definition of the merged four-momentum). If
the smallest distance is a $d_{iB}$, then *i* is removed from the list of
particles and called a jet. As with cone algorithms, one usually
considers only jets above some transverse-momentum threshold
$p_{t,\text{min}}$. The parameter *p* determines the kind of algorithm: *p* = 1
corresponds to the (inclusive-) $k_t$ algorithm [296, 428, 429], *p* = 0
defines the *Cambridge-Aachen* algorithm [430, 431], while for the
*anti*- $k_t$ algorithm *p* = -1 [432]. All these variants are infrared and
collinear safe. Whereas the former two lead to irregularly shaped
jet boundaries, the latter results in cone-like boundaries, except
in situations where there are nearby jets. The *anti*- $k_t$ algorithm
has become the de-facto standard for the LHC experiments.

In $e^+e^-$ annihilation the $k_t$ algorithm [296] uses $y_{ij} = 2 \min(E_i^2, E_j^2)(1 - \cos\theta_{ij})/Q^2$ as distance measure between two particles/partons *i* and *j* and repeatedly merges the pair with smallest $y_{ij}$, until all $y_{ij}$ distances are above some threshold $y_{cut}$, the jet resolution parameter. *Q* is a measure of the overall hardness of the event. The (pseudo)-particles that remain at this point are called the jets. Here it is $y_{cut}$ (rather than *R* and $p_{t,\text{min}}$) that should be chosen according to the needs of the analysis. The two-jet rate in the $k_t$ algorithm has the property that logarithms $\ln(1/y_{cut})$ exponentiate. This is one reason why it is preferred over the earlier JADE algorithm [298], which uses the distance measure $y_{ij} = 2E_i E_j (1 - \cos\theta_{ij})/Q^2$. Note that other variants of sequential recombination algorithms for $e^+e^-$ annihilations, using different definitions of the resolution measure $y_{ij}$, exhibit much larger sensitivities to fragmentation and hadronization effects than the $k_t$ and JADE algorithms [433]. Efficient implementations of the above algorithms are available through the *FastJet* package [434].

While building infrared (IR) safe jets is generally considered a
solved problem, this is not quite the case when one also wishes
to assign a flavor to a jet. In fact, the general experimental def-
inition considers a flavored-jet to be a jet that contains at least
one flavor tag (such as a *B* or *D* meson) above a given trans-
verse momentum threshold. Because of collinear or soft wide-
angle *g* → $q\bar{q}$ splittings, it is easy to see that such a definition is
neither collinear nor infrared safe. This problem was addressed

<footnote n="11">¹¹In the data, the difference between the use of the Midpoint algorithm and the use of the SIScone algorithm is small [426], allowing the use of the SIScone algorithm for any theory comparisons at higher order in QCD.</footnote>
in Ref. [435,436] in the context of heavy-flavor production at the Tevatron. However, the jet-algorithm proposed in that work was impractical to implement experimentally because it was based on the $k_t$-algorithm. Furthermore, it required tagging two nearby flavored hadrons and involved a rather complex beam-distance measure. Since most experimental studies instead rely on the anti-$k_t$ algorithm, recent investigations have focused on developing algorithms that preserve the anti-$k_t$ kinematics of the jets while assigning jet flavors in an infrared safe way [437–439]. It turns out that the problem is more involved than anticipated, and a formulation of infrared-safe anti-$k_t$-like jets could only be achieved by introducing an interleaved flavor neutralization procedure [440]. An unfolding procedure will be necessary to convert experimental measurements of flavor-$k_t$ jets to a form that can be directly compared to theoretical predictions.

### 9.3.1.2 Event Shapes

Event-shape variables are functions of the four momenta of the particles in the final state and characterize the topology of an event’s energy flow. They are sensitive to QCD radiation (and correspondingly to the strong coupling) insofar as gluon emission changes the shape of the energy flow.

The classic example of an event shape is the thrust [441,442] in $e^+e^-$ annihilations, defined as

$$ \hat{\tau} = \max \frac{\sum_i |\vec{p}_i \cdot \vec{n}_\tau|}{\sum_i |\vec{p}_i|}, \qquad (9.22) $$

where $\vec{p}_i$ are the momenta of the particles or the jets in the final state and the maximum is obtained for the thrust axis $\vec{n}_\tau$. In the Born limit of the production of a perfect back-to-back $q\bar{q}$ pair, the limit $\hat{\tau} \to 1$ is obtained, whereas a perfectly spherical many-particle configuration leads to $\hat{\tau} \to 1/2$. Further event shapes of similar nature have been extensively measured at LEP and at HERA, and for their definitions and reviews we refer to Refs. [1, 8, 420, 443, 444]. The energy-energy correlation function [445], namely the energy-weighted angular distribution of produced hadron pairs, and its associated asymmetry are further shape variables which have been studied in detail at $e^+e^-$ colliders. For hadron colliders the appropriate modification consists in only taking the transverse momentum component [446]. The event shape variable N-jettiness has been proposed [447], that measures the degree to which the hadrons in the final state are aligned along N jet axes or the beam direction. It vanishes in the limit of exactly N infinitely narrow jets.

Phenomenological discussions of event shapes at hadron colliders can be found in Refs. [447–451]. Measurements of hadronic event-shape distributions have been published by CDF [452], ATLAS [453–459] and CMS [460–463].

Event shapes are used for many purposes. These include measuring the strong coupling (see e.g. Ref. [459]), tuning the parameters of Monte Carlo programs, investigating analytical models of hadronization and distinguishing QCD events from events that might involve decays of new particles (giving event-shape values closer to the spherical limit).

### 9.3.1.3 Jet substructure, quark vs. gluon jets

Jet substructure, which can be resolved by finding subjets or by measuring jet shapes, is sensitive to the details of QCD radiation in the shower development inside a jet and has been extensively used to study differences in the properties of quark and gluon induced jets, strongly related to their different color charges. There is clear experimental evidence that gluon jets have a softer particle spectrum and are “broader” than (light-) quark jets (as expected from perturbative QCD) when looking at observables such as the jet shape $\Psi(r/R)$ (see e.g. Ref. [464]). This is the fraction of transverse momentum contained within a sub-cone of cone-size $r$ for jets of cone-size $R$. It is sensitive to the relative fractions of quark and gluon jets in an inclusive jet sample and receives contributions from soft-gluon initial-state radiation and the underlying event. Therefore, it has been widely employed for validation and tuning of Monte Carlo parton-shower models. Furthermore, this quantity turns out to be sensitive to the modification of the gluon radiation pattern in heavy ion collisions (see e.g. Ref. [465]).

Jet shape measurements using proton-proton collision data have been presented for inclusive jet samples [466–468] and for top-quark production [469]. Further discussions, references and summaries can be found in Refs. [444,470–472] and Sec. 4 of Ref. [473].

The use of jet substructure has also been investigated in order to distinguish QCD jets from jets that originate from hadronic decays of boosted massive particles (high-$p_t$ electroweak bosons, top quarks and hypothesized new particles). A considerable number of experimental studies have been carried out with Tevatron and LHC data, in order to investigate the performance of the proposed algorithms for resolving jet substructure and to apply them to searches for new physics, as well as to the reconstruction of boosted top quarks, vector bosons and the Higgs boson. For reviews of this rapidly growing field, see Sec. 5.3 of Ref. [423], Refs. [350, 351, 473–479] and references thereto. One convenient representation for visualizing jet substructure is through the Lund plane [480]. From a theoretical perspective, the Lund plane would be constructed using quarks and gluons, but similar observables can be constructed through the use of jets, as described in detail in Ref. [481].

Neural network techniques and deep learning methods have also been applied to jet and top physics and jet substructure, see e.g. Refs. [482–485]. Perhaps no other sub-field has benefited as much from machine learning techniques as the study of jet substructure. As a jet can have O(100) constituents each with kinematic and other information, jet substructure analysis is naturally a highly multivariate problem. Deep learning techniques can use all of the available information to study jets in their natural high dimensionality. Such techniques have not only improved discrimination between different final states/types of jets, but have also improved our understanding of perturbative QCD. See for example the review in Ref. [349].

## 9.3.2 QCD measurements at colliders

There exists a wealth of data on QCD-related measurements in $e^+e^-$, ep, pp, and $p\bar{p}$ collisions, to which a short overview like this would not be able to do any justice. Reviews of the subject have been published in Refs. [443,444] for $e^+e^-$ colliders and in Ref. [486] for ep scattering, whereas for hadron colliders overviews are given in, e.g., Refs. [424, 471] and Refs. [2, 11, 487–489].

Below we concentrate our discussion on measurements that are most sensitive to hard QCD processes with focus on jet production.

### 9.3.2.1 $e^+e^-$ colliders

Analyses of jet production in $e^+e^-$ collisions are mostly based on data from the JADE experiment at center-of-mass energies between 14 and 44 GeV, as well as on LEP collider data at the Z resonance and up to 209 GeV. The analyses cover the measurements of (differential or exclusive) jet rates (with multiplicities typically up to 4, 5 or 6 jets), the study of three-jet events and particle production between the jets, as well as four-jet production and angular correlations in four-jet events.

Event-shape distributions from $e^+e^-$ data have been an important input to the tuning of parton shower MC models, typically matched to matrix elements for three-jet production. In general these models provide good descriptions of the available, highly precise data. Especially for the large LEP data sample at the Z peak, the statistical uncertainties are mostly negligible and the experimental systematic uncertainties are at the percent level or even below. These are usually dominated by the uncertainties related to the MC model dependence of the efficiency and acceptance corrections (often referred to as “detector corrections”).

Observables measured in $e^+e^-$ collisions have been used for determinations of the strong coupling constant (cf. Section 9.4 below) and for putting constraints on the QCD color factors (cf. Sec. 9.1 for their definitions), thus probing the non-Abelian nature of QCD. Angular correlations in four-jet events are sensitive at leading order. Some sensitivity to these color factors, although only at NLO, is also obtained from event-shape distributions. Scaling violations of fragmentation functions and the different subjet structure in quark and gluon induced jets also give access to these color factors. A compilation of results [444] quotes world average values of $C_A = 2.89 \pm 0.03 (\text{stat}) \pm 0.21 (\text{syst})$
and $C_F = 1.30 \pm 0.01$ (stat) $\pm 0.09$ (syst), with a correlation co-efficient of 82%. These results are in perfect agreement with the expectations from SU(3) of $C_A = 3$ and $C_F = 4/3$.

### 9.3.2.2 DIS and photoproduction

Jet measurements in ep collisions, both in the DIS and photoproduction regimes, allow for tests of QCD factorization (as they involve only one initial state proton and thus one PDF function), and provide sensitivity to both the gluon PDF and to the strong coupling constant. Calculations are available at NNLO in both regimes [219,220]. An $N^3LO$ calculation using the Projection-to-Born method was also presented in Ref. [221,222]. Experimental uncertainties of the order of 5–10% have been achieved, whereas statistical uncertainties are negligible to a large extent. For comparison to theoretical predictions, at large jet $p_t$ the PDF uncertainty dominates the theoretical uncertainty (typically of order 5–10%, in some regions of phase space up to 20%), therefore jet observables become useful inputs for PDF fits.

In general, the data are well described by the NLO and NNLO matrix-element calculations, combined with DGLAP evolution equations, in particular at large $Q^2$ and central values of jet pseudo-rapidity. At low values of $Q^2$ and $x$, in particular for large jet pseudo-rapidities, certain features of the data have been interpreted as requiring BFKL-type evolution, though the predictions for such schemes are still limited. It is worth noting that there is no indication that the BKFL approximation is needed within the currently probed phase space in the $x, Q^2$ plane, and an alternative approach [490], which implements the merging of LO matrix-element based event generation with a parton shower (using the SHERPA framework), successfully describes the data in all kinematical regions, including the low $Q^2$, low $x$ domain. At moderately small $x$ values, it should perhaps not be surprising that the BFKL approach and fixed-order matrix-element merging with parton showers may both provide adequate descriptions of the data, because some part of the multi-parton phase space that they model is common to both approaches.

In the case of photoproduction, a wealth of measurements with low $p_t$ jets were performed in order to constrain the photon content of the proton (which is by now determined with percent accuracy thanks to the LUX approach [72,73]). A few examples of measurements can be found in Refs. [491–495] for photoproduction and in Refs. [496–505] for DIS.

### 9.3.2.3 Hadron-hadron colliders

The spectrum of observables and the number of measurements performed at hadron colliders is enormous, probing many regions of phase space and covering a huge range of cross sections, as exemplified in Fig. 9.1 for a wide class of processes by the ATLAS experiment, and specialised to top-quark related processes by the CMS experiment at the LHC. In general, the theory agreement with data is excellent for a wide variety of processes, indicating the success of perturbative QCD with the PDF and the strong coupling as inputs. For the sake of brevity, in the following only certain classes of those measurements will be discussed, which permit various aspects of the QCD studies to be addressed. Most of our discussion will focus on LHC results, which are available for center-of-mass energies of 2.76, 5, 7, 8 and 13 TeV with integrated luminosities of up to $140 \text{ fb}^{-1}$. As of writing of this update, new results at 13.6 TeV are starting to appear. Generally speaking, besides representing precision tests of the Standard Model and QCD in particular, these measurements serve several purposes, such as: (i) probing pQCD and its various approximations and implementations in MC models to quantify the order of magnitude of not yet calculated contributions and to gauge their precision when used as background predictions to searches for new physics, or (ii) extracting/constraining model parameters such as the strong coupling constant or PDFs.

The final states measured at the LHC include single, double and triple gauge boson production, top production (single top, top pair and four top production), Higgs boson production, alone and in conjunction with a W or Z boson, and with a top quark pair. Many/most of these events are accompanied by additional jets. So far only limits have been placed on double Higgs production. The volume of LHC results prohibits a comprehensive description

in this Review; hence, only a few highlights will be presented.

Among the most important cross sections measured, and the one with the largest dynamic range, is the inclusive jet spectrum as a function of the jet transverse momentum ($p_t$), for several rapidity regions and for $p_t$ up to 700 GeV at the Tevatron and $\sim 3.5$ TeV at the LHC. It is worth noting that this latter upper limit in $p_t$ corresponds to a distance scale of $\sim 10^{-19}$ m: no other experiment so far is able to directly probe smaller distance scales of nature than this measurement. The Tevatron inclusive jet measurements in Run 2 (Refs. [508–511]) were carried out with the MidPoint jet clustering algorithm (or its equivalent) and, in a few cases, with the $k_t$ jet clustering algorithm. Most of the LHC measurements use the anti-$k_t$ algorithm, with a variety of jet radii. The use of multiple jet radii in the same analysis allows a better understanding of the underlying QCD dynamics. Measurements by ALICE, ATLAS and CMS have been published in Refs. [512–522].

In general, we observe a good description of the data by the NLO and NNLO QCD predictions over about 11 orders of magnitude in cross section, as long as care is taken for the functional form of the central scale choice [523]. The experimental systematic uncertainties are dominated by the jet energy scale uncertainty, quoted to be in the range of a few percent (see for instance the review in Ref. [524]), leading to uncertainties of $\sim 5–30\%$ on the cross section, increasing with $p_t$ and rapidity. The PDF uncertainties dominate the theoretical uncertainty at large $p_t$ and rapidity. In fact, inclusive jet data are one of the most important inputs to global PDF fits, in particular for constraining the high-$x$ gluon PDF [525–527]. Constraints on the PDFs can also be obtained from ratios of inclusive cross sections at different center-of-mass energies [513,518]. Ratios of jet cross sections are a means to (at least partially) cancel the jet energy scale uncertainties and thus provide jet observables with significantly improved precision.

Dijet events are typically analyzed in terms of their invariant mass or average dijet $p_t$ and angular distributions, which allows for tests of NLO and NNLO QCD predictions (see e.g. Refs. [517, 528, 529] for recent LHC results), and for setting stringent limits on deviations from the Standard Model, such as quark compositeness or contact interactions (some examples can be found in Refs. [520, 530–536]). Furthermore, dijet azimuthal correlations between the two leading jets, normalized to the total dijet cross section, are an extremely valuable tool for studying the spectrum of gluon radiation in the event. The azimuthal separation of the two leading jets is sensitive to multi-jet production, avoiding at the same time large systematic uncertainties from the jet energy calibration. For example, results from the Tevatron [537,538] and the LHC [458,539–543] show that the LO (non-trivial) prediction for this observable, with at most three partons in the final state, is not able to describe the data for an azimuthal separation below $2\pi/3$, where NLO contributions (with 4 partons) restore the agreement with data. In addition, this observable can be employed to tune Monte Carlo predictions of soft gluon radiation. Further examples of dijet observables that probe special corners of phase space are those that involve forward (large rapidity) jets and where a large rapidity separation, possibly also a rapidity gap, is required between the two jets. Reviews of such measurements can be found in Ref. [471], showing that no single prediction is capable of describing the data in all phase-space regions. In particular, no conclusive evidence for BFKL effects in these observables has been established so far.

Beyond dijet final states, measurements of the production of three or more jets, including cross section ratios, have been performed (see Refs. [471, 544] for recent reviews), as a means of testing perturbative QCD predictions, determining the strong coupling constant, and probing/tuning MC models, in particular those combining multi-parton matrix elements with parton showers [459, 545]. The calculation of three-jet production to NNLO [264] will allow more precise predictions of the three-jet to two-jet ratio, and thus the extraction of $\alpha_s(m_Z^2)$ from this observable. This calculation also allows for the use of transverse energy-energy correlations (TEEC) for the determination of $\alpha_s(m_Z^2)$ [459].

W and Z production serve as benchmark cross sections at the
Standard Model Production Cross Section Measurements

Status: February 2022

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

All results at: http://cern.ch/go/pNj7

**Figure 9.1:** Overview of cross section measurements for a wide class of processes by the ATLAS [506] experiment, and specialised to top-quark related processes by the CMS [507] experiment at the LHC, for center-of-mass energies of 5, 7, 8, 13, and 13.6 TeV. Also shown are the theoretical predictions and their uncertainties.
LHC. The large boson mass stabilizes perturbative predictions, which results in better theoretical accuracy. In terms of experimental precision, measurements of inclusive vector boson (W, Z) production provide the most precisely determined observables at hadron colliders so far. This is because the experimental signatures are based on charged leptons which are measured much more accurately than jets or photons. At the LHC [546–553], the dominant uncertainty stems from the luminosity determination (≤2–4%), while other uncertainties (e.g. statistical or from lepton efficiencies) are controlled at the ~0.5–3% level. The uncertainty from the acceptance correction of about ~1–2% can be reduced by measuring so-called fiducial cross sections, i.e. by applying kinematic cuts also to the particle level of the theoretical predictions. Measurements combining the electron and muon final states were able to achieve a precision a few per mille level in the dilepton final state for the normalized cross section for $p_T^l$ less than 30 GeV [554, 555]. This level of precision can also be achieved by measuring cross section ratios (W/Z or $W^+/W^-$). On the theory side, as discussed earlier in this Review, the production of these color-singlet states has been calculated to N³LO [119–121]. Since the dominant theoretical uncertainty is related to the choice of PDFs, these high-precision data provide useful handles for PDF determinations.

Further insights are obtained from measurements of differential vector boson production, as a function of the invariant dilepton mass, the boson's rapidity or its transverse momentum. For example, the dilepton invariant mass distribution has been measured [556–561] for masses between 15 and 3000 GeV, covering more than 8 orders of magnitude in cross section. NNLO QCD predictions, together with modern PDF sets and including higher-order electroweak and QED final-state radiation corrections, describe the data to within 5–10% over this large range, whereas NLO predictions show larger deviations, unless matched to a parton shower.

Similar conclusions can be drawn from the observed rapidity distribution of the dilepton system (see e.g. Refs. [546, 557, 562]) or, in the case of W production, from the observed charged lepton rapidity distribution and its charge asymmetry. The latter is particularly sensitive to differences among PDF sets [546, 563–565], also thanks to the high precision achieved by the ATLAS and CMS experiments for central rapidity ranges. These measurements are extended to the very forward region, up to 4.5 in lepton rapidity, by the LHCb experiment [552, 553, 566].

An overview of these kinds of measurement can be found in Ref. [471]. There one can also find a discussion of and references to LHC results from studies of the vector boson's transverse momentum distribution, $p_t^V$ (see also Refs. [567–569]). This observable covers a wide kinematic range and probes different aspects of higher-order QCD effects. It is sensitive to jet production in association with the vector boson, without suffering from the large jet energy scale uncertainties. In the $p_t^V$ region of several tens of GeV to over 1 TeV, the NNLO predictions for V+jet can be used to predict the high-$p_t$ boson production cross section.<fnref n="12" /> The NNLO predictions agree with the data to within about 10%, and agree somewhat better at high transverse momentum than do the NLO predictions [570]. At transverse momenta below ~20 GeV, the fixed-order predictions fail and soft-gluon resummation is needed to restore the agreement with data. The soft gluon resummation can either be performed analytically, or effectively using parton showering implemented in Monte Carlo programs. While analytic approaches reach a higher perturbative precision, they typically refer to inclusive measurements without fiducial cuts.

The addition of further jets to the final state extends the kinematic range as well as increasing the complexity of the calculation/measurements. The number of results obtained both at the Tevatron and at the LHC is extensive, summaries can be found in Refs. [471, 571], and more recent results can be found in Refs. [570, 572–575]. The measurements cover a very large phase space, e.g. with jet transverse momenta between 30 GeV and ~1.5 TeV and jet rapidities up to $|y| < 4.4$ [570, 576]. Jet

multiplicities as high as seven jets accompanying the vector boson have already been probed at the LHC, together with a substantial number of other kinematical observables, such as angular correlations among the various jets or among the jets and the vector boson, or the sum of jet transverse momenta, $H_T$. Whereas the jet $p_t$ and $H_T$ distributions are dominated by jet energy scale uncertainties at levels similar to those discussed above for inclusive jet production, angular correlations and jet multiplicity ratios have been measured with a precision of ~10%, see e.g. Refs. [461, 577].

NLO calculations for up to five jets [578] in addition to the vector boson are in good agreement with the data over that phase space, where the calculations are applicable. Predictions for V+jet at NNLO improve the description of the data for distributions involving the vector boson or the leading jet. MC models that implement parton shower matching to matrix elements (either at LO or NLO) have mixed results.

The challenges get even more severe in the case of vector boson plus heavy quark (b, c) production: on the theory side because an additional scale is introduced by the heavy quark mass, and different schemes exist for the handling of heavy quarks and their mass effects in the initial and/or final state; and on the experimental side because additional uncertainties related to the heavy-flavor tagging must be considered [579] (see also the discussion regarding flavor-jet definition in Sec. 9.3.1.1). A review of heavy quark production at the LHC is presented in Ref. [580], where the di-$b$-jet $p_t$ and mass spectra are found to be well modeled within uncertainties by most generators for b-jet production with or without associated W and Z bosons. However, sizable differences between data and predictions are seen in the modeling of events with single b jets, particularly at large b-jet $p_t$, where gluon splitting processes become dominant, as also confirmed by studies of b-hadron and b-jet angular correlations.

The precision reached in photon measurements is in between that for lepton and jet measurements. The photon energy and angles can be measured at about the same precision as the lepton energy and angles in Drell-Yan production, but there are greater challenges encountered in photon reconstruction (for example isolation) and in purity determination. Note, though, that the lepton purity approaches unity as the photon $p_t$ increases. At high $p_t$, it becomes increasingly difficult for a jet to fragment into an isolated neutral electromagnetic cluster which mimics the lepton signature. The inclusive photon cross section can be measured [528, 581–584], as well as the production of a photon accompanied by one or more jets [584–586, 586–590]. The kinematic range for photon production is less than that for jet production because of the presence of the electromagnetic coupling, but still reaches about 2 TeV. Better agreement is obtained with NNLO predictions for photon production than for NLO predictions, except when the latter are matched to matrix element plus parton shower predictions. Photon production in association with a heavy-flavor jet is a useful input for the determination of the b and c quark PDFs [591].

Electroweak corrections are expected to become more and more relevant now that the TeV energy range starts to be explored, and EW corrections can now also be computed automatically [171]. For a comprehensive review on electroweak corrections see Ref. [192]. For example, such corrections were found [592, 593] to be sizable (tens of percent) when studying the ratio $(d\sigma^\gamma/dp_t)/(d\sigma^Z/dp_t)$ in $\gamma(Z)+jet$ production, $p_t$ being the boson's transverse momentum, and might account for (some of) the differences observed in a CMS measurement [594] of this quantity.

A number of interesting developments, in terms of probing higher-order QCD effects, have occurred in the sector of diboson production, in particular for the WW and $\gamma\gamma$ cases. Regarding the former, an early disagreement of about 10% between the LHC measurements and the NLO predictions had led to a number of speculations of possible new physics effects in this channel. However, more recent ATLAS and CMS measurements [595–598] are in agreement with the NNLO prediction [131]. The statistical reach of the LHC has resulted in the discovery of triple massive gauge boson production [599–601].

In the case of diphoton production, ATLAS [602–604] and CMS [605] have provided accurate measurements, in particular

<footnote n="12">¹²For these calculations, there is a requirement of the presence of a jet, but the $p_t$ cut is typically small (30 GeV) compared to the high $p_t$ region being discussed here.</footnote>
for phase-space regions that are sensitive to radiative QCD corrections (multi-jet production), such as small azimuthal photon separation. While there are large deviations between data and NLO predictions in this region, a calculation [228] at NNLO accuracy manages to mostly fill this gap. This is an interesting example where scale variations can not provide a reliable estimate of missing contributions beyond NLO, since at NNLO new channels appear in the initial state (gluon fusion in this case). These missing channels can be included in a matrix element plus parton shower calculation in which two additional jets are included at NLO. The result exhibits a similar level of agreement as that obtained at NNLO. Three-photon production has also been measured [606] and is in good agreement with NNLO theory predictions [260, 261].

In terms of heaviest particle involved, top-quark production at the LHC has become an important tool for probing higher-order QCD calculations, thanks to very impressive achievements both on the experimental and theoretical side, as extensively summarized in Ref. [607]. Regarding $t\bar{t}$ production, the most precise inclusive cross section measurements are achieved using the dilepton ($e\mu$) final state, with a total uncertainty of 4% [608–612]. This is of about the same size as the uncertainty on the most advanced theoretical predictions [128, 129, 613, 614], obtained at NNLO with additional soft-gluon resummation at NNLL accuracy [615]. There is excellent agreement between data and the QCD predictions.

The $t\bar{t}$ final state allows multiple observables to be measured. A large number of differential cross section measurements have been performed at 7, 8 and 13 TeV center-of-mass energy, studying distributions such as the top-quark $p_t$ and rapidity, the transverse momentum and invariant mass of the $t\bar{t}$ system (probing the TeV range), or the number of additional jets. These measurements have been compared to a wide range of predictions, at fixed order up to NNLO as well as using LO or NLO matrix elements matched to parton showers. Each of the observables provides information on the high $x$ gluon. However, there are tensions among the multiple observables that can lead to difficulties in PDF fits. Four top production has been measured in multi-lepton final states [616] by the ATLAS and CMS collaboration, with results consistent with the standard model prediction [617].

Thanks to both the precise measurements of, and predictions for, the inclusive top-pair cross section, which is sensitive to the strong coupling constant and the top-quark mass, this observable has been used to measure the strong coupling constant at NNLO accuracy from hadron collider data [618, 619] (cf. Section 9.4 below), as well as to obtain a measurement of the top-quark's pole mass without employing direct reconstruction methods [618, 620, 621].

The Higgs boson provides a tool for QCD studies, especially as the dominant production mechanism is gg fusion, which is subject to very large QCD corrections. Higgs boson production has been measured in the ZZ, $\gamma\gamma$, $b\bar{b}$, WW and $\tau\tau$ decay channels. A measurement of the cross section in the ZZ and $\gamma\gamma$ channels is one of the first measurements at 13.6 TeV to be published [622]. The experimental cross section is now known with a precision approaching 10% [623, 624], similar to the size of the theoretical uncertainty [625], of which the PDF+$\alpha_s$ uncertainty is the largest component. Part of this systematic is the mis-match in orders between the PDF determination (NNLO) and the cross section evaluation (N³LO), as discussed earlier. The experimental precision has allowed detailed fiducial and differential cross section measurements. For example, with the diphoton final state, the transverse momentum of the Higgs boson can be measured out to the order of 650 GeV [626–629], where top quark mass effects become important. The production of a Higgs boson with up to 4 jets has been measured [624, 626]. The experimental cross sections have been compared to NNLO predictions (for $H \ge 1$ jet), NLO for 2 and 3 jets, and NNLO+NNLL for the transverse momentum distribution. In addition, finite top quark mass effects have been taken into account at NLO. The use of the boosted $H \to b\bar{b}$ topology allows probes of Higgs boson transverse momenta on the order of 600 GeV and greater [624]. So far the agreement with the perturbative QCD corrections is good.

## 9.4 Determinations of the strong coupling constant

Beside the quark masses, the only free parameter in the QCD Lagrangian is the strong coupling constant $\alpha_s$. The coupling constant in itself is not a physical observable, but rather a quantity defined in the context of perturbation theory, which enters predictions for experimentally measurable observables, such as $R$ in Eq. (9.7). The value of the strong coupling constant must be inferred from such measurements and is subject to experimental and theoretical uncertainties. The incomplete knowledge of $\alpha_s$ propagates into uncertainties in numerous precision tests of the Standard Model. Here, we present an update of the 2022 PDG average value of $\alpha_s(m_Z^2)$ and its uncertainty [630].

Many experimental observables are used to determine $\alpha_s$. A number of recent determinations are collected in Refs. [631, 632]. Further discussions and considerations on determinations of $\alpha_s$ can also be found in Refs. [633–635]. Such considerations include:

* The observable's sensitivity to $\alpha_s$ as compared to the experimental precision. For example, for the $e^+e^-$ cross section to hadrons (cf. R in Sec. 9.2.1), QCD effects are only a small correction, since the perturbative series starts at order $\alpha_s^0$; three-jet production or event shapes in $e^+e^-$ annihilations are directly sensitive to $\alpha_s$ since they start at order $\alpha_s$; the hadronic decay width of heavy quarkonia, $\Gamma(\gamma \to \text{hadrons})$, is very sensitive to $\alpha_s$ since its leading order term is $\propto \alpha_s^3$.

* The accuracy of the perturbative prediction, or equivalently of the relation between $\alpha_s$ and the value of the observable. Several observables have been known to NNLO for quite some time. These include, for instance, inclusive observables, as well as three-jet rates and event shapes in $e^+e^-$ collisions, inclusive jet and dijet production in DIS, and inclusive jet, dijet, $t\bar{t}$, $W/Z$+jet and three-jet production cross sections in $pp$ or $p\bar{p}$ collisions. The $e^+e^-$ hadronic cross section and $\tau$, $W$ and $Z$ branching fractions to hadrons are even known to N³LO, if one denotes the LO as the first non-trivial term. In certain cases, fixed-order predictions are supplemented with resummation. The precise magnitude of the associated theory uncertainties usually is estimated as discussed in Sec. 9.2.4.

* The size of non-perturbative effects. Sufficiently inclusive quantities, like the $e^+e^-$ cross section to hadrons, have small non-perturbative contributions $\sim \Lambda^4/Q^4$. Others, such as event-shape distributions, have typically contributions $\sim \Lambda/Q$.

* The scale at which the measurement is performed. An uncertainty $\delta$ on a measurement of $\alpha_s(Q^2)$, at a scale $Q$, translates to an uncertainty $\delta' = (\alpha_s^2(m_Z^2)/\alpha_s^2(Q^2)) \cdot \delta$ on $\alpha_s(m_Z^2)$. For example, this enhances the already important impact of precise low-$Q$ measurements, such as from $\tau$ decays, in combinations performed at the $m_Z$ scale.

The selection of results from which to determine the world average value of $\alpha_s(m_Z^2)$ is restricted to those that are

- published in a peer-reviewed journal at the time of writing this report,

- based on the most complete perturbative QCD predictions of at least NNLO accuracy,

- accompanied by reliable estimates of all experimental and theoretical uncertainties.

Numerous measurements from jet production in DIS and at hadron colliders are still excluded from the average presented here, because the determination of $\alpha_s(m_Z^2)$ from those data sets has not yet been upgraded to NNLO. A few new results did appear comparing to theory at NNLO and are included in the corresponding section. NLO analyses will still be discussed in this Review, as they are important ingredients for the experimental evidence of the energy dependence of $\alpha_s$, i.e. for asymptotic freedom, one of the key features of QCD.

In order to calculate the world average value of $\alpha_s(m_Z^2)$, we apply, as in earlier editions, an intermediate step of pre-averaging results within the sub-fields now labeled “Hadronic $\tau$ decays and low $Q^2$ continuum” ($\tau$ decays and low $Q^2$), “Heavy quarkonia decays”
(QQ̃ bound states), “PDF fits” (PDF fits), “Hadronic final states of $e^{+}e^{-}$ annihilations” ($e^{+}e^{-}$ jets & shapes), “Observables from hadron-induced collisions” (hadron colliders), and “Electroweak precision fit” (electroweak) as explained in the following sections. For each sub-field, the *unweighted average* of all selected results is taken as the pre-average value of $\alpha_s(m_Z^2)$, and the *unweighted average* of the quoted total uncertainties is assigned to be the respective overall error of this pre-average. Asymmetric total uncertainties are symmetrised beforehand by adopting the larger of the two values as the (±) uncertainty. For the “Lattice QCD” (lattice) sub-field we do not perform a pre-averaging; instead, we adopt for this sub-field the $\alpha_s(m_Z^2)$ estimate and uncertainty derived by the Flavour Lattice Averaging Group (FLAG) in Ref. [636].

Assuming that the six sub-fields (excluding lattice) are largely independent of each other, we determine a non-lattice world average value using a 'χ² averaging' method. In a last step we perform an unweighted average of the values and uncertainties of $\alpha_s(m_Z^2)$ from our non-lattice result and the lattice result presented in the FLAG 2021 report [636].

### 9.4.1 Hadronic $\tau$ decays and low $Q^2$ continuum:

Based on complete N³LO predictions [34], analyses of the $\tau$ hadronic decay width and spectral functions have been performed, e.g. in Refs. [34,637-642], and lead to precise determinations of $\alpha_s$ at the energy scale of $m_\tau^2$. They are based on different approaches to treat perturbative and non-perturbative contributions, the impacts of which have been a matter of intense discussions for a long time, see e.g. Refs. [641-644]. In particular, in $\tau$ decays there is a significant difference between results obtained using fixed-order (FOPT) or contour-improved perturbation theory (CIPT), such that analyses based on CIPT generally arrive at larger values of $\alpha_s(m_\tau^2)$ than those based on FOPT.

In addition, some results show discrepancies in $\alpha_s(m_\tau^2)$ among groups of authors using the same data sets and perturbative calculations, most likely due to different treatments of the non-perturbative contributions, cf. Ref. [642] with Refs. [641, 645]. References [646, 647] question the validity of using a truncated operator product expansion (OPE) at $Q^2 = m_\tau^2$ based on the disagreement found between experimental values of the spectral moments and the theory representations based on the truncated OPE fits at $Q^2 > m_\tau^2$.

Recent developments now have shed light on the disagreement between FOPT- and CIPT-based analyses. References [648-650] argue that CIPT-based calculations require a dedicated estimation of non-perturbative effects instead of “standard” ones used so far. Otherwise an asymptotic separation between results using the two perturbative approaches would remain. Potential ways forward are suggested in Refs. [651, 652]. As a consequence we remove for the time being determinations of the strong coupling constant based on CIPT from the derivation of the central value and the uncertainty of $\alpha_s(m_Z^2)$.

We determine the pre-average value of $\alpha_s(m_Z^2)$ for this sub-field only from studies that employ FOPT expansions and remove any eventual averaging with CIPT central values or increased uncertainties due to the differences in CIPT vs. FOPT. As the results from Refs. [34, 644, 645] are not totally independent, we pre-averaged as a first step these three results in the previous edition of this Review. Lacking, however, an estimate of the theory uncertainty for the FOPT method, Ref. [645] had to be left out, such that the first entry to this category of $\alpha_s$ determinations is pre-averaged from $\alpha_s(m_Z^2) = 0.1183 \pm 0.0026$ [34] and $\alpha_s(m_Z^2) = 0.1181 \pm 0.0015$ [644] to $\alpha_s(m_Z^2) = 0.1182 \pm 0.0021$ (summarized as BP 2008-16 FO in Fig. 9.2).

Subsequently, this is combined with $\alpha_s(m_Z^2) = 0.1158 \pm 0.0022$ [653] and $\alpha_s(m_Z^2) = 0.1171 \pm 0.0010$ [654], which replaces the previous result from Ref. [642]. We also include the result from $\tau$ decay and lifetime measurements, obtained in Sec. Electroweak Model and constraints on New Physics of the 2022 edition of this Review [630], $\alpha_s(m_Z^2) = 0.1171 \pm 0.0018$. The latter result, being a global fit of $\tau$ data, involves some correlations with the other extractions of this category. However, since we perform an unweighted average of the central value and uncertainty, the effects of the potential correlations are reduced. Finally, a new

<figure><img src="image_18.png" /><figcaption>**Figure 9.2:** Summary of determinations of $\alpha_s(m_Z^2)$ with uncertainty in the seven sub-fields as discussed in the text. The yellow (light shaded) bands and dotted lines indicate the pre-average values of each sub-field. The dashed line and blue (dark shaded) band represent the final world average value of $\alpha_s(m_Z^2)$. The "*" symbol within the “hadron colliders” sub-field indicates a determination including a simultaneous fit of PDFs.</figcaption></figure>

determination reported in Ref. [655] is added to this category: $\alpha_s(m_Z^2) = 0.1183_{-0.0012}^{+0.0009}$. A recent publication evaluating data on $R_{had}$ in the continuum below the charm threshold exhibits huge experimental uncertainties and therefore has not been considered [656].

All these results are summarized in Fig. 9.2. Determining the unweighted average of the central values and their overall uncertainties, we arrive at $\alpha_s(m_Z^2) = 0.1173 \pm 0.0017$, which we will use as the first input for determining the world average value of $\alpha_s(m_Z^2)$. This corresponds to $\alpha_s(m_Z^2) = 0.314 \pm 0.014$.

### 9.4.2 Heavy quarkonia decays:

Two determinations of $\alpha_s(m_Z^2)$ have been performed [657, 658] that are based on N³LO accurate predictions. Reference [657] performs a simultaneous fit of the strong coupling and the bottom mass $m_b$, including states with principal quantum number up to $n \le 2$ in order to break the degeneracy between $\alpha_s$ and $m_b$, finding $\alpha_s(m_Z^2) = 0.1178 \pm 0.0051$. Reference [658] instead uses as input of the fit the renormalon-free combination of masses of the meson $B_c$, the bottomonium $\eta_b$ and the charmonium $\eta_c$, $M_{B_c} - M_{\eta_b}/2 - M_{\eta_c}/2$, which is weakly dependent on the heavy quark masses, but shows a good dependence on $\alpha_s$. Using this observable, they obtain $\alpha_s(m_Z^2) = 0.1195 \pm 0.0053$.
Two further values are derived at NNLO in Ref. [659,660] from mass splittings and sum rules giving $\alpha_s(m_Z^2) = 0.1183 \pm 0.0019$ and $\alpha_s(m_Z^2) = 0.1175 \pm 0.0032$ when evolved from the relevant charmonium respectively bottomonium mass scales to $m_Z^2$. Finally, by means of quarkonium sum rules, Refs. [661,662] quote $\alpha_s(m_Z^2) = 0.1168 \pm 0.0019$ and $\alpha_s(m_Z^2) = 0.1186 \pm 0.0048$ for charmonium and bottomonium respectively. These six determinations satisfy our criteria to be included in the heavy-quarkonia category of the world average. Their unweighted combination leads to the pre-average for this category of $\alpha_s(m_Z^2) = 0.1181 \pm 0.0037$.

### 9.4.3 PDF fits:

Another class of studies, analyzing structure functions at NNLO QCD (and partly beyond), provide results that serve as relevant inputs for the world average of $\alpha_s$. Some of these studies do not, however, explicitly include estimates of theoretical uncertainties when quoting fit results of $\alpha_s$. In such cases we add, in quadrature, half of the difference between the results obtained in NNLO and NLO to the quoted errors.

A combined analysis of non-singlet structure functions from DIS [663], based on QCD predictions up to N³LO in some of its parts, results in $\alpha_s(m_Z^2) = 0.1141 \pm 0.0022$ (BBG). Studies of singlet and non-singlet structure functions, based on NNLO predictions, result in $\alpha_s(m_Z^2) = 0.1162 \pm 0.0017$ [664] (JR14). The AMBP group [665,666] determined a set of parton distribution functions using data from HERA, NOMAD, CHORUS, from the Tevatron and the LHC using the Drell-Yan process and the hadro-production of single-top and top-quark pairs, and determined $\alpha_s(m_Z^2) = 0.1147 \pm 0.0024$ [665].

The MSHT group [667], also including hadron collider data, determined a set of parton density functions (MSHT20) together with $\alpha_s(m_Z^2) = 0.1174 \pm 0.0013$. Similarly, the CT group [37] determined the CT18 parton density set together with $\alpha_s(m_Z^2) = 0.1164 \pm 0.0026$. The NNPDF group [668] presented NNPDF3.1 parton distribution functions together with $\alpha_s(m_Z^2) = 0.1185 \pm 0.0012$. New in this category is the result of $\alpha_s(m_Z^2) = 0.1156 \pm 0.0031$ reported as HERAPDF2.0Jets in Ref. [669]. In addition to DIS data from HERA also jet cross sections are considered in the fit to theory at NNLO.

We note that criticism has been expressed on some of the above extractions. Among the issues raised, we mention the neglect of singlet contributions at $x \ge 0.3$ in pure non-singlet fits [670], the impact and detailed treatment of particular classes of data in the fits [670,671], possible biases due to insufficiently flexible parameterizations of the PDFs [672] and the use of a fixed-flavor number scheme [673,674]. For recent more extensive discussions see e.g. Ref. [38].

Summarizing the results from world data on structure functions, taking the unweighted average of the central values and errors of all selected results, leads to a pre-average value of $\alpha_s(m_Z^2) = 0.1161 \pm 0.0022$, see also Fig. 9.2.

### 9.4.4 Hadronic final states of $e^+e^-$ annihilation:

Re-analyses of jets and event shapes in $e^+e^-$ annihilation (j&s), measured around the Z peak and at LEP2 center-of-mass energies up to 209 GeV, using NNLO predictions matched to NLL resummation and Monte Carlo models to correct for hadronization effects, resulted in $\alpha_s(m_Z^2) = 0.1224 \pm 0.0039$ (ALEPH) [675], and in $\alpha_s(m_Z^2) = 0.1189 \pm 0.0043$ (OPAL) [676]. Similarly, an analysis of JÄDE data [677] at center-of-mass energies between 14 and 46 GeV gives $\alpha_s(m_Z^2) = 0.1172 \pm 0.0051$, with contributions from the hadronization model and from perturbative QCD uncertainties of 0.0035 and 0.0030, respectively. Precise determinations of $\alpha_s$ from three-jet production alone (3j), at NNLO, resulted in $\alpha_s(m_Z^2) = 0.1175 \pm 0.0025$ [678] from ALEPH data and in $\alpha_s(m_Z^2) = 0.1199 \pm 0.0059$ [679] from JADE. A recent determination is based on an NNLO+NNLL accurate calculation that allows to fit the region of lower three-jet rate (2j) using data collected at LEP and PETRA at different energies. This fit gives $\alpha_s(m_Z^2) = 0.1188 \pm 0.0013$ [680], where the dominant uncertainty is the hadronization uncertainty, which is estimated from Monte Carlo simulations. A fit of energy-energy-correlation (EEC), also based on an NNLO+NNLL calculation, together with

a Monte Carlo based modeling of hadronization corrections gives $\alpha_s(m_Z^2) = 0.1175 \pm 0.0029$ [681]. These results are summarized in the $e^+e^-$ sector of Fig. 9.2.

Another class of $\alpha_s$ determinations is based on analytic modeling of non-perturbative and hadronization effects, rather than on Monte Carlo models [682–685], using methods like power corrections, factorization of soft-collinear effective field theory, dispersive models and low scale QCD effective couplings. In these studies, the world data on thrust distributions (T), or the C-parameter distributions (C), are analyzed and fitted to perturbative QCD predictions at NNLO matched with resummation of leading logs up to N³LL accuracy, see Sec. 9.2.3.3. The results are $\alpha_s(m_Z^2) = 0.1135 \pm 0.0011$ [683] and $\alpha_s(m_Z^2) = 0.1134^{+0.0031}_{-0.0025}$ [684] from thrust, and $\alpha_s(m_Z^2) = 0.1123 \pm 0.0015$ [685] from C-parameter.

A long-standing question was why this latter class of determinations led systematically to rather small values of $\alpha_s(m_Z^2)$ as compared to the former one based on Monte Carlo models to correct for hadronization effects. New insights have now been gained in this regard. In a recent calculation of the leading non-perturbative contribution to the C-parameter in the three-jet symmetric limit based on an effective coupling approach it was found that it differs by a factor of two from the one in the two-jet limit [686]. Subsequently, using the same effective coupling approach, the leading $1/Q$ power correction was computed for thrust and the C-parameter in the full three-jet region under the assumption of a large (negative) $n_f$ limit [687]. Expanding on this approach, Ref. [688] extended the calculation of non-perturbative corrections to include the heavy-jet mass, the difference of jet masses, wide broadening, and, with some caveat, the three-jet resolution variable in the Durham algorithm. The key finding in Ref. [688] is that non-perturbative corrections computed in the three-jet region significantly deviate from those computed in the two-jet limit and hence the aforementioned fits based on power corrections in the two-jet limit result in smaller values of $\alpha_s(m_Z^2)$. Another important observation is that the inclusion of resummation effects introduces a relatively substantial ambiguity outside the two-jet limit. Additionally, other factors such as the choice of mass-scheme used to extend the definition of event shapes to massive hadrons can have significant effects.

These findings are inconsistent with the very small experimental, hadronization, and theoretical uncertainties of only 2, 5, and 9 per-mille, respectively, as reported in Refs. [683,685]. For these reasons, we exclude the results of Refs. [683–685] from the average. Determinations based on corrections for non-perturbative hadronization effects using QCD-inspired Monte Carlo generators have also faced criticism due to the differing nature of parton-level simulations compared to fixed-order calculations. However, these determinations typically exhibit a more conservative theoretical uncertainty.

Not included in the computation of the world average but worth mentioning are a computation of the NLO corrections to 5-jet production and comparison to the measured 5-jet rates at LEP [689], giving $\alpha_s(m_Z^2) = 0.1156^{+0.0041}_{-0.0034}$, and a computation of non-perturbative and perturbative QCD contributions to the scale evolution of quark and gluon jet multiplicities, including resummation, resulting in $\alpha_s(m_Z^2) = 0.1199 \pm 0.0026$ [690].

The unweighted average of the considered determinations as shown in the $e^+e^-$ sector of Fig. 9.2 yields $\alpha_s(m_Z^2) = 0.1189 \pm 0.0037$.

### 9.4.5 Observables from hadron-induced collisions:

Until recently, determinations of $\alpha_s$ using hadron collider data, mostly from jet or $t\bar{t}$ production processes, could be performed at NLO only. NNLO calculations have now become available for $t\bar{t}$ [128, 613, 614] and for inclusive jet, dijet, and three-jet production [253, 264, 266, 691, 692]. For $t\bar{t}$ production, in addition, logarithms to NNLL have been resummed [615]. Both should be supplemented by electroweak corrections [693–696], which become important for high-$p_T$ collisions at the LHC. Z+jet production, studied with respect to an $\alpha_s$ determination at NLO from multi-jet events in Ref. [697], is also known at NNLO for the 1-jet case [239, 698].
Determinations of $\alpha_s$ from production cross sections at hadron colliders also require a knowledge of the relevant PDFs for those $\alpha_s$ values. Two strategies are pursued for the extraction of $\alpha_s$, one using pre-determined PDFs as input and a second strategy fitting the proton PDFs together with the strong coupling constant. Each global PDF group produces PDF sets not only for a value of 0.118 that is close to the world average, but also for a wide range above and below in increments of 0.001, which can be used in such determinations. The latter technique of simultaneously fitting $\alpha_s$ and the PDFs is technically more accurate, given that the new data used in the determination of $\alpha_s$ may modify the PDFs in a manner not taken into account by the $\alpha_s$-variation PDFs provided by the fitting collaborations. The former technique may result in a bias of unknown magnitude [699]. As the LHC experiments have the ability to combine a PDF fit with the $\alpha_s$ determination, for example with tools like described in Ref. [700], we expect more experimental joint determinations of PDFs and of $\alpha_s$ for future iterations of this review.

The first determination of $\alpha_s$ at NNLO accuracy in QCD has been reported by CMS from the $t\bar{t}$ production cross section at $\sqrt{s} = 7$ TeV [618]. In former *Reviews* this opened up a new sub-field on its own. In the meantime, multiple datasets on $t\bar{t}$ production from Tevatron at $\sqrt{s} = 1.96$ TeV and from LHC at $\sqrt{s} = 7, 8$, and 13 TeV have been analyzed simultaneously to determine $\alpha_s$ [619] to $\alpha_s(m_Z^2) = 0.1177^{+0.0034}_{-0.0036}$, where the largest uncertainties are associated with missing higher orders and with PDFs, and where an additional complication arises from the top-mass dependence. Since this combined analysis contains among other things an updated measurement as compared to the dataset used by CMS, the latter is replaced in the averaging by this combined result. A second entry into this sub-field is given by an analysis of $t\bar{t}$ production data at $\sqrt{s} = 13$ TeV from the CMS Collaboration [610]. From the four values derived for different PDF sets, the unweighted average is taken: $\alpha_s(m_Z^2) = 0.1145^{+0.0036}_{-0.0031}$. A second analysis by CMS using differential distributions of $t\bar{t}$ production [611] has been performed at NLO only and is not further considered.

From collisions at HERA the $\alpha_s$ determination at NNLO using inclusive jet and dijet measurements in addition to DIS data of the H1 Collaboration [505] has been corrected for an issue with the theory prediction reported in Ref. [220]. We choose the result $\alpha_s(m_Z^2) = 0.1147 \pm 0.0025$ that is derived from a simultaneous fit of the strong coupling constant together with proton PDFs. We note that results of this section derived from such a simultaneous fit will be marked with a “*” in the corresponding figures. A second determination [280] combines multiple datasets on inclusive jet production of the H1 and ZEUS collaborations into one fit using interpolation grids at NNLO. The result of $\alpha_s(m_Z^2) = 0.1178 \pm 0.0026$ has been updated for the same issue reported above [220] and is included in the unweighted average, although some of the inclusive jet data have already been used in the previous analysis. We note that the ZEUS Collaboration has presented preliminary results from a new inclusive jet measurement [701].

A first new determination of $\alpha_s(m_Z^2)$ at NNLO from inclusive jet production at the LHC has been presented by the CMS Collaboration [522]. A simultaneous fit of the strong coupling constant and the proton PDFs to the HERA DIS and the new LHC jet data gives: $\alpha_s(m_Z^2) = 0.1166 \pm 0.0016$. The ATLAS Collaboration has published an extraction of $\alpha_s(m_Z^2)$ at NNLO from the transverse energy-energy correlation (TEEC) and its asymmetry (ATEEC) [459]. The results of this first derivation from an event shape observable requiring three-jet predictions at NNLO are $\alpha_s(m_Z^2) = 0.1175^{+0.0035}_{-0.0018}$ for the TEEC and $\alpha_s(m_Z^2) = 0.1185^{+0.0027}_{-0.0015}$ for the ATEEC, respectively. We include the unweighted average of the two numbers $\alpha_s(m_Z^2) = 0.1180^{+0.0031}_{-0.0017}$ as a new result in this category. Very recently, two preliminary determinations of the strong coupling constant, albeit at NLO, have been reported by CMS from energy correlators inside jets and from azimuthal correlations among jets [702, 703].

Finally, Ref. [704] extracts the strong coupling constant from a fit at NNLO to measurements of inclusive Z and W boson

<figure><img src="image_19.png" /><figcaption>**Figure 9.3:** Summary of determinations of $\alpha_s(m_Z^2)$ at NLO from inclusive jet, dijet, and multi-jet measurements at hadron colliders. The uncertainty is dominated by estimates of the impact of missing higher orders. The yellow (light shaded) bands and dotted lines indicate average values for each sub-field. The dashed line and blue (dark shaded) band represent the final world average value of $\alpha_s(m_Z^2)$. The “*” symbol indicates determinations including a simultaneous fit of PDFs.</figcaption></figure>

production by experiments at the LHC and Tevatron colliders. From the four values quoted for different PDF sets we include the unweighted average for the CT14 and MMHT14 PDFs $\alpha_s(m_Z^2) = 0.1180^{+0.0017}_{-0.0015}$. The other results suffer from either a bad fit quality or $\alpha_s(m_Z^2)$ values outside the range of validity for the PDF set. Furthermore, two determinations using measurements of the Z boson's recoil at low $p_T$ have been submitted using data of the CDF respectively ATLAS experiments [705,706]. The very small size of the estimated uncertainties has stimulated some discussion.

As unweighted pre-average for this sub-field we now obtain $\alpha_s(m_Z^2) = 0.1168 \pm 0.0027$. If the stricter requirement of simultaneous fits with PDFs is imposed, then only the H1 and CMS results are left giving $\alpha_s(m_Z^2) = 0.1157 \pm 0.0021$ for this sub-field. In sect. 9.4.8 below, we will also report the outcome for this choice.

Many further $\alpha_s$ determinations from jet measurements have not yet been advanced to NNLO accuracy. A selection of results from inclusive jet [504, 518, 706–712], dijet [529], and multi-jet measurements [455, 457, 458, 504, 713–717] is presented in Fig. 9.3, where the uncertainty in most cases is dominated by the impact of missing higher orders estimated through scale variations. From the CMS Collaboration we quote for the inclusive jet production at $\sqrt{s} = 7$ and 8 TeV, and for dijet production at 8 TeV the values that have been derived in a simultaneous fit with the PDFs and marked with “*” in the figure. The last point of the inclusive jet sub-field from Ref. [712] is derived from a simultaneous fit to six datasets from different experiments and partially includes data used already for the other data points, e.g. the CMS result at 7 TeV.

The multi-jet $\alpha_s$ determinations are based on three-jet cross sections (m3j), three- to two-jet cross-section ratios (R32), di-jet angular decorrelations (RdR, RdPhi), and the aforementioned transverse energy-energy-correlations, but at $\sqrt{s} = 7$ and 8 TeV only. The H1 result is extracted from a fit to inclusive one-, two-, and three-jet cross sections (nj) simultaneously.

All NLO results are within their large uncertainties in agreement with the world average and the associated analyses provide valuable new values for the scale dependence of $\alpha_s$ at energy scales now extending beyond 2.0 TeV as shown in Fig. 9.5.
### 9.4.6 Electroweak precision fit:

For this category, we take the global electroweak fit of Ref. [718], which includes kinematic top quark and W boson mass measurements from the LHC, determinations of the effective leptonic electroweak mixing angles from the Tevatron, a Higgs mass measurement from ATLAS and CMS, and an evaluation of the hadronic contribution to the running of the electromagnetic coupling at the Z-boson mass. We also consider the fit of $\alpha_s(m_Z^2)$ from the global fit to electroweak data presented in the supplemental material of Ref. [719]. We choose to include the result of their “conservative scenario”<fnref n="13" /> in accounting for the uncertainties of $m_t$ and $m_W$, which avoids a potential bias from inconsistencies between the world average of the W boson mass and the new measurement reported by the CDF Collaboration in Ref. [721]. In addition, we use the newer results of the electroweak fit at the Z mass pole from LEP and SLC data presented in Sec. Electroweak Model and constraints on New Physics of the 2022 edition of this Review. All three determinations, $\alpha_s(m_Z^2) = 0.1194 \pm 0.0029$ [718], $\alpha_s(m_Z^2) = 0.1193 \pm 0.0028$ [719], and $\alpha_s(m_Z^2) = 0.1221 \pm 0.0027$ [630], are also in perfect agreement with the original result obtained from LEP and SLD data [722]. Our pre-averaging gives $\alpha_s(m_Z^2) = 0.1203 \pm 0.0028$.

We note, however, that results from electroweak precision data strongly depend on the strict validity of Standard Model predictions and the existence of the minimal Higgs mechanism to implement electroweak symmetry breaking. Any - even small - deviation of nature from this model could strongly influence this extraction of $\alpha_s$.

### 9.4.7 Lattice QCD:

Several methods exist to extract the strong coupling constant from lattice QCD, as reviewed also in Sec. Lattice QCD of this Review and in Ref. [723]. The Flavour Lattice Averaging Group has considered the most up-to-date determinations and combined them to produce an updated $\alpha_s(m_Z^2)$ estimate [636]. Their final result is obtained by considering a multitude of possible input calculations and by retaining in their final estimate only those that fulfill their predefined quality criteria, detailed in the Sec. 9.2.1 of Ref. [636]. In summary, a determination of $\alpha_s$ needs to satisfy the following requirements:

* The determination of $\alpha_s$ should be based on a comparison of a short-distance quantity $Q$ at scale $\mu$ with a well-defined continuum limit, which does not involve UV or IR divergences, when the quantity is expressed by using an expansion in terms of a short distance definition of the strong coupling (e.g. in the $\overline{MS}$ scheme);

* The scale $\mu$, at which the determination is carried out, has to be sufficiently large such that the error associated with the perturbative truncation remains small enough;

* If $Q$ is defined by physical quantities in infinite volume, one needs to satisfy the constraints $L/a \gg \mu/\Lambda_{QCD}$, where $L$ is the lattice size and $a$ the lattice spacing;

* Only results for $n_f \ge 3$ are considered.

* As is the case for the PDG average, a calculation must be published in a peer-reviewed journal to be eligible for inclusion into the FLAG estimate.

We note that the FLAG criteria applied now are unchanged compared to the FLAG 2019 Review [724] and are considered to be relatively loose. More stringent criteria have already been formulated by FLAG, and it is likely that in future their estimate includes only those results that satisfy these stricter criteria.

Altogether, from a large number of old and new results six from the 2019 Review also enter the final FLAG 2021 estimate [725–729]<fnref n="14" /> and three new ones qualify for inclusion [730–732]. These determinations, together with their uncertainties, are displayed in Fig. 9.4. The yellow (light shaded) band and dotted line indicate the FLAG 2021 estimate, while the dashed line and blue (dark

<figure><img src="image_20.png" /><figcaption>Figure 9.4: Lattice determinations that enter the FLAG 2021 estimate. The yellow (light shaded) band and dotted line indicate the *unweighted average* value for this sub-field. The dashed line and blue (dark shaded) band represent our final world average value of $\alpha_s(m_Z^2)$.</figcaption></figure>

shaded) band represent the world average (see later). The level of agreement of individual results to the world average, or to the non-lattice world average is very similar.

Similarly to what is done here, the FLAG Collaboration built pre-averages of results within the various classes. The five categories that contribute to their final estimate are: step-scaling methods ($\alpha_s(m_Z^2) = 0.1184 \pm 0.00081$), the potential at short distances ($\alpha_s(m_Z^2) = 0.11782 \pm 0.00165$), Wilson loops ($\alpha_s(m_Z^2) = 0.11871 \pm 0.00128$), heavy-quark current two-point functions ($\alpha_s(m_Z^2) = 0.11826 \pm 0.00200$) and, for the first time, light quark vacuum polarization ($\alpha_s(m_Z^2) = 0.11863 \pm 0.00360$). Other categories such as the calculation of QCD vertices, or of the eigenvalue spectrum of the Dirac operator have not yet published results that fulfill all requirements.

We note that, in addition to presenting new results and improvements on previous works, the ALPHA Collaboration [733] has introduced a novel approach to non-perturbative renormalization called decoupling. This strategy shifts the perspective on results involving unphysical flavor numbers, particularly $n_f = 0$. By performing a non-perturbative matching calculation, these results can be non-perturbatively related to results with $n_f > 0$. Consequently, obtaining precise and controlled $n_f = 0$ results becomes of great importance, with significant implications for future FLAG reports.

The final value is obtained by performing a *weighted average* of the pre-averages. The final uncertainty however is not the combined uncertainty of the pre-averages (which is 0.0006), since the errors on almost all determinations are dominated by the perturbative truncation error. Instead, the error on the pre-range for $\alpha_s$ from the step-scaling method is taken, since perturbative truncation errors are sub-dominant in this method. The final FLAG 2021 estimate (rounded to four digits) is

$$ \alpha_s(m_Z^2) = 0.1184 \pm 0.0008 \quad (\text{FLAG 2021 estimate}), \quad (9.23) $$

which is fully compatible with the FLAG 2019 result of $\alpha_s(m_Z^2) = 0.1182 \pm 0.0008$.

We believe that this result expresses to a large extent the consensus of the lattice community and that the imposed criteria and the rigorous assessment of systematic uncertainties qualify for a direct inclusion of this FLAG estimate here. As in the previous review, we therefore adopt the FLAG estimate with its uncertainty as our value of $\alpha_s$ for the lattice category. Moreover, this lattice result will not be directly combined with any other sub-field

<footnote n="13">¹³We note that this result agrees with the previous fit in the “conservative scenario” of Ref. [720], which precedes the latest CDF W mass determination.</footnote>

<footnote n="14">¹⁴Reference [727] contains two results.</footnote>
<table><caption>**Table 9.1:** Unweighted and weighted pre-averages of $\alpha_s(m_Z^2)$ for each sub-field in columns two and three. The bottom line corresponds to the combined result (without lattice gauge theory) using the $\chi^2$ averaging method. The same $\chi^2$ averaging is used for column four combining all unweighted averages except for the sub-field of column one. See text for more details.</caption><thead><tr><th>averages per sub-field</th><th>unweighted</th><th>weighted</th><th>unweighted without subfield</th></tr></thead><tbody><tr><td>τ decays &amp; low Q<sup>2</sup></td><td>0.1173 ± 0.0017</td><td>0.1174 ± 0.0009</td><td>0.1177 ± 0.0013</td></tr><tr><td>QQ bound states</td><td>0.1181 ± 0.0037</td><td>0.1177 ± 0.0011</td><td>0.1175 ± 0.0011</td></tr><tr><td>PDF fits</td><td>0.1161 ± 0.0022</td><td>0.1168 ± 0.0014</td><td>0.1179 ± 0.0011</td></tr><tr><td>e<sup>+</sup>e<sup>−</sup> jets &amp; shapes</td><td>0.1189 ± 0.0037</td><td>0.1187 ± 0.0017</td><td>0.1174 ± 0.0011</td></tr><tr><td>hadron colliders</td><td>0.1168 ± 0.0027</td><td>0.1169 ± 0.0014</td><td>0.1177 ± 0.0011</td></tr><tr><td>electroweak</td><td>0.1203 ± 0.0028</td><td>0.1203 ± 0.0016</td><td>0.1171 ± 0.0011</td></tr><tr><td>PDG 2023 (without lattice)</td><td>0.1175 ± 0.0010</td><td>0.1178 ± 0.0005</td><td>n/a</td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>**Figure 9.5:** Summary of determinations of $\alpha_s$ as a function of the energy scale $Q$ compared to the running of the coupling computed at five loops taking as an input the current PDG average, $\alpha_s(m_Z^2) = 0.1180 \pm 0.0009$. Compared to the previous edition, numerous points have been updated or added.</figcaption></figure>

average, but with our non-lattice average to give our final world average value for $\alpha_s$.

**9.4.8 Determination of the world average value of $\alpha_s(m_Z^2)$:** Obtaining a world average value for $\alpha_s(m_Z^2)$ is a non-trivial exercise. A certain arbitrariness and subjective component is inevitable because of the choice of measurements to be included in the average, the treatment of (non-Gaussian) systematic uncertainties of mostly theoretical nature, as well as the treatment of correlations among the various inputs, of theoretical as well as experimental origin.

We have chosen to determine pre-averages for sub-fields of measurements that are considered to exhibit a maximum degree of independence among each other, considering experimental as well as theoretical issues. The six pre-averages, illustrated also in Fig. 9.2, are listed in column two of Table 9.1. We recall that these are exclusively obtained from extractions that are based on (at least) NNLO QCD predictions, and are published in peer-reviewed journals at the time of completing this Review. To obtain our final world average, we first combine these six pre-averages, excluding the lattice result, using a $\chi^2$ averaging method. This gives

$$ \alpha_s(m_Z^2) = 0.1175 \pm 0.0010 \quad (\text{PDG 2023 without lattice}). \qquad (9.24) $$

This result is fully compatible with the lattice estimate Eq. (9.23) and has a comparable error. To avoid a possible over-reduction, we combine these two numbers using an unweighted average and take as an uncertainty the average between these two uncertainties. This gives our final world average value

$$ \alpha_s(m_Z^2) = 0.1180 \pm 0.0009 \quad (\text{PDG 2023 average}). \qquad (9.25) $$

If for the sub-field of hadron colliders we are more restrictive and instead only accept results from a simultaneous fit of PDFs, we arrive at $0.1157 \pm 0.0021$ for this sub-field leading to $0.1172 \pm$

$0.0010$ (without lattice) and $\alpha_s(m_Z^2) = 0.1178 \pm 0.0009$ for the final average. Both the new world average value and the restricted result are compatible with each other and changed only marginally as compared to the values reported in the last edition of this Review.

It also stands to question whether the sub-fields of PDF fits and hadron colliders still are as independent as originally assumed. To test the potential impact, the fit has been repeated while grouping all $\alpha_s(m_Z^2)$ determinations of both sub-fields into a common one. We obtain $0.1164 \pm 0.0024$ for the new larger sub-field and $0.1178 \pm 0.0011$ for the combination with all other sub-fields except lattice, which is well within the estimated uncertainties. Moreover, we present in column four of Table 9.1 the combined result for $\alpha_s(m_Z^2)$ when the respective sub-field is omitted from the combination. The variation in values obtained for $\alpha_s(m_Z^2)$ is $\pm 0.0004$, which is less than half our estimated uncertainty.

Since two long-standing issues causing offsets among determinations of $\alpha_s(m_Z^2)$ in the sub-fields of $\tau$ decays and low $Q^2$ and $e^+e^-$ jets & shapes have been resolved, it may be argued that a weighted fit among the non-lattice sub-fields may be warranted. We compare the outcome of weighted fits with our standard procedure in columns two and three of Table 9.1. We observe that the weighted averages are rather close to the unweighted ones. However, the uncertainties become significantly smaller. This approach may be too aggressive as it ignores the correlations among the data, methods, and theory ingredients of the various determinations. We feel that the uncertainty of $\pm 0.0005$ is an underestimation of the true error. We also note that in the unweighted combination the estimated uncertainty for each sub-field is larger than the spread of the results as given by the standard deviation. In the weighted fit this crosscheck fails in four out of six cases.

The last several years have seen clarification of some persistent concerns and a wealth of new results at NNLO, providing not only a rather precise and reasonably stable world average value of $\alpha_s(m_Z^2)$, but also a clear signature and proof of the energy dependence of $\alpha_s$ in full agreement with the QCD prediction of asymptotic freedom. This is demonstrated in Fig. 9.5, where results of $\alpha_s(Q^2)$ obtained at discrete energy scales $Q$ are summarized, which now mostly include those based on NNLO QCD.<fnref n="15" /> Thanks to the results from the LHC, the energy scales, at which $\alpha_s$ is determined, now extend even beyond 2 TeV.<fnref n="16" /> The points in this plot are extracted from Refs. [280, 458, 459, 518, 619, 630, 653, 660, 675–677, 679, 709, 715, 716].

In this combination, as in past combinations, we have considered lattice QCD calculations of $\alpha_s$ independently of experimental/phenomenological determinations. In the future, when the lattice continuum extrapolations are under better control, it may be useful to group lattice QCD determinations of $\alpha_s$ with experimental determinations of $\alpha_s$ that have systematics of similar

<footnote n="15">¹⁵The uncertainties of the HERA jets points were, by mistake, shown at only half their size. The uncertainties of the ALEPH points of the $e^+e^-$ jets/shapes sub-field now correspond to the total uncertainty.</footnote>

<footnote n="16">¹⁶We note, however, that the relevant energy scale of a measurement is not uniquely defined. In addition to being multiplied by factors of e.g. 1/2 or 2, for instance in studies of the ratio of three- to two-jet cross sections at the LHC, the relevant scale can be taken to be the average of the transverse momenta of the two leading jets [714], or alternatively might be chosen to be the transverse momentum of the 3rd jet.</footnote>
origin, in a similar manner as we currently group, for example,
hadron collider results together [723].

References

[1] R. K. Ellis, W. J. Stirling and B. R. Webber, Camb. Monogr.
Part. Phys. Nucl. Phys. Cosmol. **8**, 1 (1996).

[2] J. Campbell, J. Huston, F. Krauss "The Black Book
of Quantum Chromodynamics, a Primer for the QCD
Era," Oxford University Press, UK (2017).

[3] C. A. Baker et al., Phys. Rev. Lett. **97**, 131801 (2006), [hep-ex/0602020].

[4] J. M. Pendlebury et al., Phys. Rev. D92, 9, 092003 (2015), [arXiv:1509.04411].

[5] B. Graner et al., Phys. Rev. Lett. 116, 16, 161601 (2016), [Erratum: Phys. Rev. Lett.119,no.11,119901(2017)], [arXiv:1601.04339].

[6] https://www PSI .ch/en/nedm/edms-world-wide.

[7] J. E. Kim and G. Carosi, Rev. Mod. Phys. **82**, 557 (2010),
[arXiv:0807.3125].

[8] G. Dissertori, I. G. Knowles and M. Schmelling, *High energy experiments and theory*, Oxford, UK: Clarendon (2003).

[9] R. Brock et al. (CTEQ), Rev. Mod. Phys. **67**, 157 (1995).

[10] K. Melnikov, CERN Yellow Rep. School Proc. 3, 37 (2018).

[11] F. Gross et al. (2022), [arXiv:2212.11107].

[12] T. van Ritbergen, J. A. M. Vermaseren and S. A. Larin,
Phys. Lett. B400, 379 (1997), [hep-ph/9701390].

[13] M. Czakon, Nucl. Phys. B710, 485 (2005), [hep-ph/0411261].

[14] P. A. Baikov, K. G. Chetyrkin and J. H. Kühn, Phys. Rev. Lett. 118, 8, 082002 (2017), [arXiv:1606.08659].

[15] T. Luthe et al., JHEP 07, 127 (2016), [arXiv:1606.08662].

[16] F. Herzog et al., JHEP 02, 090 (2017), [arXiv:1701.01404].

[17] T. Luthe et al., JHEP 10, 166 (2017), [arXiv:1709.07718].

[18] K. G. Chetyrkin et al., JHEP 10, 179 (2017), [Addendum: JHEP12,006(2017)], [arXiv:1709.08541].

[19] W. A. Bardeen et al., Phys. Rev. D18, 3998 (1978).

[20] D. J. Gross and F. Wilczek, Phys. Rev. Lett. 30, 1343
(1973), [,271(1973)].

[21] H. D. Politzer, Phys. Rev. Lett. 30, 1346 (1973), [,274(1973)].

[22] Y. Schroder and M. Steinhauser, JHEP 01, 051 (2006),
[hep-ph/0512058].

[23] K. G. Chetyrkin, J. H. Kuhn and C. Sturm, Nucl. Phys. B744, 121 (2006), [hep-ph/0512060].

[24] A. G. Grozin et al., JHEP 09, 066 (2011), [arXiv:1107.5970].

[25] M. Dalla Brida et al. (ALPHA), Phys. Rev. Lett. 117, 18, 182001 (2016), [arXiv:1604.06193].

[26] M. Beneke, Phys. Rept. 317, 1 (1999), [hep-ph/9807443].

[27] M. Beneke et al., Phys. Lett. B775, 63 (2017), [arXiv:1605.03609].

[28] A. H. Hoang, C. Lepenik and M. Preisser, JHEP 09, 099 (2017), [arXiv:1706.08526].

[29] P. Marquard et al., Phys. Rev. Lett. 114, 14, 142002 (2015), [arXiv:1502.01030].

[30] P. A. Baikov et al., Phys. Lett. B714, 62 (2012), [arXiv:1206.1288].

[31] K. G. Chetyrkin, J. H. Kuhn and A. Kwiatkowski (1996),
[Phys. Rept.277,189(1996)], [hep-ph/9503396].

[32] Y. Kiyo et al., Nucl. Phys. B823, 269 (2009), [arXiv:0907.2120].

[33] P. A. Baikov et al., Phys. Rev. Lett. 108, 222003 (2012), [arXiv:1201.5804].

[34] P. A. Baikov, K. G. Chetyrkin and J. H. Kuhn, Phys. Rev. Lett. 101, 012002 (2008), [arXiv:0801.1821].

[35] F. Herzog et al., JHEP 08, 113 (2017), [arXiv:1707.01044].

[36] V. A. Novikov et al., Nucl. Phys. B174, 378 (1980).

[37] T.-J. Hou et al., Phys. Rev. D 103, 1, 014013 (2021),
[arXiv:1912.10053].

[38] R. D. Ball et al. (PDF4LHC Working Group), J. Phys. G 49, 8, 080501 (2022), [arXiv:2203.05506].

[39] R. Abdul Khalek et al. (NNPDF), Eur. Phys. J. C 79:838 (2019), [arXiv:1905.04311].

[40] R. Abdul Khalek et al. (NNPDF), Eur. Phys. J. C 79, 11, 931 (2019), [arXiv:1906.10698].

[41] R. D. Ball and R. L. Pearson, Eur. Phys. J. C 81, 9, 830 (2021), [arXiv:2105.05114].

[42] Z. Kassabov, M. Ubiali and C. Voisey, JHEP 03, 148 (2023), [arXiv:2207.07616].

[43] L. A. Harland-Lang and R. S. Thorne, Eur. Phys. J. C 79, 3, 225 (2019), [arXiv:1811.08434].

[44] J. McGowan et al., Eur. Phys. J. C 83, 3, 185 (2023), [Erratum: Eur.Phys.J.C 83, 302 (2023)], [arXiv:2207.04739].

[45] X. Ji, Phys. Rev. Lett. 110, 262002 (2013), [arXiv:1305.1539].

[46] G. C. Rossi and M. Testa, Phys. Rev. D96, 1, 014507 (2017), [arXiv:1706.04428].

[47] G. Rossi and M. Testa, Phys. Rev. D 98, 5, 054028 (2018), [arXiv:1806.00808].

[48] L. Del Debbio, T. Giani and C. J. Monahan, JHEP 09, 021 (2020), [arXiv:2007.02131].

[49] H.-W. Lin et al., Phys. Rev. D91, 054510 (2015), [arXiv:1402.1462]; C. Alexandrou et al., Phys. Rev. D92, 014502 (2015), [arXiv:1504.07455].

[50] H.-W. Lin et al., Prog. Part. Nucl. Phys. 100, 107 (2018), [arXiv:1711.07916].

[51] C. Alexandrou et al., Phys. Rev. Lett. 121, 11, 112001 (2018), [arXiv:1803.02685].

[52] J.-W. Chen et al. (2018), [arXiv:1803.04393].

[53] K. Cichy, L. Del Debbio and T. Giani, JHEP 10, 137 (2019), [arXiv:1907.06037].

[54] K. Cichy and M. Constantinou, Adv. High Energy Phys. 2019, 3036904 (2019), [arXiv:1811.07248].

[55] M. A. Ebert, I. W. Stewart and Y. Zhao, JHEP 09, 037 (2019), [arXiv:1901.03685].

[56] M. Constantinou, Eur. Phys. J. A 57, 2, 77 (2021), [arXiv:2010.02445].

[57] M. Constantinou et al., Prog. Part. Nucl. Phys. 121, 103908 (2021), [arXiv:2006.08636].

[58] L. Del Debbio et al., JHEP 02, 138 (2021), [arXiv:2010.03996].

[59] Z. Fan, R. Zhang and H.-W. Lin, Int. J. Mod. Phys. A 36, 13, 2150080 (2021), [arXiv:2007.16113].

[60] T. Khan et al. (HadStruc), Phys. Rev. D 104, 9, 094516 (2021), [arXiv:2107.08960].

[61] Z. Fan, W. Good and H.-W. Lin (2022), [arXiv:2210.09985].

[62] J. Gao, L. Harland-Lang and J. Rojo, Phys. Rept. 742, 1 (2018), [arXiv:1709.04922].

[63] K. Kovařík, P. M. Nadolsky and D. E. Soper, Rev. Mod. Phys. 92, 4, 045003 (2020), [arXiv:1905.06957].

[64] T.-J. Hou et al., Phys. Rev. D 107, 7, 076018 (2023), [arXiv:2211.11064].

[65] J. D. Bjorken and E. A. Paschos, Phys. Rev. 185, 1975 (1969).

[66] J. A. M. Vermaseren, A. Vogt and S. Moch, Nucl. Phys. B724, 3 (2005), [hep-ph/0504242].

[67] S. Moch, J. A. M. Vermaseren and A. Vogt, Nucl. Phys. B813, 220 (2009), [arXiv:0812.4168].
[68] J. Davies et al., PoS DIS2016, 059 (2016), [arXiv:1606.08907].

[69] E. Laenen et al., Nucl. Phys. B392, 162 (1993); S. Riemersma, J. Smith and W. L. van Neerven, Phys. Lett. B347, 143 (1995), [hep-ph/9411431].

[70] J. Blümlein et al. (2019), [arXiv:1903.06155].

[71] J. Blümlein (2023), [arXiv:2306.01362].

[72] A. Manohar et al., Phys. Rev. Lett. 117, 24, 242002 (2016), [arXiv:1607.04266].

[73] A. V. Manohar et al., JHEP 12, 046 (2017), [arXiv:1708.01256].

[74] L. Buonocore et al., JHEP 08, 08, 019 (2020), [arXiv:2005.06477].

[75] J. C. Collins, D. E. Soper and G. F. Sterman, Adv. Ser. Direct. High Energy Phys. 5, 1 (1989), [hep-ph/0409313].

[76] J. Collins and J.-W. Qiu, Phys. Rev. D 75, 114014 (2007), [arXiv:0705.2141].

[77] G. Sterman, in "Snowmass 2021," (2022), [arXiv:2207.06507].

[78] R. Boussarie et al. (2023), [arXiv:2304.03302].

[79] V. N. Gribov and L. N. Lipatov, Sov. J. Nucl. Phys. 15, 438 (1972), [Yad. Fiz.15,781(1972)]; L. N. Lipatov, Sov. J. Nucl. Phys. 20, 94 (1975), [Yad. Fiz.20,181(1974)]; G. Altarelli and G. Parisi, Nucl. Phys. B126, 298 (1977); Y. L. Dokshitzer, Sov. Phys. JETP 46, 641 (1977), [Zh. Eksp. Teor. Fiz.73,1216(1977)].

[80] G. Curci, W. Furmanski and R. Petronzio, Nucl. Phys. B175, 27 (1980); W. Furmanski and R. Petronzio, Phys. Lett. 97B, 437 (1980).

[81] A. Vogt, S. Moch and J. A. M. Vermaseren, Nucl. Phys. B691, 129 (2004), [hep-ph/0404111]; S. Moch, J. A. M. Vermaseren and A. Vogt, Nucl. Phys. B688, 101 (2004), [hep-ph/0403192].

[82] J. Davies et al., Nucl. Phys. B 915, 335 (2017), [arXiv:1610.07477].

[83] S. Moch et al., JHEP 10, 041 (2017), [arXiv:1707.08315].

[84] A. Vogt et al., PoS RADCOR2017, 046 (2018), [arXiv:1801.06085].

[85] A. Vogt et al., PoS LL2018, 050 (2018), [arXiv:1808.08981].

[86] G. Falcioni et al., Phys. Lett. B 842, 137944 (2023), [arXiv:2302.07593].

[87] S. Moch, J. A. M. Vermaseren and A. Vogt, Nucl. Phys. B889, 351 (2014), [arXiv:1409.5131].

[88] D. de Florian, G. F. R. Sborlini and G. Rodrigo, Eur. Phys. J. C76, 5, 282 (2016), [arXiv:1512.00612]; D. de Florian, G. F. R. Sborlini and G. Rodrigo, JHEP 10, 056 (2016), [arXiv:1606.02887].

[89] D. de Florian and L. P. Conte (2023), [arXiv:2305.14144].

[90] R. S. Thorne, Phys. Rev. D73, 054019 (2006), [hep-ph/0601245].

[91] S. Forte et al., Nucl. Phys. B834, 116 (2010), [arXiv:1001.2312].

[92] M. Guzzi et al., Phys. Rev. D86, 053005 (2012), [arXiv:1108.5112].

[93] V. S. Fadin, E. A. Kuraev and L. N. Lipatov, Phys. Lett. 60B, 50 (1975).

[94] I. I. Balitsky and L. N. Lipatov, Sov. J. Nucl. Phys. 28, 822 (1978), [Yad. Fiz.28,1597(1978)].

[95] R. D. Ball et al., Eur. Phys. J. C78, 4, 321 (2018), [arXiv:1710.05935].

[96] H. Abramowicz et al. (H1, ZEUS), Eur. Phys. J. C75, 12, 580 (2015), [arXiv:1506.06042].

[97] L. N. Lipatov, Sov. J. Nucl. Phys. 23, 338 (1976), [Yad. Fiz.23,642(1976)].

[98] E. A. Kuraev, L. N. Lipatov and V. S. Fadin, Sov. Phys. JETP 45, 199 (1977), [Zh. Eksp. Teor. Fiz.72,377(1977)].

[99] V. S. Fadin and L. N. Lipatov, Phys. Lett. B429, 127 (1998), [hep-ph/9802290].

[100] M. Ciafaloni and G. Camici, Phys. Lett. B430, 349 (1998), [hep-ph/9803389].

[101] G. Altarelli, R. D. Ball and S. Forte, Nucl. Phys. B799, 199 (2008), [arXiv:0802.0032].

[102] M. Ciafaloni et al., JHEP 08, 046 (2007), [arXiv:0707.1453].

[103] C. D. White and R. S. Thorne, Phys. Rev. D75, 034005 (2007), [hep-ph/0611204].

[104] E. Iancu et al., Phys. Lett. B744, 293 (2015), [arXiv:1502.05642].

[105] N. Gromov, F. Levkovich-Maslyuk and G. Sizov, Phys. Rev. Lett. 115, 25, 251601 (2015), [arXiv:1507.04010]; V. N. Velizhanin (2015), [arXiv:1508.02857]; S. Caron-Huot and M. Herranen, JHEP 02, 058 (2018), [arXiv:1604.07417]; V. Fadin, Chapter 4: BFKL — Past and Future, 63–90 (2021), [arXiv:2012.11931].

[106] I. Balitsky, Nucl. Phys. B463, 99 (1996), [hep-ph/9509348].

[107] Y. V. Kovchegov, Phys. Rev. D60, 034008 (1999), [hep-ph/9901281].

[108] M. Davier et al., Eur. Phys. J. C 56, 305 (2008), [arXiv:0803.0979].

[109] J. C. Collins, D. E. Soper and G. F. Sterman, Nucl. Phys. B261, 104 (1985).

[110] C. Anastasiou et al., Phys. Rev. Lett. 114, 212001 (2015), [arXiv:1503.06056]; C. Anastasiou et al., JHEP 05, 058 (2016), [arXiv:1602.00695].

[111] B. Mistlberger, JHEP 05, 028 (2018), [arXiv:1802.00833].

[112] L.-B. Chen et al., Phys. Lett. B 803, 135292 (2020), [arXiv:1909.06808].

[113] F. Dulat, B. Mistlberger and A. Pelloni, Phys. Rev. D99, 3, 034004 (2019), [arXiv:1810.09462].

[114] L. Cieri et al., JHEP 02, 096 (2019), [arXiv:1807.11501].

[115] X. Chen et al. (2021), [arXiv:2102.07607].

[116] J. Baglio et al., JHEP 12, 066 (2022), [arXiv:2209.06138].

[117] F. A. Dreyer and A. Karlberg, Phys. Rev. Lett. 117, 7, 072001 (2016), [arXiv:1606.00840].

[118] F. A. Dreyer and A. Karlberg, Phys. Rev. D 98, 11, 114016 (2018), [arXiv:1811.07906].

[119] C. Duhr, F. Dulat and B. Mistlberger, Phys. Rev. Lett. 125, 17, 172001 (2020), [arXiv:2001.07717].

[120] X. Chen et al. (2021), [arXiv:2107.09085].

[121] C. Duhr, F. Dulat and B. Mistlberger, JHEP 11, 143 (2020), [arXiv:2007.13313].

[122] M. Bonvini et al., JHEP 08, 105 (2016), [arXiv:1603.08000].

[123] F. Dulat, A. Lazopoulos and B. Mistlberger, Comput. Phys. Commun. 233, 243 (2018), [arXiv:1802.00827].

[124] R. V. Harlander, S. Liebler and H. Mantler, Comput. Phys. Commun. 212, 239 (2017), [arXiv:1605.03190].

[125] O. Brein, A. Djouadi and R. Harlander, Phys. Lett. B579, 149 (2004), [hep-ph/0307206].

[126] D. de Florian and J. Mazzitelli, Phys. Rev. Lett. 111, 201801 (2013), [arXiv:1309.6594].

[127] S. Borowka et al., Phys. Rev. Lett. 117, 1, 012001 (2016), [Erratum: Phys. Rev. Lett.117,no.7,079901(2016)], [arXiv:1604.06447].

[128] M. Czakon, P. Fiedler and A. Mitov, Phys. Rev. Lett. 110, 252004 (2013), [arXiv:1303.6254].

[129] S. Catani et al., Phys. Rev. D99, 5, 051501 (2019), [arXiv:1901.04005].

[130] S. Catani et al., JHEP 03, 029 (2021), [arXiv:2010.11906].
[131] T. Gehrmann et al., Phys. Rev. Lett. 113, 21, 212001 (2014), [arXiv:1408.5243].

[132] F. Cascioli et al., Phys. Lett. B735, 311 (2014), [arXiv:1405.2219].

[133] M. Grazzini, S. Kallweit and M. Wiesemann, Eur. Phys. J. C78, 7, 537 (2018), [arXiv:1711.06631].

[134] G. Heinrich, Phys. Rept. 922, 1 (2021), [arXiv:2009.00516].

[135] A. Huss et al., J. Phys. G 50, 4, 043001 (2023), [arXiv:2207.02122].

[136] M. Greco and A. Vicini, Nucl. Phys. B415, 386 (1994).

[137] A. Hebecker, Phys. Rept. 331, 1 (2000), [hep-ph/9905226].

[138] A. V. Belitsky and A. V. Radyushkin, Phys. Rept. 418, 1 (2005), [hep-ph/0504030].

[139] E. Boos et al. (CompHEP), Nucl. Instrum. Meth. A534, 250 (2004), [hep-ph/0403113].

[140] J. Alwall et al., JHEP 07, 079 (2014), [arXiv:1405.0301]; https://launchpad.net/mg5amcnlo.

[141] M. L. Mangano et al., JHEP 07, 001 (2003), [hep-ph/0206293]; http://cern.ch/mlm/alpgen/.

[142] T. Gleisberg and S. Hoeche, JHEP 12, 039 (2008), [arXiv:0808.3674]; https://sherpa-team.gitlab.io/.

[143] A. Cafarella, C. G. Papadopoulos and M. Worek, Comput. Phys. Commun. 180, 1941 (2009), [arXiv:0710.2427]; http://cern.ch/helac-phegas/.

[144] F. A. Berends and W. T. Giele, Nucl. Phys. B306, 759 (1988).

[145] L. J. Dixon, in "QCD and beyond. Proceedings, Theoretical Advanced Study Institute in Elementary Particle Physics, TASI-95, Boulder, USA, June 4-30, 1995," 539-584 (1996), [hep-ph/9601359], URL http://www-public.slac.stanford.edu/sciDoc/docMeta.aspx?slacPubNumber=SLAC-PUB-7106.

[146] R. Britto, F. Cachazo and B. Feng, Nucl. Phys. B715, 499 (2005), [hep-th/0412308].

[147] F. Cachazo, P. Svrcek and E. Witten, JHEP 09, 006 (2004), [hep-th/0403047].

[148] S. Catani and M. H. Seymour, Nucl. Phys. B485, 291 (1997), [Erratum: Nucl. Phys.B510,503(1998)], [hep-ph/9605323].

[149] S. Frixione, Z. Kunszt and A. Signer, Nucl. Phys. B467, 399 (1996), [hep-ph/9512328].

[150] D. A. Kosower, Phys. Rev. D57, 5410 (1998), [hep-ph/9710213]; J. M. Campbell, M. A. Cullen and E. W. N. Glover, Eur. Phys. J. C9, 245 (1999), [hep-ph/9809429]; D. A. Kosower, Phys. Rev. D71, 045016 (2005), [hep-ph/0311272].

[151] G. Ossola, C. G. Papadopoulos and R. Pittau, Nucl. Phys. B763, 147 (2007), [hep-ph/0609007].

[152] R. Britto, F. Cachazo and B. Feng, Nucl. Phys. B725, 275 (2005), [hep-th/0412103].

[153] R. K. Ellis et al., Nucl. Phys. B822, 270 (2009), [arXiv:0806.3467].

[154] C. F. Berger and D. Forde, Ann. Rev. Nucl. Part. Sci. 60, 181 (2010), [arXiv:0912.3534].

[155] F. Cascioli, P. Maierhofer and S. Pozzorini, Phys. Rev. Lett. 108, 111601 (2012), [arXiv:1111.5206].

[156] Z. Bern, L. J. Dixon and D. A. Kosower, Annals Phys. 322, 1587 (2007), [arXiv:0704.2798].

[157] R. K. Ellis et al., Phys. Rept. 518, 141 (2012), [arXiv:1105.4319].

[158] G. Bevilacqua et al., Comput. Phys. Commun. 184, 986 (2013), [arXiv:1110.1499]; http://cern.ch/helac-phegas/.

[159] G. Cullen et al., Eur. Phys. J. C74, 8, 3001 (2014), [arXiv:1404.7096]; http://gosam.hepforge.org/.

[160] S. Badger et al., Comput. Phys. Commun. 184, 1981 (2013), [arXiv:1209.0100]; https://bitbucket.org/njet/wiki/Home/.

[161] F. Buccioni, S. Pozzorini and M. Zoller, Eur. Phys. J. C78, 1, 70 (2018), [arXiv:1710.11452]; https://openloops.hepforge.org/.

[162] S. Actis et al., Comput. Phys. Commun. 214, 140 (2017), [arXiv:1605.01090].

[163] Z. Nagy, Phys. Rev. D68, 094002 (2003), [hep-ph/0307268]; http://www.desy.de/~znagy/Site/NLOJet++.html.

[164] J. M. Campbell and R. K. Ellis, Phys. Rev. D62, 114012 (2000), [hep-ph/0006304].

[165] J. Baglio et al. (2011), [arXiv:1107.4038]; http://www-itp.particle.uni-karlsruhe.de/~vbfnloweb.

[166] T. Binoth et al., Eur. Phys. J. C16, 311 (2000), [hep-ph/9911340]; http://lapth.in2p3.fr/PHOX FAMILY/.

[167] Z. Bern et al., PoS LL2012, 018 (2012), [arXiv:1210.6684].

[168] G. Cullen, N. Greiner and G. Heinrich, Eur. Phys. J. C73, 4, 2388 (2013), [arXiv:1212.5154].

[169] S. Kallweit et al., JHEP 04, 012 (2015), [arXiv:1412.5157].

[170] A. Denner et al., JHEP 04, 018 (2015), [arXiv:1412.7421].

[171] M. Chiesa, N. Greiner and F. Tramontano, J. Phys. G43, 1, 013002 (2016), [arXiv:1507.08579].

[172] S. Frixione et al., JHEP 06, 184 (2015), [arXiv:1504.03446].

[173] B. Biedermann et al., Eur. Phys. J. C77, 492 (2017), [arXiv:1704.05783].

[174] R. Frederix et al., JHEP 07, 185 (2018), [Erratum: JHEP 11, 085 (2021)], [arXiv:1804.10017].

[175] C. Anastasiou, R. Boughezal and F. Petriello, JHEP 04, 003 (2009), [arXiv:0811.3458].

[176] S. Dittmaier, A. Huss and C. Schwinn, Nucl. Phys. B885, 318 (2014), [arXiv:1403.3216].

[177] S. Dittmaier, A. Huss and C. Schwinn, Nucl. Phys. B904, 216 (2016), [arXiv:1511.08016].

[178] D. de Florian, M. Der and I. Fabre, Phys. Rev. D98, 9, 094008 (2018), [arXiv:1805.12214].

[179] M. Bonetti, K. Melnikov and L. Tancredi, Phys. Rev. D97, 3, 034004 (2018), [arXiv:1711.11113].

[180] M. Bonetti, K. Melnikov and L. Tancredi, Phys. Rev. D97, 5, 056017 (2018), [Erratum: Phys. Rev.D97,no.9,099906(2018)], [arXiv:1801.10403].

[181] C. Anastasiou et al., JHEP 03, 162 (2019), [arXiv:1811.1211].

[182] V. Hirschi, S. Lionetti and A. Schweitzer, JHEP 05, 002 (2019), [arXiv:1902.10167].

[183] M. Becchetti et al., Phys. Rev. D 103, 5, 054037 (2021), [arXiv:2010.09451].

[184] M. Bonetti et al., JHEP 11, 045 (2020), [arXiv:2007.09813].

[185] A. Behring et al., Phys. Rev. D 103, 1, 013008 (2021), [arXiv:2009.10386].

[186] S. Dittmaier, T. Schmidt and J. Schwarz, JHEP 12, 201 (2020), [arXiv:2009.02229].

[187] F. Buccioni et al., Phys. Lett. B 811, 135969 (2020), [arXiv:2005.10221].

[188] L. Buonocore et al., Phys. Rev. D 103, 114012 (2021), [arXiv:2102.12539].

[189] R. Bonciani et al., Phys. Rev. Lett. 128, 1, 012002 (2022), [arXiv:2106.11953].

[190] T. Armadillo et al., JHEP 05, 072 (2022), [arXiv:2201.01754].

[191] F. Buccioni et al., JHEP 06, 022 (2022), [arXiv:2203.11237].

[192] A. Denner and S. Dittmaier, Phys. Rept. 864, 1 (2020), [arXiv:1912.06823].
[193] Z. Bern et al., Nucl. Phys. **B425**, 217 (1994), [hep-ph/9403226].

[194] J. M. Campbell and E. W. N. Glover, Nucl. Phys. **B527**, 264 (1998), [hep-ph/9710255].

[195] S. Catani and M. Grazzini, Phys. Lett. **B446**, 143 (1999), [hep-ph/9810389].

[196] T. Binoth and G. Heinrich, Nucl. Phys. **B585**, 741 (2000), [hep-ph/0004013].

[197] C. Anastasiou, K. Melnikov and F. Petriello, Phys. Rev. **D69**, 076010 (2004), [hep-ph/0311311].

[198] A. Gehrmann-De Ridder, T. Gehrmann and E. W. N. Glover, JHEP **09**, 056 (2005), [hep-ph/0505111].

[199] G. Somogyi, Z. Trocsanyi and V. Del Duca, JHEP **01**, 070 (2007), [hep-ph/0609042].

[200] M. Czakon, Phys. Lett. **B693**, 259 (2010), [arXiv:1005.0274].

[201] F. Caola, K. Melnikov and R. Röntsch, Eur. Phys. J. C **77**, 4, 248 (2017), [arXiv:1702.01352].

[202] L. Magnea et al., JHEP **12**, 107 (2018), [Erratum: JHEP 06, 013 (2019)], [arXiv:1806.09570].

[203] R. M. Prisco and F. Tramontano, JHEP **06**, 089 (2021), [arXiv:2012.05012].

[204] G. Bertolotti et al. (2022), [arXiv:2212.11190].

[205] C. Anastasiou and G. Sterman, JHEP **05**, 242 (2023), [arXiv:2212.12162].

[206] S. Catani and M. Grazzini, Phys. Rev. Lett. **98**, 222002 (2007), [hep-ph/0703012]; http://theory.fi.infn.it/grazzini/codes.html.

[207] R. Bonciani et al., Eur. Phys. J. C **75**, 12, 581 (2015), [arXiv:1508.03585].

[208] R. Boughezal et al., Phys. Rev. Lett. **115**, 6, 062002 (2015), [arXiv:1504.02131].

[209] R. Boughezal, X. Liu and F. Petriello, Phys. Rev. D **91**, 9, 094035 (2015), [arXiv:1504.02540].

[210] J. Gaunt et al., JHEP **09**, 058 (2015), [arXiv:1505.04794].

[211] M. Cacciari et al., Phys. Rev. Lett. **115**, 8, 082002 (2015), [Erratum: Phys. Rev. Lett.120.no.13,139901(2018)], [arXiv:1506.02660].

[212] Z. Capatti, V. Hirschi and B. Ruijl, JHEP **10**, 120 (2022), [arXiv:2203.11038].

[213] L. Cieri et al., JHEP **09**, 155 (2020), [arXiv:2005.01315].

[214] W. J. Torres Bobadilla et al., Eur. Phys. J. C **81**, 3, 250 (2021), [arXiv:2012.02567].

[215] A. Gehrmann-De Ridder et al., Phys. Rev. Lett. **99**, 132002 (2007), [arXiv:0707.1285]; A. Gehrmann-De Ridder et al., JHEP **12**, 094 (2007), [arXiv:0711.4711]; A. Gehrmann-De Ridder et al., Phys. Rev. Lett. **100**, 172001 (2008), [arXiv:0802.0813].

[216] A. Gehrmann-De Ridder et al., Comput. Phys. Commun. **185**, 3331 (2014), [arXiv:1402.4140]; https://eerad3.hepforge.org/.

[217] S. Weinzierl, Phys. Rev. Lett. **101**, 162001 (2008), [arXiv:0807.3241]; S. Weinzierl, JHEP **06**, 041 (2009), [arXiv:0904.1077].

[218] V. Del Duca et al., Phys. Rev. D **94**, 7, 074019 (2016), [arXiv:1606.03453]; V. Del Duca et al., Phys. Rev. Lett. **117**, 15, 152004 (2016), [arXiv:1603.08927].

[219] J. Currie, T. Gehrmann and J. Niehues, Phys. Rev. Lett. **117**, 4, 042001 (2016), [arXiv:1606.03991].

[220] J. Currie et al., JHEP **07**, 018 (2017), [Erratum: JHEP 12, 042 (2020)], [arXiv:1703.05977].

[221] J. Currie et al., JHEP **05**, 209 (2018), [arXiv:1803.09973].

[222] T. Gehrmann et al., Phys. Lett. **B792**, 182 (2019), [arXiv:1812.06104].

[223] K. Melnikov and F. Petriello, Phys. Rev. **D74**, 114017 (2006), [hep-ph/0609070]; http://gate.hep.anl.gov/fpetriello/FEWZ.html.

[224] S. Catani et al., Phys. Rev. Lett. **103**, 082001 (2009), [arXiv:0903.2120]; http://theory.fi.infn.it/gazzani/dy.html.

[225] C. Anastasiou, K. Melnikov and F. Petriello, Nucl. Phys. **B724**, 197 (2005), [hep-ph/0501130]; http://www.phys.ethz.ch/~pheno/fehipro/.

[226] M. Czakon et al. (2021), [arXiv:2105.04436].

[227] M. Grazzini et al., JHEP **05**, 139 (2017), [arXiv:1703.09065].

[228] S. Catani et al., Phys. Rev. Lett. **108**, 072001 (2012), [Erratum: Phys. Rev. Lett.117,no.8,089901(2016)], [arXiv:1110.2375].

[229] J. M. Campbell et al., JHEP **07**, 148 (2016), [arXiv:1603.02663].

[230] M. Grazzini et al., Phys. Lett. **B731**, 204 (2014), [arXiv:1309.7000].

[231] M. Grazzini, S. Kallweit and D. Rathlev, JHEP **07**, 085 (2015), [arXiv:1504.01330].

[232] J. M. Campbell and R. K. Ellis, Phys. Rev. D **60**, 113006 (1999), [hep-ph/9905386].

[233] R. Boughezal et al., Eur. Phys. J. **C77**, 1, 7 (2017), [arXiv:1605.08011].

[234] J. M. Campbell, R. K. Ellis and C. Williams, Phys. Rev. Lett. **118**, 22, 222001 (2017), [arXiv:1612.04333].

[235] X. Chen et al., Submitted to: J. High Energy Phys. (2019), [arXiv:1904.01044].

[236] J. M. Campbell, R. K. Ellis and C. Williams, Phys. Rev. D **96**, 1, 014037 (2017), [arXiv:1703.10109].

[237] A. Gehrmann-De Ridder et al., Phys. Rev. Lett. **120**, 12, 122001 (2018), [arXiv:1712.07543].

[238] A. Gehrmann-De Ridder et al., Phys. Rev. Lett. **117**, 2, 022001 (2016), [arXiv:1507.02850].

[239] R. Boughezal et al., Phys. Rev. Lett. **116**, 15, 152001 (2016), [arXiv:1512.01291].

[240] R. Boughezal et al., Phys. Rev. Lett. **115**, 8, 082003 (2015), [arXiv:1504.07922].

[241] R. Boughezal et al., Phys. Lett. **B748**, 5 (2015), [arXiv:1505.03893].

[242] F. Caola, K. Melnikov and M. Schulze, Phys. Rev. **D92**, 7, 074032 (2015), [arXiv:1508.02684].

[243] X. Chen et al., JHEP **10**, 066 (2016), [arXiv:1607.08817].

[244] G. Ferrera, M. Grazzini and F. Tramontano, Phys. Rev. Lett. **107**, 152003 (2011), [arXiv:1107.1164].

[245] G. Ferrera, M. Grazzini and F. Tramontano, Phys. Lett. **B740**, 51 (2015), [arXiv:1407.4747].

[246] Z. L. Liu and J. Gao, Phys. Rev. D **98**, 7, 071501 (2018), [arXiv:1807.03835].

[247] M. Brucherseifer, F. Caola and K. Melnikov, Phys. Lett. **B736**, 58 (2014), [arXiv:1404.7116].

[248] E. L. Berger et al., Phys. Rev. **D94**, 7, 071501 (2016), [arXiv:1606.08463].

[249] E. L. Berger, J. Gao and H. X. Zhu, JHEP **11**, 158 (2017), [arXiv:1708.09405].

[250] J. Campbell, T. Neumann and Z. Sullivan, JHEP **02**, 040 (2021), [arXiv:2012.01574].

[251] M. Czakon, P. Fiedler and A. Mitov, Phys. Rev. Lett. **115**, 5, 052001 (2015), [arXiv:1411.3007].

[252] M. Czakon et al., JHEP **05**, 034 (2016), [arXiv:1601.05375].

[253] J. Currie, E. W. N. Glover and J. Pires, Phys. Rev. Lett. **118**, 7, 072002 (2017), [arXiv:1611.01460].

[254] X. Chen et al., JHEP **09**, 025 (2022), [arXiv:2204.10173].
[255] M. Czakon et al. (2020), [arXiv:2011.01011].

[256] M. Czakon et al., JHEP 02, 241 (2023), [arXiv:2212.00467].

[257] M. L. Czakon et al. (2021), [arXiv:2102.08267].

[258] J. Cruz-Martinez et al., Phys. Lett. B781, 672 (2018), [arXiv:1802.02445].

[259] T. Liu, K. Melnikov and A. A. Penin (2019), [arXiv:1906.10899].

[260] H. A. Chawdhry et al., JHEP 02, 057 (2020), [arXiv:1911.00479].

[261] S. Kallweit, V. Sotnikov and M. Wiesemann, Phys. Lett. B 812, 136013 (2021), [arXiv:2010.04681].

[262] H. A. Chawdhry et al., JHEP 09, 093 (2021), [arXiv:2105.06940].

[263] S. Badger et al. (2023), [arXiv:2304.06682].

[264] M. Czakon, A. Mitov and R. Poncelet, Phys. Rev. Lett. 127, 15, 152001 (2021), [Erratum: Phys.Rev.Lett. 129, 119901 (2022)], [arXiv:2106.05331].

[265] H. B. Hartanto et al., Phys. Rev. D 106, 7, 074016 (2022), [arXiv:2205.01687].

[266] M. Alvarez et al., JHEP 03, 129 (2023), [arXiv:2301.01086].

[267] G. P. Salam and E. Slade, JHEP 11, 220 (2021), [arXiv:2106.08329].

[268] M. A. Ebert et al., JHEP 04, 102 (2021), [arXiv:2006.11382].

[269] S. Amoroso et al., in "11th Les Houches Workshop on Physics at TeV Colliders: PhysTeV Les Houches," (2020), [arXiv:2003.01700].

[270] S. Catani et al., Phys. Rev. Lett. 130, 11, 111902 (2023), [arXiv:2210.07846].

[271] S. Camarda et al., Eur. Phys. J. C 80, 3, 251 (2020), [Erratum: Eur.Phys.J.C 80, 440 (2020)], [arXiv:1910.07049].

[272] R. Gavin et al., Comput. Phys. Commun. 182, 2388 (2011), [arXiv:1011.3540].

[273] S. Camarda et al., Eur. Phys. J. C 82, 5, 492 (2022), [arXiv:2202.10343].

[274] R. V. Harlander, S. Liebler and H. Mantler, Comput. Phys. Commun. 184, 1605 (2013), [arXiv:1212.3249].

[275] T. Carli et al., Eur. Phys. J. C 66, 503 (2010), [arXiv:0911.2985].

[276] T. Kluge, K. Rabbertz and M. Wobisch, in "14th International Workshop on Deep Inelastic Scattering (DIS 2006)," 483, Tsukuba, Japan, April 20-24 (2006), [hep-ph/0609285].

[277] D. Britzger et al., in "Proceedings, XX. International Workshop on Deep-Inelastic Scattering and Related Subjects (DIS 2012)," 217, Bonn, Germany, March 26-30 (2012), [arXiv:1208.3641].

[278] S. Carrazza et al., JHEP 12, 108 (2020), [arXiv:2008.12789].

[279] M. Czakon, D. Heymes and A. Mitov (2017), [arXiv:1704.08551].

[280] D. Britzger et al., Eur. Phys. J. C 79, 10, 845 (2019), [Erratum: Eur.Phys.J.C 81, 957 (2021)], [arXiv:1906.05303].

[281] D. Maître, J. Phys. Conf. Ser. 1525, 1, 012014 (2020).

[282] D. Maître, G. Heinrich and M. Johnson, PoS LL2016, 016 (2016), [arXiv:1607.06259].

[283] Y. L. Dokshitzer, D. Diakonov and S. I. Troian, Phys. Rept. 58, 269 (1980).

[284] G. Parisi and R. Petronzio, Nucl. Phys. B154, 427 (1979).

[285] G. Curci, M. Greco and Y. Srivastava, Nucl. Phys. B159, 451 (1979).

[286] A. Bassetto, M. Ciafaloni and G. Marchesini, Nucl. Phys. B163, 477 (1980).

[287] J. C. Collins and D. E. Soper, Nucl. Phys. B193, 381 (1981), [Erratum: Nucl. Phys.B213,545(1983)].

[288] J. C. Collins and D. E. Soper, Nucl. Phys. B197, 446 (1982).

[289] J. Kodaira and L. Trentadue, Phys. Lett. 112B, 66 (1982).

[290] J. Kodaira and L. Trentadue, Phys. Lett. 123B, 335 (1983).

[291] J. C. Collins, D. E. Soper and G. F. Sterman, Nucl. Phys. B250, 199 (1985).

[292] S. Catani et al., Nucl. Phys. B407, 3 (1993).

[293] C. W. Bauer et al., Phys. Rev. D63, 114020 (2001), [hep-ph/0011336].

[294] C. W. Bauer, D. Pirjol and I. W. Stewart, Phys. Rev. D65, 054022 (2002), [hep-ph/0109045].

[295] T. Becher, A. Broggio and A. Ferroglia, Lect. Notes Phys. 896, pp.1 (2015), [arXiv:1410.1892].

[296] S. Catani et al., Phys. Lett. B269, 432 (1991).

[297] N. Brown and W. J. Stirling, Phys. Lett. B252, 657 (1990).

[298] W. Bartel et al. (JADE), Z. Phys. C33, 23 (1986), [e53(1986)].

[299] N. Kidonakis, G. Oderda and G. F. Sterman, Nucl. Phys. B531, 365 (1998), [hep-ph/9803241].

[300] R. Bonciani et al., Phys. Lett. B575, 268 (2003), [hep-ph/0307035].

[301] A. Banfi, G. P. Salam and G. Zanderighi, JHEP 03, 073 (2005), [hep-ph/0407286].

[302] D. de Florian and M. Grazzini, Phys. Rev. Lett. 85, 4678 (2000), [hep-ph/0008152].

[303] G. Bozzi et al., Nucl. Phys. B737, 73 (2006), [hep-ph/0508068]; http://theory.fi.infn.it/grazzini/codes.html.

[304] G. Bozzi et al., Phys. Lett. B696, 207 (2011), [arXiv:1007.2351].

[305] T. Becher and M. Neubert, Eur. Phys. J. C71, 1665 (2011), [arXiv:1007.4005].

[306] http://cute.hepforge.org/.

[307] D. de Florian et al., JHEP 06, 132 (2012), [arXiv:1203.6321]; http://theory.fi.infn.it/grazzini/codes.html.

[308] C. Balazs and C. P. Yuan, Phys. Rev. D56, 5558 (1997), [hep-ph/9704258].

[309] S. Catani et al., JHEP 12, 047 (2015), [arXiv:1507.06937].

[310] A. Banfi et al., Phys. Lett. B715, 152 (2012), [arXiv:1205.4760].

[311] M. Grazzini et al., JHEP 08, 154 (2015), [arXiv:1507.02565].

[312] J. M. Campbell et al., JHEP 03, 080 (2023), [arXiv:2210.10724].

[313] D. de Florian and M. Grazzini, Nucl. Phys. B704, 387 (2005), [hep-ph/0407241].

[314] T. Becher and G. Bell, JHEP 11, 126 (2012), [arXiv:1210.0580].

[315] A. Banfi et al., Phys. Rev. Lett. 109, 202001 (2012), [arXiv:1206.4998]; T. Becher, M. Neubert and L. Rothen, JHEP 10, 125 (2013), [arXiv:1307.0025].

[316] I. W. Stewart et al., Phys. Rev. D 89, 5, 054001 (2014), [arXiv:1307.1808].

[317] J. M. Campbell et al., JHEP 04, 106 (2023), [arXiv:2301.11768].

[318] I. W. Stewart, F. J. Tackmann and W. J. Waalewijn, Phys. Rev. Lett. 106, 032001 (2011), [arXiv:1005.4060].

[319] Y.-T. Chien et al., Phys. Rev. D87, 1, 014010 (2013), [arXiv:1208.0010]; T. T. Jouttenus et al., Phys. Rev. D88, 5, 054031 (2013), [arXiv:1302.0846].

[320] M. Dasgupta et al., JHEP 10, 126 (2012), [arXiv:1207.1640].

[321] V. Ahrens et al., JHEP 09, 097 (2010), [arXiv:1003.5827].