Here, QN+1 is the hardness scale related to the emission of the (N + 1)th particle. As advertised before, this scale defines two different regimes, namely the parton shower region QN+1 < Qcut and the matrix-element region QN+1 > Qcut. While in the parton shower only expression, Eq. (5.304), such a differentiation is not being made, in contrast here the description of the emission of the (N + 1)th particle is susceptible to this difference.

 Taking a closer look at the square bracket in the first line it becomes apparent that it does not integrate to unity any more, due to the phase-space constraint — encoded through Θ(Qcut − QN+1) — present in the second term, the emission term. The missing hard emissions are of course supplied through the second line, explicit through the complementary constraint Θ(QN+1 − Qcut). There are two slight mismatches which prevent this term to also account for the pieces lacking in the first line to integrate to unity. First of all, there is a mismatch in the form of the emission term: the exact matrix element for the (N + 1)-particle state is different from the exact matrix element for the N-particle state convoluted with the parton shower kernel. This does not come as a surprise, as this was the very reason multijet merging has been introduced in the first place.

On top of this, the phase space available for this additional emission, the (N + 1)th particle, differs in both lines. In the first line the upper limit on the available phase space is given by the relevant parton shower starting or resummation scale defined through the N particle kinematics, µN , while in the second line it is the potentially different (N + 1)-particle kinematics, which defines the corresponding scale µN+1. While in processes with a fixed hardest scale, such as electron–positron annihilations to jets it is safe to assume that typically µN and µN+1 are very similar or even identical, this is not true for processes at hadron colliders such as the production of lepton pairs (the Drell-Yan process) in association with jets.

There, the emission of additional jets typically may offer phase-space regions associated with larger scales than processes without these additional jets — in the example of Drell–Yan-type processes this would correspond to jet emissions taking place at transverse momentum scales above the invariant mass of the lepton pair. Taken together, this leads to the fact that the combined contributions from the hard and soft first emissions will not exactly combine with the no-emission term to yield unity. This has been, in a somewhat sloppy use of language been coined "unitarity violation". It consequently leads to a variation of the total cross-section related to the inclusive sample produced with respect to the Born result. Postponing details to later parts of this section, it should be noted that this effect, while present in some of the more widely used multijet merging implementations, has been taken care of by the UMEPS algorithm [733, 806], which conserves inclusive cross-sections by suitably reshuffling different contributions.

Taking into account first emissions only, but combining many Born matrix elements up to a maximal number Nmax of external legs promotes the simple expression of Eq. (5.349) to

$$\begin{split} \mathsf{d}\sigma &= \sum\_{n=N}^{N\_{\max}-1} \left\{ \mathsf{d}\Phi\_{n}\mathcal{B}\_{n} \overbrace{\left[\prod\_{j=N}^{n} \Theta(Q\_{j}-Q\_{\mathrm{cut}})\right]}^{n} \overbrace{\left[\prod\_{j=N}^{n-1} \Delta\_{j}^{(\mathsf{K})}(t\_{j},t\_{j+1})\right]}^{n-1 \text{ emissions of internal lines}} \right\} \\ & \times \left[ \underbrace{\Delta\_{N}^{(\mathsf{K})}(t\_{N},t\_{\mathsf{c}})}\_{\mathsf{no emission}} + \underbrace{\int\_{t\_{\mathsf{c}}}^{t\_{\mathsf{N}}} \mathrm{d}\Phi\_{1}\,\mathsf{K}\_{N}\Delta\_{N}^{(\mathsf{K})}(t\_{N},t\_{N+1})\Theta(Q\_{\mathrm{cut}}-Q\_{N+1})}\_{\mathsf{next emission} \text{ no jet} \&\ \mathsf{k} \text{ HE emission}} \right] \right\} \\ &+ \quad \mathrm{d}\Phi\_{N\_{\max}}\mathcal{B}\_{N\_{\max}} \left[ \prod\_{j=N}^{N\_{\max}} \Theta(Q\_{j+1}-Q\_{\mathrm{cut}}) \right] \left[ \prod\_{j=N}^{N\_{\max}-1} \Delta\_{j}^{(\mathsf{K})}(t\_{j},t\_{j+1}) \right] \\ &= \sum\_{t\_{N\_{\max}}} \end{split}$$

$$\begin{split} \times \left[ \Delta\_{N\_{\text{max}}}^{\langle \mathcal{K} \rangle} (t\_{N\_{\text{max}}}, t\_{\text{c}}) + \int\_{t\_{\text{c}}}^{t\_{N\_{\text{max}}}} \mathrm{d}\Phi\_{1} \,\mathcal{K}\_{N\_{\text{max}}} \Delta\_{N\_{\text{max}}}^{\langle \mathcal{K} \rangle} (t\_{N\_{\text{max}}}, t\_{N\_{\text{max}}+1}) \right] \\ \cdot \Theta(Q\_{N\_{\text{max}}} - Q\_{N\_{\text{max}}+1}) \Big]. \end{split} \tag{5.350}$$

Note that in the contribution from the Nmax configuration, the phase space filled by the parton shower is not constrained by the jet cut Qcut but by the jet measure of the last emission filled by the matrix element. This allows the parton shower to account — within its limitations — for even higher jet multiplicities.

At this point it should of course be stressed that further emissions are based on splitting kernels supplemented with the phase-space veto:

$$\mathcal{K}\_N(\Phi\_1) \stackrel{\text{MEPS}}{\longrightarrow} \mathcal{K}\_N^{$$

Taking this into account and applying it to the subsequent emissions encoded in the parton shower evolution operator, cf. Eq. (5.310), E (K) N becomes

$$\begin{split} \mathcal{E}\_{N}^{(\mathcal{K},$$

Inserting this into the merging equation Eq. (5.350) above results in

Multijet merging of parton showers and matrix elements 385

$$\begin{split} \mathrm{d}\sigma &= \sum\_{n=N}^{N\_{\mathrm{max}}-1} \mathrm{d}\Phi\_{n} \mathcal{B}\_{n} \, \Theta(Q\_{n} - Q\_{\mathrm{cut}}) \, \mathcal{E}\_{n}^{(\mathcal{K}, < Q\_{\mathrm{cut}})} (\mu\_{n}^{2}, t\_{c}) \\ &+ \mathrm{d}\Phi\_{N\_{\mathrm{max}}} \mathcal{B}\_{N\_{\mathrm{max}}} \, \Theta(Q\_{N\_{\mathrm{max}}} - Q\_{\mathrm{cut}}) \, \mathcal{E}\_{N\_{\mathrm{max}}}^{(\mathcal{K}, < Q\_{N\_{\mathrm{max}}})} (\mu\_{N\_{\mathrm{max}}}^{2}, t\_{c}) \end{split} \tag{5.353}$$

This improved description has been implemented for a variety of parton showers, most of which use an ordering parameter related to the transverse momentum of the emissions and thus do not exhibit any real mismatch of evolution and hardness parameters. For a more in-depth discussion of the various realizations the reader is referred to the literature [632, 732]. Multijet merging with an angular-ordered parton shower in the framework of HERWIG++ has been explored in [612], where the impact of not employing truncated showering in the merging has been analysed, too.

#### <IP_ADDRESS> Alternative: The MLM method

The MLM prescription for multijet merging [147, 744] was introduced and further refined in parallel to the method discussed so far. It is based on slightly different ideas: most notably it aims at using traditional parton shower routines like the ones implemented in PYTHIA [853] (virtuality or p⊥-ordered shower) or in HERWIG [414] (angular-ordered shower) which are accessed through a standardized interface between general matrix elements and parton showers [145] without any direct alteration in the code. The structure of the resulting algorithm implies that vetoes on unwanted emissions can only be applied a posteriori, i.e. after the parton shower evolution has finished. In practice this is realized by reclustering the partons after the parton shower, i.e. at the hadronization scale tc, into jets and by comparing these jets with the original ones, from the parton level at scale Qcut. If an additional jet with respect to the original ones has been produced, or, conversely, if a jet has been lost, the event then is rejected. This algorithm of course incorporates a slightly different way of treating emissions off intermediate legs, which here are re-interpreted as radiation off external legs, by a suitable definition of parton shower starting conditions.12

In its original version, the MLM prescription uses a simple cone definition as the jet criterion to generate the parton configurations at the matrix-element level. The scales of strong couplings are reconstructed by using a backward clustering with a k⊥ measure. Then the accepted configurations are passed on to the parton-shower routines. They in turn typically reconstruct the parton-shower starting scales of a multi-parton configuration by directly inspecting colour connections, without any backward clustering and therefore partially neglecting intermediate legs and the radiation originating from them. Having the starting conditions at hand, the parton shower is invoked without any constraint. After it has terminated at the hadronization scale tc, the original partons stemming from the matrix element are exactly matched to the jets present at parton-shower level, again defined by the cone algorithm with parameters, which in principle may slightly differ from the ones applied at matrix-element level. If such a

12In such a treatment, the full capture of truncated showering effects is not guaranteed, and as a consequence some residual sub-leading logarithmic terms may be left out [632], which would be present in the original parton shower algorithm.

#### 386 QCD to All Orders

one-to-one match is not possible, either due to extra, unwanted jets being produced in the parton shower or due to "losing" jets in the parton shower, the event is rejected. This means that the Sudakov rejection factors that are applied locally, either through analytic Sudakov form factors or by using the shower, as explained above, are applied inclusively over the full parton configuration. As a by-product of this treatment, the MLM approach penalizes "losing" jets which is not the case in the original merging prescriptions. This algorithm has originally been implemented in ALPGEN [743, 744]; a variant of it using the Durham k⊥-algorithm has later been provided in the MADGRAPH framework [146].

Despite the subtle differences in the different approaches, by far and large a good agreement of predictions obtained between both methods can be observed. Respective results have been reported for the case of W+jets production at the TEVATRON and the LHC for example in [147].

#### <IP_ADDRESS> Extensions: Dealing with photons

The extension of the merging algorithm to also include photons is fairly straightforward: treating the emission of photons "democratically", i.e. on the same footing as QCD emissions allows to embed them into a multijet merging such that both the number of QCD particles and photons will vary. Evaluating matrix elements at the tree-level with NQCD QCD particles and Nγ photons is not a problem and can be dealt with standard technology. There is also no difficulty in supplementing the jet definition encoded in QJ and a corresponding cut Qcut with an isolation criterion Qγ and a corresponding cut Qiso for the emission of photons. From the parton shower side, things are similarly trivial and basically amount to supplementing q → qγ splitting kernels, which typically can be obtained from the q → qg ones by suitably replacing coupling factors. In principle it is also possible — and probably even desirable — to also include γ → f ¯f splittings and to have different infrared cut-offs for the QCD and the QED part of the parton shower, in other words to supplement tc with a tc,(QED). The latter typically would be chosen on scales of the order of the π 0 mass or similar, which is of course possible due to the absence of a Landau pole in QED in the soft regime.

Of course, the same logic could also be applied to extend this treatment to, e.g., the emission of the weak vector bosons, the W± and Z bosons. There is one caveat, though, related to the fact that the coupling of, say, W bosons to fermions is chiral and therefore highly sensitive to the spin of the fermions and, consequently, introduces non-trivial spin correlations in the parton shower. Including them in a systematic way would imply that the easily implemented structure of more or less independent emissions in a simple probabilistic fashion would have to be augmented with spincorrelation matrices of the kind discussed in [825] or a treatment similar to the one discussed in [779]. However, for the emission of one boson only some matrix-element corrections could be applied which would by far and large capture such effects. This has been studied in more detail in the framework of the PYTHIA event generator [396].

#### <IP_ADDRESS> Aside: unitarization of multijet merging: UMEPS

The idea underlying the "unitarization" of multijet merging is to conserve the fixedorder cross-section of the inclusive process with lowest multiplicity — in the MEPS master equation, Eq. (5.353), this would be the cross-section of the N-particle process,

$$\int \mathrm{d}\sigma^{(LO)} \, \, = \int \mathrm{d}\Phi\_N \mathcal{B}\_N \, \stackrel{!}{=} \int \mathrm{d}\sigma^{(\text{UMEPS})} \,. \tag{5.354}$$

This is in principle achieved by shuffling contributions from higher multiplicities to lower ones, thus compensating for the mismatch of fixed-order real emission crosssections and their parton-shower approximation above Qcut. Such a mismatch manifested itself in Eq. (5.349) as a difference between BN+1 and BN ⊗ KN .

To see in more detail how this works, consider first the expression for the Sudakov form factor, Eq. (5.251), which can be decomposed as

$$\begin{split} \Delta\_{N}^{(\mathcal{K})}(T,t) &= \Delta\_{N}^{(\mathcal{K})}(T,t;Q\_{N+1}\geq Q\_{\text{cut}}) \times \Delta\_{N}^{(\mathcal{K})}(T,t;Q\_{N+1}$$

Using this in Eq. (5.349) allows to rewrite the first line as

$$\begin{split} \mathrm{d}\Phi\_{N}\mathcal{B}\_{N}\left[\Delta\_{N}^{(\mathcal{K})}(\mu\_{N}^{2},t\_{c}) + \int\_{t\_{c}}^{\mu\_{N}^{2}} \mathrm{d}\Phi\_{1}\mathcal{K}\_{N}\Delta\_{N}^{(\mathcal{K})}(\mu\_{N}^{2},t\_{N+1})\Theta(Q\_{\mathrm{cut}} - Q\_{N+1})\right] \\ = \mathrm{d}\Phi\_{N}\mathcal{B}\_{N}\Delta\_{N}^{\geq Q\_{\mathrm{cut}}(\mathcal{K})}(\mu\_{N}^{2},t\_{c}) \\ \quad \times \left[\Delta\_{N}^{$$

Again, because of the identical kernels in the individual terms and the identical phasespace constraints, the square bracket integrates to unity. Using the probabilistic nature of the parton shower, the Sudakov form factor, being interpreted as the no-emission probability, can be recast as unity minus the emission probability. In other words, in

$$
\Delta\_N^{\geq Q\_{\rm cut}, \langle \mathcal{K} \rangle} (\mu\_N^2, t\_c) = 1 - \int\_{t\_c}^{\mu\_N^2} \mathrm{d}\Phi\_1 \, \mathsf{K}^{\geq Q\_{\rm cut}} \, \Delta\_N^{\geq Q\_{\rm cut}, \langle \mathcal{K} \rangle} (\mu\_N^2, t) \tag{5.357}
$$

the integrated parton shower emission rate above Qcut is given by the second term.

#### 388 QCD to All Orders

The unitarization of the cross-section is now achieved by suitably replacing the kernel in the integrand of Eq. (5.357),

$$\Delta\_N^{\geq Q\_{\mathrm{cut}}, (\mathbb{K})} (\mu\_N^2, t\_c) = 1 - \int\_{t\_c}^{\mu\_N^2} \mathrm{d}\Phi\_1 \, \frac{\mathcal{B}\_{N+1}}{\mathcal{B}\_N} \, \Theta(Q\_{N+1} - Q\_{\mathrm{cut}}) \, \Delta\_N^{\geq Q\_{\mathrm{cut}}, (\mathbb{K})} (\mu\_N^2, t) \tag{5.358}$$

Merging now the first emission through the corresponding Born matrix element yields

$$\begin{split} \mathrm{d}\sigma &= \, \mathrm{d}\Phi\_{N}\mathcal{B}\_{N}\, \bar{\Delta}\_{N}^{\geq Q\_{\mathrm{cut}}(\mathcal{K})}(\mu\_{N}^{2},t\_{\mathrm{c}}) \\ & \times \left[\Delta\_{N}^{$$

Similar replacements are applied in Eq. (5.350) to all brackets encoding the parton shower evolution, with the exception of the Nmax-term. Consequently, all N-parton inclusive cross-sections are governed by the input cross-sections associated with the BN . Some comments are in order here. First of all, it is possible that a parton shower cannot produce all particles or phase-space configurations entering these cross-sections. In such a case the corresponding fixed-order cross-sections must be regarded as genuine corrections to the parton shower and therefore they will be just added. This is also true for the impact of the phase-space constraints given by the Θ functions. It is not always guaranteed that integrating over the one-particle phase space will yield a lower-multiplicity state that passes the Qcut criterion — sometimes one loses more than one jet. Handling such contributions is ambiguous, since they could be regarded as genuine fixed-order corrections as discussed in [806] or as contributions to even lower multiplicity states as in [733].

#### 5.5.2 Hybrids of merging and matching: MENLOPS

In a first step towards a full multijet merging on the basis of NLO multijet matrix elements, and with NLO matching and LO merging methods well established it is worth discussing a hybrid between these two methods. By now this has been established as the MENLOPS method, and it has been introduced and implemented for various processes in [608, 629]. It combines matching, either with the POWHEG or the MC@NLO method, of QCD NLO matrix elements for the lowest-multiplicity final states with an LO multijet merging of higher-multiplicity matrix elements. This implies that the NLO matching has to be performed in such a way that the real-emission part of the NLO correction does not produce any additional jets; fully in line with the LO merging paradigm therefore the emission phase space has to be constrained through a jet measure QJ and a corresponding cut Qcut. This requirement is realized by suitably multiplying with a Sudakov form factor. In addition, the no-jet emission constraint must act on the phase space available for both the real-emission correction encoded through an exact matrix element and on the parton shower parts.

While this is, by far and large, not too hard to implement, there is another somewhat more nagging problem. The jet-exclusive cross-section for the lowest-multiplicity final state will be modified by the local K-factor discussed in Sections 5.4.3 and 5.4.4, while the higher multiplicities will not be corrected and be at Born level only. This can lead to discontinuities in the radiation pattern, in particular of the first, hardest emission, and especially in those cases, where the K factors are fairly different from unity, like, for instance, in the case of Higgs production through gluon fusion. In order to remedy this situation, the higher-multiplicity part could be multiplied by an interpolating K-factor, kN , capturing the NLO correction to the lowest multiplicity N. For example, for the MENLOPS method, where the NLO matching part is realized through the MC@NLO algorithm, this interpolating K-factor may be given by something like

$$k\_N(\Phi\_{N+1}) = \frac{\tilde{\mathcal{B}}\_N}{\mathcal{B}\_N} \left(1 - \frac{\mathcal{H}\_N}{\mathcal{B}\_{N+1}}\right) + \frac{\mathcal{H}\_N}{\mathcal{B}\_{N+1}} \longrightarrow \begin{cases} \tilde{\mathcal{B}}\_N/\mathcal{B}\_N \text{ for soft emissions} \\ 1 \qquad \text{for hard emissions.} \end{cases} \tag{5.360}$$

Such an interpolation works fairly well, since in the soft limit the real correction term and its subtraction become very similar and therefore the hard remainder is given by

$$\mathcal{H}\_N = \mathcal{R}\_N - \mathcal{S}\_N \stackrel{\text{soft}}{\longrightarrow} 0,\tag{5.361}$$

while in the hard limit the subtraction term is not very prominent and therefore

$$\mathcal{H}\_N = \mathcal{R}\_N - \mathcal{S}\_N \stackrel{\text{hard}}{\longrightarrow} \mathcal{R}\_N = \mathcal{B}\_{N+1} \,. \tag{5.362}$$

Concentrating again on contributions up to the first emissions only, the differential cross-section for such a sample is given by

dσ = dΦN Θ(QN − Qcut) B˜N × " ∆ (K) N (µ 2 N , tcut) + µ 2 ZN tcut dΦ1 KN ∆ (K) N (µ 2 N , tN+1)Θ(Qcut − QN+1) # + dΦN+1 Θ(QN − Qcut)Θ(Qcut − QN+1) HN ∆ (K) N (µ 2 N+1, tN+1) + dΦN+1 kN (ΦN+1) Θ(QN+1 − Qcut) BN+1∆ (K) N (µ 2 N+1, tN+1) × " ∆ (K) N+1(tN+1, tcut) + tZN+1 tcut dΦ1 KN+1∆ (K) N+1(tN+1, tN+2)Θ(Qcut − QN+2) # + dΦN+2 kN (ΦN+2) Θ(QN+2 − Qcut) BN+2∆ (K) N+1(µ 2 N+2, tN+2)∆(K) N (tN+2, tN+1) × " ∆ (K) N+1(tN+2, tcut) + . . . . (5.363)

#### 390 QCD to All Orders

Here, the first three lines correspond to an MC@NLO simulation, where the phase space for the (N + 1)th particle is constrained such that it does not yield another jet — this is what the Θ–functions Θ(Qcut −QN+1) are there for. In order to make the jet counting more explicit, another Θ-function, Θ(QN −Qcut) has been added to highlight that all other QCD particles in the N-particle Born-level final state should be jets. Similarly, the fourth and fifth lines encapsulate the first additional tree-level matrix element, merged into the sample, and supplemented with the interpolating K factor kN . Again, no further jet emission is allowed. This could now continue like indicated by the sixth line, to include higher and higher jet multiplicities. In any case, any further emissions, as discussed already for the multijet merging at leading order, cannot result in additional unwanted jets.

#### 5.5.3 Multijet merging at next-to-leading order: MEPS@NLO

#### <IP_ADDRESS> Basic idea

The idea underlying multijet merging with next-to-leading order matrix elements, MEPS@NLO, is identical to multijet merging at leading order: towers of matrix elements with increasing jet multiplicities are combined into one inclusive sample in such a way that no double-counting occurs. Of course, as before, it is important to maintain the accuracy of both the matrix elements, i.e. the cross-sections related to the processes as well as the fixed-order accuracy of the first emission, and the parton shower, i.e. the resummation of leading and the next-to-leading logarithms encoded in the showering.

This has been first achieved in [558, 631], where the first missing terms have explicitly been shown to be of order α 2 sL 3/N2 c — the colour-suppressed sub-leading logarithms beyond shower accuracy. Alternative methods and implementations have been presented in [535], and in [737, 806]. The latter two are closely related to one another, and both guarantee the proper "unitarization" of the emissions in line with the leading-order treatment by the same authors [733, 806].

#### <IP_ADDRESS> First emission(s), once more

Following the original presentation of a multijet-merging method of next-to leadingorder matrix elements in [558, 631], this MEPS@NLO method can be understood as a merging of individual MC@NLO simulations. Some additional terms need to be included, though, to guarantee the fixed-order and logarithmic correctness. Concentrating first on the merging of matrix elements for the production of an N and an (N + 1)-particle state only, and taking into account the first emissions only, a naive addition of two MC@NLO simulations for both would yield the (wrong) cross-section

$$\begin{split}d\sigma^{\text{(wong)}} &= \quad \mathrm{d}\Phi\_{N}\Theta(Q\_{N}-Q\_{\text{cut}})\tilde{B}\_{N} \\ & \quad \times \left[\Delta\_{N}^{(\mathcal{K})}(\mu\_{N}^{2},t\_{\text{c}}) + \int\_{t\_{\text{c}}}^{\mu\_{N}^{(\mathcal{K})}} \mathrm{d}\Phi\_{N}\mathcal{K}\_{N}^{(\mathcal{K})}(\mu\_{N}^{2},t\_{N+1})\Theta(Q\_{\text{cut}}-Q\_{N+1})\right] \\ & \quad + \,\mathrm{d}\Phi\_{N+1}\Theta(Q\_{N}-Q\_{\text{cut}})\Theta(Q\_{\text{cut}}-Q\_{N+1})\,\mathcal{H}\_{N}\Delta\_{N}^{(\mathcal{K})}(\mu\_{N}^{2},t\_{N+1}) \\ & \quad + \,\mathrm{d}\Phi\_{N+1}\Theta(Q\_{N+1}-Q\_{\text{cut}})\,\tilde{B}\_{N+1} \\ & \quad \times \left[\Delta\_{N+1}^{(\mathcal{K})}(t\_{N+1},t\_{\text{c}}) + \int\_{t\_{\text{c}}}^{t\_{N+1}} \mathrm{d}\Phi\_{1}\,\mathcal{K}\_{N+1}\Delta\_{N+1}^{(\mathcal{K})}(t\_{N+1},t\_{N+2})\right] \\ & \quad + \,\mathrm{d}\Phi\_{N+2}\Theta(Q\_{N+1}-Q\_{\text{cut}})\,\mathcal{H}\_{N+1}\Delta\_{N}^{(\mathcal{K})}(\mu\_{N}^{2},t\_{N+1})\Delta\_{N+1}^{(\mathcal{K})}(t\_{N+1},t\_{N+2}). \end{split} \tag{5.344}$$

The first three lines would correspond to an MC@NLO simulation, where, in complete agreement with Eq. (5.363), the phase space for the (N + 1)th particle is constrained such that it does not yield another jet, realized by Θ(Qcut−QN+1) in the two emission terms, soft and hard. The next three lines would then stand for the next MC@NLO simulation, for a process with one more particle in the final state.

Naively, it seems as if there is no double counting present here. However, this is not entirely true. To see this, consider the emission term in the first, soft part of the lowest order MC@NLO simulation, the term KN ∆ (K) N (µ 2 N , tN+1) in the first square bracket and the corresponding hard emission part in the line below. Closer inspection reveals that at first order in αs, there is some unwanted contribution to the phase space of the next MC@NLO simulation. It stems from the expansion of the Sudakov form factor accounting for no emissions harder than tN+1 in the combined soft and hard radiation pattern, ranging over the full emission phase space from tN+1 up to µ 2 N :

$$\begin{aligned} & \left[ \mathrm{d}\Phi\_N \int \mathrm{d}\Phi\_1 \, \tilde{\mathcal{B}}\_N \mathcal{K}\_N + \mathrm{d}\Phi\_{N+1} \mathcal{H}\_N \right] \Theta(Q\_{\mathrm{cut}} - Q\_{N+1}) \, \Delta\_N^{(\mathbb{K})} (\mu\_N^2, t\_{N+1}) \\ &= \left[ \mathrm{d}\Phi\_N \int \mathrm{d}\Phi\_1 \, \left( \mathcal{B}\_N \otimes \mathcal{K}\_N + \mathcal{H}\_N \right) + \mathcal{O} \left( \alpha\_s^2 \right) \right] \Theta(Q\_{\mathrm{cut}} - Q\_{N+1}) \, \Delta\_N^{(\mathbb{K})} (\mu\_N^2, t\_{N+1}) \\ &= \left. \mathrm{d}\Phi\_{N+1} \mathcal{B}\_{N+1} \, \Theta(Q\_{\mathrm{cut}} - Q\_{N+1}) \left[ 1 - \int\_{t\_{N+1}}^{\mu\_N^2} \mathrm{d}\Phi\_1 \mathcal{K}\_N + \mathcal{O} \left( \alpha\_s^2 \right) \right] \right. \end{aligned} \tag{5.365}$$

The second part of the bracket in the last line therefore interferes with the emissions of the MC@NLO simulation of the incremented multiplicity. This double-counting of emissions must obviously be avoided. There are various way of achieving this, the simplest one is by adding the second term in the square bracket above to the (N + 1)-

#### 392 QCD to All Orders

MC@NLO simulation:

$$\begin{split} \mathsf{d}\sigma &= \mathsf{d}\Phi\_{N}\Theta(Q\_{N} - Q\_{\mathrm{cut}})\vec{\mathcal{B}}\_{N} \\ &\times \left[\Delta\_{N}^{(\mathsf{K})}(\mu\_{N}^{2}, t\_{\mathrm{c}}) + \int\_{t\_{\mathrm{c}}}^{\mu\_{N}^{2}} \mathsf{d}\Phi\_{1}\mathcal{K}\_{N}\Delta\_{N}^{(\mathsf{K})}(\mu\_{N}^{2}, t\_{N+1})\Theta(Q\_{\mathrm{cut}} - Q\_{N+1})\right] \\ &+ \mathsf{d}\Phi\_{N+1}\Theta(Q\_{N} - Q\_{\mathrm{cut}})\Theta(Q\_{\mathrm{cut}} - Q\_{N+1})\ \mathcal{H}\_{N}\Delta\_{N}^{(\mathsf{K})}(\mu\_{N}^{2}, t\_{N+1}) \end{split}$$

$$\begin{split} &+ \operatorname{d}\Phi\_{N+1}\Theta(Q\_{N+1} - Q\_{\text{cut}})\tilde{\mathcal{B}}\_{N+1}\Delta\_{N}^{(\mathcal{K})}(\mu\_{N+1}^{2}, t\_{N+1}) \Big(1 + \frac{\mathcal{B}\_{N+1}}{\tilde{\mathcal{B}}\_{N+1}} \int\_{t\_{N+1}}^{\mu\_{N}^{2}} \operatorname{d}\Phi\_{1}\mathcal{K}\_{N} \Big) \\ & \qquad \times \Big[\Delta\_{N+1}^{(\mathcal{K})}(t\_{N+1}, t\_{\mathcal{c}}) + \int\_{t\_{\mathcal{c}}}^{t\_{N+1}} \operatorname{d}\Phi\_{1}\mathcal{K}\_{N+1}\Delta\_{N+1}^{(\mathcal{K})}(t\_{N+1}, t\_{N+2}) \Big] \\ &+ \operatorname{d}\Phi\_{N+2}\Theta(Q\_{N+1} - Q\_{\text{cut}})\Theta(Q\_{\text{cut}} - Q\_{N+2})\mathcal{H}\_{N+1} \\ & \qquad \times \Delta\_{N}^{(\mathcal{K})}(\mu\_{N+1}^{2}, t\_{N+1})\Delta\_{N+1}^{(\mathcal{K})}(t\_{N+1}, t\_{N+2}) . \end{split}$$

(5.366)

As in the MENLOPS method described in Section 5.5.2, the first three lines merely describe an MC@NLO simulation of the lowest multiplicity sample with an additional jet veto applied to all emissions. In a similar way, the next lines describe an MC@NLO simulation for the next higher multiplicity, modified by the term in round brackets in the fourth line compensating for a double-counting of terms stemming from the lowest multiplicity MC@NLO.

As already worked out, this term compensates contributions stemming from the Sudakov form factor encoding the veto of hard emissions from the lower multiplicity, the two first lines. Such a compensation must be introduced into the formalism in order to maintain the NLO accuracy of the overall procedure. The fact that this double-counting of NLO terms is introduced by the parton shower (or analytic Sudakov form factors encoding the jet veto) renders these contributions potentially hard to understand at first. Realizing, however, that Sudakov form factors encode some of the NLO corrections in a leading-logarithmic approximation, this should not come as a surprise. In fact, the treatment here actually is fully analogous to the one of the Sudakov form factor in its interaction with higher-order matrix elements, presented later in Section 5.6.1, and made explicit in Eq. (5.372) there.

Naively, such terms also seem to be hard to implement when generating the Sudakov rejection directly from the parton shower rather than analytically. This turns out not to be entirely true. In close analogy to the discussion of vetoed emissions, cf. Eq. (5.348), it can be seen that these terms correspond to a vetoed emission off the N-particle state only. Instead of vetoing an event when an unwanted emission takes place, it is merely the hard emission itself that should be vetoed.

Of course, the very same logic could be iterated to include even higher multiplicities. In addition, further leading-order matrix elements could be added following the same reasoning with an interpolating K-factor as in the MENLOPS algorithm.

#### <IP_ADDRESS> Scale setting in MEPS@NLO

One remaining problem is the definition of scales, and in particular the nodal values for the individual emissions and, correspondingly, the renormalization scales of the strong coupling in the various terms. Simply put, the same logic as in the LO case, Section 5.5.1, will apply. Starting from Born-type configurations, the partons are clustered backwards, yielding, as before, a parton shower history for the configuration. Its nodal values, ti will be used in the Sudakov rejection factors, and, together with the rest of the kinematics of the ith splitting the corresponding renormalization scale µ (i) R can be deduced. As in Eq. (5.345), then the overall scale for a process with M powers of the strong coupling in the 2 → 2 core process and with m additional QCD emissions will be given as

$$\alpha\_{\rm s}^{M+m}(\mu\_R^2) = \alpha\_{\rm s}^n(\mu\_R^2) = \alpha\_{\rm s}^M(\mu\_{R\text{(core)}}^2) \cdot \prod\_{i \in m} \alpha\_{\rm s}(\mu\_{R,\langle i\rangle}^2) \,. \tag{5.367}$$

This scale is also used for the strong coupling related to loop corrections or to the soft emissions below the jet-threshold encoded for instance in the terms HN . They are effectively ignored in the construction of the parton shower history, which relies on tree-level configurations only.

When performing a scale variation of the renormalization scale in order to assess the corresponding uncertainty, one must therefore modify the Born matrix element by a factor such that

$$\alpha\_{\rm s}^{n}(\mu\_{R}^{2}) \longrightarrow \alpha\_{\rm s}^{n}(\tilde{\mu}\_{R}^{2}) \left(1 - \frac{\alpha\_{\rm s}(\tilde{\mu}\_{R}^{2})}{2\pi} \beta\_{0} \sum\_{i=1}^{n} \log \frac{\mu\_{i}^{2}}{\tilde{\mu}\_{R}^{2}}\right),\tag{5.368}$$

while all higher order terms are just evaluated with the coupling taken at the new scale ˜µ 2 R. This is slightly different from the scale setting prescription in MINLO, where the value of αs chosen for the additional real or virtual correction is fixed to be the average of all other values of the strong coupling, cf. Eq. (5.375).

In a similar way, a new factorization scale can be chosen for the matrix elements, but must be compensated for by a term of the form

$$\mathcal{B}\_N \log \frac{\tilde{\mu}\_F^2}{\mu\_F^2} \left[ \sum\_{c=q,g} \int \frac{\mathrm{d}z}{z} \mathcal{P}\_{ac}(z) \, f\_{c/h\_a} \left( \frac{x\_a}{z}, \tilde{\mu}\_F^2 \right) + \sum\_{d=q,g} \int \frac{\mathrm{d}z}{z} \mathcal{P}\_{bd}(z) \, f\_{d/h\_b} \left( \frac{x\_b}{z}, \tilde{\mu}\_F^2 \right) \right]. \tag{5.369}$$

## 5.6 NNLO and parton showers

Building on the algorithms and technologies discussed so far, the first implementations of combinations of calculations at next-to-next-to leading order (NNLO) in the strong

#### 394 QCD to All Orders

coupling for simple processes — essentially the production of simple colour–singlet systems — and parton showers have emerged. The first of them, dubbed MINLO, builds on a scale–setting prescription similar to the one employed in MEPS@NLO; its alternative, UNNLOPS in turn employs technology imported from the UNLOPS approach to multijet merging of NLO-QCD matrix elements. Both will briefly be sketched in the remainder of this chapter.

#### 5.6.1 Multijet merging without scales: MINLO

#### <IP_ADDRESS> Underlying ideas

The MINLO ideas have been worked out in [611] (MINLO-1) and [609] (MINLO-2).

The first instance focuses on improving the behaviour of fixed-order calculations in the soft region of parton emission, where they usually show unphysical, i.e. diverging, behaviour. This is achieved by a combination of clever choices for the scale of the strong coupling and the modification of the matrix elements by Sudakov form factors, similar to the programme already introduced in multijet merging at leading order, cf. Section 5.5.1. The main complication, though, that needed to be addressed stems from the fact that in contrast to the leading-order matrix elements, on which the CKKW method [351], operated, MINLO-1 aims at modifying next-to-leading order matrix elements.

MINLO-2 pushes the ideas in MINLO-1 one step further, by reweighting the matrix elements with Sudakov form factors not at next-to-leading logarithmic accuracy (NLL) but actually at the next-to-next-to-leading logarithmic accuracy (NNLL). This ingenious idea allows the transverse momentum of the softest parton of the Bornlevel configuration to include scales of the order the parton shower cut-off tc, around 1 GeV. In turn MINLO-2 feeds such configurations directly into the parton shower through the POWHEG method, thereby essentially taking into account NLO corrections for processes such as Drell–Yan or Higgs production in association with one potentially very soft parton in a full hadron-level simulation.

#### <IP_ADDRESS> MINLO-1

As already stated, the primary aim of MINLO-1 is to improve the behaviour of single fixed–order calculations at NLO accuracy in the strong coupling in the Sudakov region, without upsetting their formal NLO accuracy. This is an important step towards improved stability of NLO calculations where a heavy system such as, e.g., a gauge or a Higgs boson is produced in association with light jets. In such cases, the stability of the calculation suffers with decreasing transverse momenta of the heavy system, or, correspondingly, of the light jets, due to the emergence of increasingly large logarithms, which must be resummed to all orders. As already encountered in previous sections, a convenient way to achieve this stabilization is through inclusion of now familiar Sudakov form factors. It should thus not be a big surprise that successful algorithms achieving this aim, such as the MINLO-1 method, are further developing ideas imported from multijet merging, presented in Section 5.5.3, but usually without the additional intricacies related to the combination of various NLO calculations for increasing multiplicities.

There are two ingredients in MINLO-1, already familiar from multijet merging. First of all, the renormalization and factorization scales are chosen to capture better the interplay of various scales in the process, and, secondly, the introduction of Sudakov weights which tame the instabilities emerging in the regime of small transverse momenta. This proceeds by constructing a parton shower history for the underlying Born-level configuration, omitting the softest emission in the real-emission contribution. As before, a strong ordering of the N nodal values in the hardness scale for the emissions is assumed,

Q0 ≥ Q1 ≥ Q2 . . . ≥ QN−1 ≥ QN = Qcut , (5.370)

where Q0 denotes the hardest scale of the core process. For instance, in the case of the production of a heavy singlet in association with a number of jets, Q0 typically is the mass of the singlet. The weight corresponding to the individual phase-space point is multiplied with an overall suppression factor, given by a product of Sudakov form factors of all internal lines i of the Born configuration and a product of all Sudakov form factors of all outgoing partons k of the Born configuration, which have been produced through a splitting at the scale Qk,

$$\mathcal{S} = \prod\_{i=1}^{N} \Delta\_i(Q\_{i-1}^2, Q\_i^2) \prod\_k \Delta\_k(Q\_k^2, Q\_{\text{cut}}^2) \,. \tag{5.371}$$

Here, the subscripts i and k in ∆i and ∆k denote the flavour of the internal and external lines and the analytic Sudakov form factors from Eq. (2.183) are employed. The two external partons emerging at the last branching, i.e. at scale QN , will be associated with a Sudakov factor of ∆k(Q2 cut, Q2 cut) = 1.

In order to compensate for higher-order effects induced by the Sudakov form factors, the Born-term is modified by a factor which captures the first order in the αsexpansion of the analytic Sudakov form factors above. For each of the Sudakov form factors, the corresponding first-order term in its expansion is denoted by ∆(1). With this notation, the correction factor is given by

$$\begin{split} 1 &= \sum\_{i} \Delta^{(1)}(Q\_{i-1}^{2}, Q\_{i}^{2}) - \sum\_{k} \Delta\_{k}^{(1)}(Q\_{k}^{2}, Q\_{\text{cut}}^{2}) \\ &= 1 + \sum\_{i} \int\_{Q\_{i}^{2}} \frac{\mathrm{d}q\_{\perp}^{2}}{q\_{\perp}^{2}} \frac{\alpha\_{s}(q\_{\perp}^{2})}{2\pi} \Gamma\_{i}(Q\_{i-1}^{2}, q\_{\perp}^{2}) + \sum\_{k} \int\_{Q\_{\text{cut}}^{2}} \frac{\mathrm{d}q\_{\perp}^{2}}{q\_{\perp}^{2}} \frac{\alpha\_{s}(q\_{\perp}^{2})}{2\pi} \Gamma\_{i}(Q\_{k}^{2}, q\_{\perp}^{2}), \end{split} \tag{5.372}$$

where the Γi are the integrated splitting functions introduced in Eq. (2.178), and

$$\Delta\_k^{(1)}(Q^2, Q\_0^2) = -\int\_{Q\_{\rm cut}^2}^{Q\_k^2} \frac{\mathrm{d}q\_\perp^2}{q\_\perp^2} \frac{\alpha\_s(t)}{2\pi} \Gamma\_i(Q\_k^2, q\_\perp^2) \,. \tag{5.373}$$

#### 396 QCD to All Orders

Here, the first sum ranges over all internal lines i, while the second sum includes all outgoing parton of the Born configuration. This factor is in complete analogy to the one in Eq. (5.365), see also Section 5.5.3.

To see how the scale setting works out for the MINLO method, consider again a Born-level configuration with M powers of the strong coupling related to the core process and with m QCD emissions with corresponding scales µR,(i) . In tree-level multijet merging methods, the overall scale µR is implicitly given by

$$
\alpha\_{\rm s}^{M+N}(\mu\_R^2) = \alpha\_{\rm s}^M(\mu\_{R,\text{(core)}}) \cdot \prod\_{i \in N} \alpha\_{\rm s}(\mu\_{R,\{i\}}) \,, \tag{5.374}
$$

cf. Eq. (5.345).

It may seem natural to use this scale also as the argument of the strong coupling in the real or virtual correction, and in fact, this is the way multijet merging methods at NLO choose the scale, cf. Section 5.5.3. In slight contrast, in MINLO-1 the average value of αs,

$$\alpha\_{\rm s}^{(M+N+1)} = \frac{1}{M+n} \left[ M \alpha\_{\rm s} (\mu\_{R, \text{(core)}}^2) + \sum\_{i=1}^{N} \alpha\_{\rm s} (\mu\_{R, \text{(i)}}) \right],\tag{5.375}$$

is employed for them; the rationale for this choice is detailed in the publication [611].

#### <IP_ADDRESS> MINLO-2: first implementation of NNLOPS

Up to now, the MINLO-2 method has been formulated for the production of singlet systems only, for reason that will become obvious. The idea in it is to basically feed the next-to-leading order expressions in MINLO-1 for the production of such a singlet system S in association with a parton j into the POWHEG formalism, thereby providing the link to the parton shower. By allowing the parton j to become as soft or collinear as the parton shower cut-off tc allows, the full one-parton emission phase space is filled at NLO accuracy, and the second parton is distributed according to LO accuracy. This is possible by reweighting the Born-level singlet plus jet configuration with a Sudakov form factor at O α 2 s accuracy, i.e. including terms A2 and B2 from QT -resummation, thereby compensating all dangerous logarithms of the low scale. The catch then is to reweight the emergent event sample to a suitable distribution of the singlet at NNLO accuracy, to achieve the overall NNLO+PS accuracy.

Taking as an illustrative example for this procedure the production of a Higgs boson in conjunction with a parton, the B¯ term for H + j production can be written in a similar way as before, including the Born term plus the real and virtual corrections. Following the reasoning of MINLO-1 above, however, Sudakov form factors are invoked to account for higher-order corrections which become large for small transverse momenta, and their interplay with the genuine higher-order terms must be accounted for by subtracting out their respective first order expansion. Taken together, and making all factors of αs explicit,

NNLO and parton showers 397

$$
\begin{split}
\bar{\mathcal{B}}(\Phi\_{\mathcal{B}}) = \alpha\_{\text{s}}(m\_{H}^{2})\,\alpha\_{\text{s}}(q\_{\perp})\,\Delta\_{\text{g}}^{2}(m\_{H}^{2},Q\_{\perp}^{2})\cdot\left\{\mathcal{B}(\Phi\_{\mathcal{B}})\left[1-2\Delta^{(1)}(m\_{H}^{2},Q\_{\perp}^{2})\right]\right. \\ \\ \left. + \alpha\_{\text{s}}(Q\_{\perp}^{2})\left[\tilde{\mathcal{V}}(\Phi\_{\mathcal{B}})+\int\mathrm{d}\Phi\_{1}\,\mathcal{R}(\Phi\_{\mathcal{B}}\times\Phi\_{1})\right]\right\}, \end{split} \tag{5.376}
$$

where Q⊥ is the transverse momentum of the Higgs boson (and therefore, at Bornlevel, the extra parton) with Q2 ⊥ ≥ tc.

Integrating Eq. (5.376) over the full phase space of the Born-level H +j configuration will therefore, by construction, yield the cross-section for H +j at NLO accuracy, including all terms that are singular in the limit Q⊥ → 0, giving rise to logarithms of Q⊥. All singular terms, formally up to O α 2 s with respect to the inclusive Higgs boson production, can of course also be obtained through the QT -resummation formalism at NNLL accuracy, which recovers all these logarithms. These are all terms of the form αsL 2 and αsL as well as α 2 sL 4 , α 2 sL 3 , α 2 sL 2 , and α 2 sL, where the shorthand notation L = log(Q2/Q2 ⊥) has been used. In QT resummation the resummation scale Q typically is identified to be of order of the singlet mass, in this case therefore Q = O (mH).

Therefore, including also the A2 and B2 terms into the reweighting Sudakov form factors above, and fixing the scale of αs in the real and virtual contributions to Q2 ⊥ guarantees that the result not only is NLO-accurate for H + j, but also for inclusive H production.

There is one minor problem remaining, though, namely that the Sudakov form factors and coefficients given in Section 5.2.1, Eq. (5.65), Eq. (5.66), and Eq. (5.72), are for QT resummation in the conjugate b⊥-space. Here, however, the Sudakov form factors are meant to be directly applied in transverse momentum space, and therefore a translation between both must be applied. As noted in [609], this actually has been worked out in [505] and essentially leads to adding a term

$$
\Delta B\_{2,b\perp \to q\perp}^{(q,g)} = 4\zeta(3) \left( A\_1^{(q,g)} \right)^2,\tag{5.377}
$$

with A (q,g) 1 the corresponding soft coefficient, while the A2 remain unaltered. The B2 coefficients in the case of resummation directly in QT -space therefore are given by

$$B\_{2,Q\_T}^{(q,g)} = B\_{2,b\_\perp}^{(q,g)} + \Delta B\_2^{(q,g)}.\tag{5.378}$$

These terms are used in the analytic NNLL Sudakov form factors, multiplied to the matrix element, of the MINLO method such that

$$\Delta\_k^{\text{(NNLL)}}(Q^2, Q\_0^2) = \exp\left\{-\int \frac{q^2}{q\_\perp^2} \left[A(q\_\perp^2)\log\frac{Q^2}{q\_\perp^2} + B(q\_\perp^2)\right] \right\},\tag{5.379}$$

where

398 QCD to All Orders

$$\begin{split} A(q\_{\perp}^{2}) &= \frac{\alpha\_{\text{s}}(q\_{\perp}^{2})}{2\pi} A\_{1} + \left(\frac{\alpha\_{\text{s}}(q\_{\perp}^{2})}{2\pi}\right)^{2} A\_{2} \\ B(q\_{\perp}^{2}) &= \frac{\alpha\_{\text{s}}(q\_{\perp}^{2})}{2\pi} B\_{1} + \left(\frac{\alpha\_{\text{s}}(q\_{\perp}^{2})}{2\pi}\right)^{2} B\_{2,Q\_{T}}, \end{split} \tag{5.380}$$

with the coefficients A from Eq. (5.65), B1 from Eq. (5.67), and the term B2 from Eq. (5.378) with the B2 term in b⊥-space given in Eq. (5.72). Note that the latter also include finite loop-correction terms, making them process-dependent.

Finally, with both the inclusive process — in the example, inclusive Higgs production — and the process with an additional jet — H+ jet production — NLO correct, a seamless merging of the two multiplicities has been successfully achieved. This happened, and it is important to stress this, without any jet cut like in the usual multijet merging. Instead, by carefully adjusting Sudakov weights and scales in αs it was possible to modify the H+ jet part in such a way that it automatically also accounts for the inclusive part.

In [610] it was then further realized that this could be turned into a full partonlevel simulation, but accurate not only at NLO for both H and H+ jet production, but in fact accurate at NNLO in the inclusive cross-section. In order to achieve this, the rapidity distribution of the Higgs boson in the seamless merged sample merely has to be reweighted to the NNLO distribution. The same algorithm was also applied in a follow-up study, concerning Drell–Yan production [652]. There, the reweighting proceeds in three dimensions, thereby capturing the dynamics of the lepton system.

#### 5.6.2 An alternative NNLO+PS implementation: UNNLOPS

To see how the UNNLOPS method [633, 634] works, consider first the simpler NLO case, known as UNLOPS. The idea there is to decompose the radiation pattern of all emissions into one part, where no emission happens and a complementary part capturing all emissions, and in particular at least one. Omitting the by-now obvious arguments of phase space for the different fixed-order contributions, the expectation value of an observable O is given by

$$
\begin{split}
\langle\mathcal{O}\rangle &= \left\{ \int \mathrm{d}\Phi\_{\mathcal{B}} \left[ \tilde{\mathcal{B}}\_{N} - \int \mathrm{d}\Phi\_{1} \mathcal{R}\_{N} \right] + \int\_{t\_{c}} \mathrm{d}\Phi\_{\mathcal{R}} \left[ 1 - \Delta\_{N}(t\_{1}, \mu\_{Q}^{2}) \right] \, \mathcal{R}\_{N}(\Phi\_{\mathcal{R}}) \right\} \mathcal{O}(\Phi\_{\mathcal{B}}), \\ &\quad + \int\_{t\_{c}} \mathrm{d}\Phi\_{\mathcal{R}} \, \Delta\_{N}(t\_{1}, \mu\_{Q}^{2}) \, \mathcal{R}\_{N} \mathcal{E}\_{N}^{(\mathcal{K})}(t\_{1}, t\_{c}; \mathcal{O}) \,.
\end{split}
$$

Here, B¯ is the differential NLO cross-section for Born-level kinematics from Eq. (5.326),

$$\mathcal{B}\_N(\Phi\_\mathcal{B}) = \mathcal{B}(\Phi\_\mathcal{B}) + \tilde{\mathcal{V}}\_N(\Phi\_\mathcal{B}) + \int \mathrm{d}\Phi\_1 \left[ \mathcal{R}\_N(\Phi\_\mathcal{B} \otimes \Phi\_1) - \mathcal{S}\_N(\Phi\_\mathcal{B} \otimes \Phi\_1) \right] . \tag{5.381}$$

The argument of the observable indicates the particle composition and phase space from which it is evaluated. In view of this, the first line contains the observable after no emission has taken place — the first term is the fixed-order result for a Born-level configuration where all fixed-order parton emissions above tc in the real correction have been subtracted, such that the one potentially remaining emission is below the parton shower cut-off. The second term in the first line refers to those events, which, starting from a real correction with a parton above the parton shower cut-off where vetoed due to an unwanted emission harder than t1. In this case the parton configuration is projected onto the underlying Born configuration. The second line refers to those events which starting from a real-emission configuration experienced a parton shower evolution without unwanted emissions. This is signalled by the parton shower evolution operator E (K) n (t, tc; O), defined in full analogy to the Eq. (5.310),

$$\begin{split} \mathcal{E}\_{n}^{(\mathsf{K})}(t,\mathsf{t}\_{\mathsf{c}};\mathcal{O}) &= \Delta\_{n}^{(\mathsf{K})}(t,\mathsf{t}\_{\mathsf{c}})\mathcal{O}(\Phi\_{n}) \\ &+ \int\_{t\_{\mathsf{c}}}^{t} \mathrm{d}\Phi\_{1} \left[ \mathcal{K}\_{n}(\Phi\_{1}) \, \Delta\_{n}^{(\mathsf{K})}(t,\mathsf{t}(\Phi\_{1})) \otimes \mathcal{E}\_{n+1}^{(\mathsf{K})}(t(\Phi\_{1}),\mathsf{t}\_{\mathsf{c}};\mathcal{O}) \right] . \end{split}$$

While the expression in Eq. (5.381) captures all relevant terms to maintain fixedorder accuracy, it does not account for the dependence of the real-emission contribution RN on the renormalization and factorization scales. This can trivially remedied, however, by suitable weights, see also the original literature.

The same logic has been extended to match the parton shower to NNLO calculations for the production of colour singlets, performed with qT subtraction. There, contributions with exactly zero, exactly one, and more than one emissions above tc have been identified, leading to a somewhat lengthy expression, for which we refer to the literature.

# 6 Parton Distribution Functions

Parton Distribution Functions (PDFs) are a necessary ingredient in the calculation of particle cross-sections at collider experiments with hadron beams.

In Chapter 2, the master formula Eq. (2.52)

$$\sigma \ll = \sum\_{a,b} \int\_0^1 \mathrm{d}x\_a \mathrm{d}x\_b \int f\_{a/h\_1}(x\_a, \mu\_F) f\_{b/h\_2}(x\_b, \mu\_F) \,\mathrm{d}\hat{\sigma}\_{ab \to n}(\mu\_F, \mu\_R) \tag{6.1}$$

was introduced to describe the calculation of the production of n-parton final states in processes of the type h1h2 → n + X with incoming hadrons h1 and h2 – usually protons. The calculation of such a hard (parton-level) cross-section within perturbative QCD relies on also using partons in the inital state and thus requires the knowledge of the distribution of partons i in the hadrons h. This is most often achieved by using a scheme called collinear factorization in which this distribution depends on the longitudinal momentum fraction x of the partons with respect to the hadron, taken at the factorization scale µF . The resulting distributions fi/h(x, µF ) are known as parton distribution functions (PDFs), and parameterize the transition of incident hadrons to incident partons, thereby absorbing all emissions at scales below µF , cf. Chapter 2. In turn, the PDFs obviously depend upon this measure of the hardness of the parton–level interaction, which will be of the same order as the renormalization scale µR. As a consequence of the factorization of possible parton emissions into a soft and collinear part, below µF , and a hard part, above µF , the latter are explicitly described by the matrix element. Factorization theorems, proven for deep–inelastic lepton–proton scattering and for Drell-Yan production of gauge bosons in hadron collisions, assert that the PDFs in collinear factorization are process–independent. There are however terms of higher dimension – so called higher–twist contributions – that are process dependent. These terms are ignored throughout this book, since they are suppressed by factors of the type m2 p/µ2 F . Note that in this chapter both µR and µF will often be represented just as Q.

Pictorially, as was shown in Fig. 2.5 (left), more of the quantum fluctuations can be resolved as the hardness scale, corresponding to the inverse of the time scale, increases. Low–scale processes can only resolve longer time intervals, and with increasing scale,

*The Black Book of Quantum Chromodynamics: A Primer for the LHC Era.* John Campbell, Joey Huston, and Frank Krauss. © John Campbell, Joey Huston, and Frank Krauss 2018. Published in 2018 by Oxford University Press. DOI 10.1093/oso/9780199652747.001.0001

smaller time intervals are being probed, and more of the quantum fluctuations inside the hadron are resolved, which transfer momentum from high–x partons to low–x partons. Thus, at higher µF one expects to find more of the momentum of the proton given to gluons and sea quarks with relatively low values of the momentum fraction x, as the quantum fluctuations creating these partons can be resolved. Conversely, the population of partons containing a large fraction of the parent hadron's momentum will decrease.

The PDFs describing this dynamics are treated as being universal, i.e. they can be determined by one set of processes and used to predict the cross-section for any process, using the factorized form shown above. These parton distribution functions currently can not be calculated perturbatively; ultimately, however, it may become possible in the future to calculate these PDFs non-perturbatively, using lattice gauge theory [310]. On the other hand, the evolution of PDFs with Q2 can be calculated with a perturbative treatment using the DGLAP evolution equations, as sketched in Section 2.1.3 and discussed in more detail in Section 6.1.

PDFs have been determined by global fits to data from a plethora of different data. Modern global PDF fits use of the order of 3000 data points from a number of processes from fixed target experiments, the TEVATRON, HERA, and the LHC. In particular, data from deep-inelastic scattering (DIS), Drell-Yan (DY) and jet production processes, have played a dominant role in the past. With the advent of the LHC and its huge amount of data, also processes such as γ+jet, W+jets, and heavy flavour production (including tt¯) have become increasingly important as the statistical and systematic errors of the data sets have improved [827]. Many of these processes have recently been calculated to NNLO, allowing the data to be used in PDF fits at that order. At the moment this endeavour is actively pursued by a number of different groups: ABM [136], CTEQ/CT [489, 551], HERAPDF [115, 413, 819], JR [649], MSTW/MMHT [614, 755] and NNPDF [192, 194], which provide semi-regular updates to their fits of parton distributions, when new data and/or theoretical developments become available. Similarities and differences in the fitting procedure performed by the different groups will be discussed in some detail in Section 6.2.

The most commonly used are those PDFs which perform a global analysis on a broad variety of data using a variable flavour number scheme; the most recent updates (CT14, MMHT2014, NNPDF3.0) are given in [194, 489, 614]. The resulting PDFs are available at leading order (LO), next-to-leading order (NLO), and next-to-next-toleading order (NNLO) in the strong coupling constant αS, depending on the order(s) at which the global PDF fits have been carried out. Some PDFs have also been produced at what has been termed modified leading order [715, 846] (LO∗ or similar), in an attempt to reduce some of the problems that result from the use of LO PDFs in parton shower Monte Carlo programs. These PDFs are no longer in wide use, however. The choices of parameterization for these PDFs are then discussed, along with the impact of the use of particular renormalization and factorization schemes. Given the wide kinematic coverage of the data, the parameterization of the PDFs must be flexible enough to describe the parton distributions over a wide range of x and Q2 , and to not introduce artificial correlations between different x regions. Modern PDFs take into account the finite charm and bottom quark masses (using a variety of heavy quark

#### 402 Parton Distribution Functions

mass schemes) in their fits to data, particularly to data from deep-inelastic scattering. This has consequences for low-mass cross-sections at the LHC that will be examined. Technical aspects related to the different orders, schemes and parameterizations are highlighted in Section 6.2.2.

There are two different classes of technique employed for the PDF determination: those based on the Hessian approach, and those using a Monte Carlo approach. The background behind both classes will be discussed in this chapter, and some of the actual procedures employed in the fitting of the PDFs and the determination of the respective errors will be highlighted. An important consequence of the method chosen is the way uncertainties in the fits are handled. These uncertainties also have implications on the accuracy of overall cross-section calculations. This is followed by a similar discussion of the choice of the value of αs(mZ) in the global fits. The data in the global fits are mostly from strong interaction physics, and they are thus sensitive to the exact value of αs(mZ). As a consequence, the global fit itself can be used to determine it. Since αs(mZ) is a universal parameter, however, another approach is to assume the world average value, and to not let αs be a free parameter in the fit. Both approaches can be and have been used. Recently, most groups have provided PDFs at the world average value of αs(mZ), along with PDFs at alternate values of αs(mZ). Issues related to the fitting technology, impacting on estimates of the intrinisic uncertainties of PDFs and their impact on cross-sections are discussed in Section 6.3.

The next section is devoted to a discussion of PDF correlations. The examination and understanding of correlations between two fitted PDFs, or between a PDF and a cross-section, or between two cross-sections is crucial to a detailed appreciation of PDFs at the LHC. In adition, such correlations can be used to decrease the PDF uncertainties for the ratios of such quantities at the LHC. In Section 6.4, the resultant PDFs at LO (modified LO), NLO, and NNLO are presented. Parton luminosities are defined and are shown for several important initial states at the LHC in Section 6.5. PDF uncertainties using inputs from several PDF groups are then defined, using the PDF4LHC accords. Correlations among LHC processes are examined using these PDF combinations, and finally, in Section 6.6 several useful PDF tools that are currently available are described.

## 6.1 PDF evolution: the DGLAP equation revisited

#### 6.1.1 PDF evolution at leading order

An essential feature of the PDFs is the fact that their evolution with energy can be calculated within the framework of perturbative QCD, as has already been outlined in Chapter 2. In that chapter the equation governing this evolution – the QCD DGLAP equation – was first introduced. At the leading order it reads,

$$\begin{split} &\frac{\partial}{\partial\log Q^2} \left( \frac{f\_{q/h}(x,Q^2)}{f\_{g/h}(x,Q^2)} \right) \\ &= \frac{\alpha\_s(Q^2)}{2\pi} \int\_x^1 \frac{dz}{z} \left( \frac{\mathcal{P}\_{qq}^{(1)}\left(\frac{x}{z}\right)}{\mathcal{P}\_{gg}^{(1)}\left(\frac{x}{z}\right)} \frac{\mathcal{P}\_{gg}^{(1)}\left(\frac{x}{z}\right)}{\mathcal{P}\_{gg}^{(1)}\left(\frac{x}{z}\right)} \right) \begin{pmatrix} f\_{q/h}(z,Q^2) \\ f\_{g/h}(z,Q^2) \end{pmatrix}, \end{split} \tag{6.2}$$

The leading order splitting functions P (1) ij have previously been given in Eq. (2.33). They represent the kernel of this evolution equation and clearly couple different PDFs together. Note that the symmetries of QCD mean that, at this order, splitting functions involving anti-quarks are simply equal to their quark counterparts. Before turning to the form of the DGLAP equation and the splitting kernels at higher orders in the next section, the calculation of the kernels Pij at leading order will be sketched below.

#### <IP_ADDRESS> Factorization of initial state emissions in the collinear limit

The starting point for the calculation of the splitting kernels is the collinear limit of particle emission off initial-state partons. In order to see how this works, assume the production of a system X with mass MX in the scattering of two incident partons j and j 0 with momenta p and p 0 ,

$$p^{\mu} = \frac{\sqrt{\hat{s}}}{2} \begin{pmatrix} 1, 0, 0, 1 \end{pmatrix} \quad \text{and} \quad p'^{\mu} = \frac{\sqrt{\hat{s}}}{2} \begin{pmatrix} 1, 0, 0, -1 \end{pmatrix} \,,\tag{6.3}$$

where ˆs = (p + p 0 ) 2 . The differential parton–level cross-section for this process, at leading order, is

$$\mathrm{d}\hat{\sigma}^{(0)}\_{jj'\to X}(p,p') = \frac{1}{2\hat{s}}\,\mathrm{d}\Phi\_X\left|\mathcal{M}^{(0)}\_{jj'\to X}(p,p')\right|^2. \tag{6.4}$$

Assume now that parton j emits a parton k with momentum k into the final state, turning into a parton i with reduced four momentum p − k, with (p − k) 2 < 0, which together with j 0 produces X. The momentum k is then given by

$$k^{\mu} = \left(1 - x\right)p^{\mu} + \beta \left.p^{\prime \mu} + k\_{\perp}^{\mu}\right.\tag{6.5}$$

where k⊥ is perpendicular to both p and p 0 . β is fixed by the on-shell condition for k, k 2 = 0, and hence

$$k^{\mu} = \left(1 - x\right)p^{\mu} + \frac{\mathbf{k}\_{\perp}^{2}}{(1 - x)\hat{s}}p^{\prime \mu} + k\_{\perp}^{\mu}.\tag{6.6}$$

Here k 2 ⊥ denotes the positive square of the transverse momentum. Similarly,

$$q^2 = (p-k)^2 = \frac{\mathbf{k}\_{\perp}^2}{1-x},\tag{6.7}$$

showing that the propagator of i is proportional to k 2 ⊥. As a consequence, the corresponding matrix element for the process under emission of k diverges with k⊥ → 0. This is the collinear limit, where the transverse momentum of the emitted particle vanishes. This singular behaviour also triggers divergences in the phase space integration over k. All these divergences however cancel, according to the Kinoshita–Lee-Nauenberg theorem, and as seen previously, when virtual corrections are added.

#### 404 Parton Distribution Functions

To proceed, the phase space integral over particle k has to be cast into a form useful for further consideration. With a bit of algebra it is easy to show that

$$\mathrm{d}\Phi\_k = \frac{\mathrm{d}^4 k}{(2\pi)^4} (2\pi) \delta(k^2) = \frac{1}{(2\pi)^3} \mathrm{d}x \mathrm{d}\beta \mathrm{d}^2 \mathbf{k}\_\perp \,\delta\left( (1-x)\beta - \frac{\mathbf{k}\_\perp^2}{\hat{s}} \right)$$

$$= \frac{\mathrm{d}x \mathrm{d} \mathbf{k}\_\perp^2}{16\pi^2 (1-x)}.\tag{6.8}$$

Adding in the emission of the final state parton k, and ignoring all other possible emissions, yields a higher–order contribution to the cross-section of X production, and Eq. (6.4) becomes

$$\begin{split} \mathrm{d}\hat{\sigma}\_{jj' \to X}(p, p') &= \quad \frac{1}{2\hat{s}} \mathrm{d}\Phi\_X \left| \mathcal{M}\_{jj' \to X}^{(0)}(p, p') \right|^2 \\ &+ \frac{1}{2\hat{s}} \frac{1}{16\pi^2} \mathrm{d}\Phi\_X \int\_0^1 \frac{\mathrm{d}x}{1-x} \int \mathrm{d}^2 \mathbf{k}\_\perp^2 \left| \mathcal{M}\_{jj' \to ikj' \to kX}^{(0)}(p, p') \right|^2 . \end{split} \tag{6.9}$$

The leading-order splitting functions P (1) ji (x) in collinear factorization are then defined as the part of the second line of the equation above that diverges in the collinear limit:

$$\frac{1-x}{x}\frac{\alpha\_s}{2\pi}P\_{ji}^{(1)}(x)\left|\mathcal{M}\_{ij'\to X}^{(0)}(xp,p')\right|^2 = \frac{1}{16\pi^2}\lim\_{\mathbf{k}\to\pm 0} \left[\mathbf{k}\_{\perp}^2 \left|\mathcal{M}\_{jj'\to ikj'\to kX}^{(0)}(p,p')\right|^2\right].\tag{6.10}$$

The x in the denominator of the left-hand side of the equation arises from the fact that the matrix element for ij0 → X has an incoming flux given by 2xsˆ instead of the original 2ˆs.

#### <IP_ADDRESS> Calculating the kernels: P (1) qq

To calculate the splitting function P (1) qq , i and j are quark lines and k is a gluon. Clearly, gauge invariance of the overall matrix element for the associated production of X and the gluon k is only guaranteed if amplitudes for the emission of the gluon off all coloured lines are included. The question is whether these contributions exhibit collinear divergent behaviour.

It can be shown that this is not the case when choosing a smart physical gauge where the gluon's polarization vector µ is transverse to both p 0 and k. To see this, an additional axis n⊥ ⊥ k⊥ with n 2 ⊥ = 1 is introduced, which allows the explicit construction of the polarization vectors as

$$
\epsilon^{\mu}\_{\pm} = \frac{\sqrt{2}\,\mathbf{k}\_{\perp}}{(1-x)\hat{s}} \, p^{\prime \mu} + \frac{1}{\sqrt{2}\,\mathbf{k}\_{\perp}} \, k^{\mu}\_{\perp} \pm i n^{\mu}\_{\perp} \,. \tag{6.11}
$$

In addition to the orthogonality requirements underlying their construction,

PDF evolution: the DGLAP equation revisited 405

$$
\epsilon \cdot p' = \epsilon \cdot k = 0 \,. \tag{6.12}
$$

they satisfy

$$
\epsilon\_{\pm} \cdot k\_{\perp} = \frac{\mathbf{k}\_{\perp}}{\sqrt{2}} \quad \text{and} \quad \epsilon\_{\pm} \cdot p = \frac{\mathbf{k}\_{\perp}}{\sqrt{2}(1-x)}.\tag{6.13}
$$

While the square of the graph with the propagator of i boasts two propagators and is thus proportional to 1/k4 ⊥, all other graphs have at best one of these propagators. And since in the squared amplitude there are two scalar products involving the polarization vector, the numerator of the squared amplitude is proportional to k 2 ⊥. This leaves no overall divergence in 1/k2 ⊥ the squared amplitude apart from the one stemming from the graph with the i-propagator.

Analysing the structure of the reduced amplitude M(0) ij0→X, it becomes clear that its square can be written as

$$\left|\mathcal{M}\_{ij'\to X}^{(0)}\right|^2 = \bar{u}(xp,\lambda)\left[\gamma\_\mu M^\mu\right]u(xp,\lambda)\,,\tag{6.14}$$

where the spinors refer to the intermediate quark i. Because it is a massless quark, and since incoming and outgoing spinors come with the same helicity, this is the only allowed structure. Decomposing M according to

$$M^{\mu} = ap^{\mu} + bp^{\prime \mu} + m\_{\perp}^{\mu} \tag{6.15}$$

with m⊥ a vector in the same transverse plane as k⊥ and using the Dirac equatoin elimiates the term proportional to p. Also, since the helicities are the same, the spinflip operator γ · m vanishes when sandwiched between the two spinors. In the end, therefore

$$\left|\mathcal{M}\_{ij'\rightarrow X}^{(0)}\right|^2 = \bar{u}(xp,\lambda)\left[\gamma\_\mu bp^{\prime\mu}\right]u(xp,\lambda)\,. \tag{6.16}$$

realizing that

$$\left|\bar{u}(xp,\lambda)\left(\gamma \cdot p'\right)u(xp,\lambda)\right|^2 = \left(2xpp'\right)^2 = x^2\hat{s}^2\tag{6.17}$$

identifies b and therefore

$$M^{\mu} = \frac{p^{\prime \mu}}{2x\hat{s}} \left| \mathcal{M}\_{ij^{\prime} \to X}^{(0)} \right|^2. \tag{6.18}$$

These identities can now be rolled out to the calculation of the emission matrix element in the collinear limit. Including factors of 1/2 and 1/Nc for the average over incoming quark helicities and colours,

$$\begin{split} & \left| \mathcal{M}\_{jj' \to ikj' \to kX}^{(0)} (p, p') \right|^2 \\ &= \frac{g^2 \text{Tr} \left[ T^a T^a \right]}{2N\_c} \sum\_{\lambda, \pm} \bar{u}(p, \lambda) \not{e}\_{\pm} \frac{\not{p} - \not{k}}{(p - k)^2} \left[ \frac{\not{p}'}{x\hat{s}} \left| \mathcal{M}\_{ij' \to X}^{(0)} \right|^2 \right] \frac{\not{p} - \not{k}}{(p - k)^2} \not{\epsilon}\_{\pm}^\* u(p, \lambda) \\ &= \frac{g^2 \text{Tr} \left[ T^a T^a \right]}{2N\_c x \hat{s} (-2p \cdot k)^2} \sum\_{\pm} \text{Tr} \left[ \not{p} \not{\epsilon}\_{\pm} \left( \not{p} - \not{k} \right) \not{p}' \left( \not{p} - \not{k} \right) \not{\epsilon}\_{\pm}^\* \right] \left| \mathcal{M}\_{ij' \to X}^{(0)} \right|^2 \end{split}$$

406 Parton Distribution Functions

$$=\frac{2g^2C\_F}{x\mathbf{k}\_{\perp}^2}\left(1+x^2\right)\left|\mathcal{M}\_{ij'\to X}^{(0)}\right|^2. \tag{6.19}$$

Here, the following identities have been used:

$$2kp = \frac{\mathbf{k}\_{\perp}^{2}}{1 - x}$$

$$2kp' = \hat{s}(1 - x)$$

$$\sum\_{\pm} \epsilon\_{\mu}^{\pm} \epsilon\_{\nu}^{\*\pm} = -g\_{\mu\nu} + \frac{k\_{\mu}p'\_{\nu} + k\_{\nu}p'\_{\mu}}{kp'}.\tag{6.20}$$

Inserting this into Eq. (6.10) yields

$$\frac{1-x}{x} \frac{\alpha\_{\rm s}}{2\pi} P\_{qq}^{(1)}(x) = \frac{1}{16\pi^2} \times \frac{2g^2 C\_F}{x} \left(1+x^2\right) = \frac{8\pi\alpha\_{\rm s}C\_F(1+x^2)}{16\pi^2 x} \tag{6.21}$$

and therefore

$$P\_{qq}^{(1)}(x) = C\_F \frac{1+x^2}{1-x} \,, \tag{6.22}$$

as expected from Eq. (2.33).

#### <IP_ADDRESS> Virtual contributions: The "+"-function and all that

Having arrived at the expression in Eq. (6.22) is not the end of the story, though. Taking a closer look at the expression for P (1) qq shows the by now well–anticipated divergent behaviour for x → 1, or, for vanishing gluon energies. In this limit of course also the transverse momentum vanishes, leading, ultimately, to the usual double soft and collinear divergence. Of course, these divergences will cancel, when taking into account the virtual contributions, a feature that consistently appeared throughout the book. Formalizing this idea results in

$$\text{finite} = \int\_0^1 \text{d}x \left[ \frac{\alpha\_s}{2\pi} P\_{qq}^{(1)}(x) \left| \mathcal{M}\_{qj' \to X}(xp) \right|^2 + (1+V) \left| \mathcal{M}\_{qj' \to X}(xp) \right|^2 \right]$$

$$= \int\_0^1 \text{d}x \left[ \frac{\alpha\_s}{2\pi} P\_{qq}^{(1)}(x) + \delta(1-x)(1+V) \right] \left| \mathcal{M}\_{qj' \to X}(xp) \right|^2,\tag{6.23}$$

where V is the virtual contribution at order αs. This contribution can be obtained by invoking a simple probabilistic idea, similar to what has been seen in the construction of parton showers. To order αs, the probability of a gluon emission by the quark plus the probability of no emission must add to unity. Hence,

$$\int\_0^1 \mathrm{d}x \left[ \frac{\alpha\_s}{2\pi} P\_{qq}^{(1)}(x) + \delta(1-x) \left(1+V\right) \right] = 1\tag{6.24}$$

and therefore the virtual part must compensate the real emission encoded in the splitting function:

$$\begin{split} V &= -\frac{\alpha\_s}{2\pi} \lim\_{\delta \to 0} \int\_0^{1-\delta} \mathrm{d}x \, P\_{qq}^{(1)}(x) \\ &= \frac{\alpha\_s C\_F}{2\pi} \lim\_{\delta \to 0} \int\_0^{1-\delta} \mathrm{d}x \left( \frac{2}{1-x} - (1+x) \right) = \frac{\alpha\_s C\_F}{2\pi} \left( \frac{3}{2} + 2 \log \delta \right) . \end{split} \tag{6.25}$$

This implies that the complete splitting function is given by

$$\mathcal{P}\_{qq}(x) = C\_F \lim\_{\delta \to 0} \left[ \frac{1+x^2}{1-x} \Theta(1-x-\delta) + \delta(1-x) \left(\frac{3}{2} + 2\log\delta\right) \right] \tag{6.26}$$

This looks quite cumbersome. A solution, however, presents itself by realizing that the "splitting functions" Pij actually are not really functions but in fact distributions that will always act on some other functions f(x) such as PDFs or similar. These functions are typically regular at x = 1 and it is therefore meaningful to define a way to regulate the divergent behaviour of the Pij at x = 1 independent of the limitprocedure above. Introducing the "+" prescription, applicable for distributions that diverge at x = 1, such as 1/(1 − x),

$$\left(\frac{1}{1-x}\right)\_+ f(x) \stackrel{!}{=} \frac{f(x) - f(1)}{1-x},\tag{6.27}$$

cf. Eq. (2.10). This prescription will take care of the log δ-term, located at x = 1 by virtue of the δ-function and therefore, finally,

$$\mathcal{P}\_{qq}(x) = C\_F \left[ \left( \frac{1+x^2}{1-x} \right)\_+ + \frac{3}{2} \delta(1-x) \right]. \tag{6.28}$$

This actually implies that the x–integral over Pij vanishes exactly,

$$\int\_0^1 \mathrm{d}x \, \mathcal{P}\_{qq}(x) = \begin{array}{c} \text{0} \\ \text{1} \end{array} \tag{6.29}$$

thereby incorporating the exact cancellation of real and virtual contributions.

#### 6.1.2 PDF evolution at higher orders

PDF evolution at the next order in perturbation theory is achieved by expanding the splitting functions as a series in the strong coupling,

$$\mathcal{P}\_{ij}(x) = \mathcal{P}\_{ij}^{(1)} + \frac{\alpha\_{\text{s}}}{2\pi} \mathcal{P}\_{ij}^{(2)} + \mathcal{O}(\alpha\_{\text{s}}^2) \tag{6.30}$$

Fig. 6.1 Processes corresponding to the NLO splitting functions P (2) qiqj (left) and P (2) qiq¯j (right). Unobserved partons are indicated in the figure by parentheses.

The inclusion of higher–order terms, P (2) ij , introduces several subtleties that did not arise at leading order. The first is that, at leading order, the flavour structure of the evolution is trivial: for two flavours of quark, qi and qj , the splitting function trivially vanishes unless i = j, i.e. P (1) qiqj ∝ δij . This is no longer true at the next order, as indicated in Fig. 6.1 (a). Contributions originating from the splitting of a virtual gluon into a quark-antiquark pair of different flavour, qi → qj (qiq¯j ), give rise to kernels that are no longer diagonal in flavour space. Moreover, the same process also gives rise to kernels that directly couple quarks and antiquarks, qi → q¯j (qiqj ), cf. Fig. 6.1 (b).

Therefore Eq. (6.2) must be generalized in order to account for effects at higher orders. The evolution can be written as,

$$\begin{split} &\frac{\partial}{\partial\log Q^{2}} \begin{pmatrix} f\_{\bar{q}\_{i}/h}(x,Q^{2}) \\ f\_{\bar{q}\_{i}/h}(x,Q^{2}) \\ f\_{g/h}(x,Q^{2}) \end{pmatrix} \\ &= \frac{\alpha\_{s}(Q^{2})}{2\pi} \int\_{x}^{1} \frac{\mathrm{d}z}{z} \begin{pmatrix} \mathcal{P}\_{q\_{i}q\_{k}}\left(\frac{x}{z}\right) \mathcal{P}\_{q\_{i}\bar{q}\_{k}}\left(\frac{x}{z}\right) \mathcal{P}\_{qg}\left(\frac{x}{z}\right) \\ \mathcal{P}\_{\bar{q}\_{i}q\_{k}}\left(\frac{x}{z}\right) \mathcal{P}\_{\bar{q}\_{i}\bar{q}\_{k}}\left(\frac{x}{z}\right) \mathcal{P}\_{\bar{q}\_{g}}\left(\frac{x}{z}\right) \\ \mathcal{P}\_{gq\_{k}}\left(\frac{x}{z}\right) \mathcal{P}\_{g\bar{q}\_{k}}\left(\frac{x}{z}\right) \mathcal{P}\_{gg}\left(\frac{x}{z}\right) \end{pmatrix} \begin{pmatrix} f\_{q\_{k}/h}(z,Q^{2}) \\ \bar{f}\_{\bar{q}\_{k}/h}(z,Q^{2}) \\ \bar{f}\_{g/h}(z,Q^{2}) \end{pmatrix}, \end{split} (6.31)$$

with all of the distributions coupled together. The presence of some amount of antiquarks in the parent quark means that there is no physically-meaningful separation into "valence" and "sea" contributions beyond leading order. However it is possible to identify combinations of quark PDFs that can be identified with the valence contribution as follows. The nature of the leading order splitting functions suggests that they be decomposed according to,

$$\begin{aligned} \mathcal{P}\_{q\_i q\_j} &= \delta\_{ij} \mathcal{P}\_{qq}^V + \mathcal{P}\_{qq}^S \, , \\ \mathcal{P}\_{q\_i \bar{q}\_j} &= \delta\_{ij} \mathcal{P}\_{q\bar{q}}^V + \mathcal{P}\_{q\bar{q}}^S \, . \end{aligned} \tag{6.32}$$

In order to solve the coupled evolution equations, Eq. (6.31), it is useful to introduce the following combinations of the splitting functions [550],

$$\begin{aligned} \mathcal{P}\_{(\pm)} &= \mathcal{P}\_{qq}^V \pm \mathcal{P}\_{q\bar{q}}^V \,, \\ \mathcal{P}\_{QQ} &= \mathcal{P}\_{(+)} + 2n\_f \mathcal{P}\_{qq}^S \end{aligned}$$

PDF evolution: the DGLAP equation revisited 409

$$\begin{aligned} \mathcal{P}\_{Qg} &= 2n\_f \mathcal{P}\_{qg} \\ \mathcal{P}\_{gQ} &= \mathcal{P}\_{gq} \end{aligned} \tag{6.33}$$

for nf flavours of light quarks. Similarly it is convenient to express the parton distributions in terms of the following quantities,

$$f\_{q\_i^{(\pm)}/h}(x,Q^2) = f\_{q\_i/h}(x,Q^2) \pm f\_{\bar{q}\_i/h}(x,Q^2) \,, \tag{6.34}$$

$$f\_{Q^{(\pm)}/h}(x,Q^2) = \sum\_{i=1}^{n\_f} f\_{q\_i^{(\pm)}/h}(x,Q^2) \tag{6.35}$$

that is, into either a sum or difference of quark and antiquark distributions. The evolution expressed in Eq. (6.31) can now be written in a much simpler form in which the components decouple significantly. Two combinations satisfy a particularly simple form of the evolution,

$$\frac{\partial}{\partial \log Q^2} f\_{q\_i^{(-)} / h}(x, Q^2) = \frac{\alpha\_s(Q^2)}{2\pi} \int\_x^1 \frac{dz}{z} \mathcal{P}\_{(-)} f\_{q\_i^{(-)} / h}(z, Q^2) \tag{6.36}$$

$$\frac{\partial}{\partial \log Q^2} S\_i(z, Q^2) = \frac{\alpha\_s(Q^2)}{2\pi} \int\_x^1 \frac{dz}{z} \mathcal{P}\_{(+)} S\_i(z, Q^2) \,, \tag{6.37}$$

where the quantity Si(x, Q2 ) is given by,

$$S\_i(x, Q^2) = n\_f f\_{q\_i^{(+)} / h}(x, Q^2) - f\_{Q^{(+)} / h}(x, Q^2) \,. \tag{6.38}$$

These evolution equations can be solved separately; the solutions for f q (−) i /h(x, Q2 ) and Si(x, Q2 ) are referred to as non-singlet contributions. The remaining quark singlet combination, fQ(+)/h(x, Q2 ), remains coupled to the gluon contributions, but in a fashion akin to the leading order evolution (cf. Eq. (6.2)),

$$\begin{split} & \frac{\partial}{\partial \log Q^2} \begin{pmatrix} f\_{Q^{(+)}/h}(x,Q^2) \\ f\_{g/h}(x,Q^2) \end{pmatrix} \\ &= \frac{\alpha\_s(Q^2)}{2\pi} \int\_x^1 \frac{dz}{z} \begin{pmatrix} \mathcal{P}\_{QQ}\left(\frac{x}{z}\right) \mathcal{P}\_{Qg}\left(\frac{x}{z}\right) \\ \mathcal{P}\_{gQ}\left(\frac{x}{z}\right) \mathcal{P}\_{gg}\left(\frac{x}{z}\right) \end{pmatrix} \begin{pmatrix} f\_{Q^{(+)}/h}(z,Q^2) \\ f\_{g/h}(z,Q^2) \end{pmatrix} . \end{split} \tag{6.39}$$

Eqs. (6.36), (6.37), and (6.39) are sufficient to determine the evolution of all of the individual PDFs at a given order, if the corresponding splitting functions are known to the same accuracy.

The first corrections to the leading-order splitting functions were originally computed in Refs. [421, 549]. Their forms are reproduced here for completeness. The twoloop gluon splitting function receives contributions from two different colour structures, 410 Parton Distribution Functions

$$\mathcal{P}\_{gg}^{(2)}(z) = C\_A \mathcal{P}\_{ggA}^{(2)}(z) + T\_R n\_f \mathcal{P}\_{ggF}^{(2)}(z) \,, \tag{6.40}$$

that are given by,

$$\begin{split} \mathcal{P}\_{ggA}^{(2)}(z) &= \frac{\Gamma\_1}{8} \frac{1}{C\_A} \left[ P\_{gg}^{(1)}(z) \right]\_+ + \delta(1-z) \left[ C\_A(-1+3\zeta\_3) + \beta\_0 \right] \\ &+ \left[ P\_{gg}^{(1)}(z) \right]\_+ \left( -2\ln(1-z) + \frac{1}{2}\ln z \right) \ln z + \left[ P\_{gg}^{(1)}(-z) \right]\_+ \left( S\_2(z) + \frac{1}{2}\ln^2 z \right) \\ &+ C\_A \left( 4(1+z)\ln^2 z - \frac{4(9+11z^2)}{3}\ln z - \frac{277}{18z} + 19(1-z) + \frac{277}{18}z^2 \right) \\ &+ \beta\_0 \left( \frac{13}{6z} - \frac{3}{2}(1-z) - \frac{13}{6}z^2 + (1+z)\ln z \right), \\ \mathcal{P}\_{ggF}^{(2)}(z) &= C\_F \left[ -\delta(1-z) + \frac{4}{3z} - 16 + 8z + \frac{20}{3}z^2 - 2(1+z)\ln^2 z - 2(3+5z)\ln z \right], \end{split} \tag{6.41}$$

where the plus-part of the first-order splitting functions has been defined previously in Eq. (2.33). This equation also introduces both the two-loop cusp anomalous dimension,

$$
\Gamma\_1 = \frac{4}{3} \left( C\_A (4 - \pi^2) + 5\beta\_0 \right) \tag{6.42}
$$

as well as the function,

$$\operatorname{Li}\_2(z) = -2\operatorname{Li}\_2(-z) - 2\ln(1+z)\ln z - \frac{\pi^2}{6} \, . \tag{6.43}$$

The corresponding result for P (2) gq (z) is,

$$\begin{split} \mathcal{P}\_{qq}^{(2)}(z) &= C\_{A} \left\{ P\_{qq}^{(1)}(z) \left[ \ln^2(1-z) - 2\ln(1-z)\ln z - \frac{101}{18} - \frac{\pi^2}{6} \right] + P\_{qq}^{(1)}(-z)S\_2(z) \\ &+ C\_{F} \left( 2z\ln(1-z) + (2+z)\ln^2 z - \frac{36+15z+8z^2}{3}\ln z + \frac{56-z+88z^2}{18} \right) \right\} \\ &- C\_{F} \left\{ P\_{qq}^{(1)}(z)\ln^2(1-z) + [3P\_{qq}^{(1)}(z) + 2zC\_{F}]\ln(1-z) \\ &+ C\_{F} \left( \frac{2-z}{2}\ln^2 z - \frac{4+7z}{2}\ln z + \frac{5+7z}{2} \right) \right\} \\ &+ \beta\_0 \left\{ \mathcal{P}\_{qq}^{(1)}(z) \left[ \ln(1-z) + \frac{5}{3} \right] + z \right\}. \end{split} \tag{6.44}$$

The quark splitting functions employ the decomposition of Eq. (6.32). The components are,

$$\mathcal{P}\_{qqV}^{(2)}(z) = \frac{\Gamma\_1}{8} C\_F \frac{(1+z^2)}{(1-z)\_+} $$

#### Fitting parton distribution functions 411

$$\begin{aligned} &+\delta(1-z)C\_{F}\left[C\_{F}\left(\frac{3}{8}-\frac{\pi^{2}}{2}+6\zeta\_{3}\right)+C\_{A}\left(\frac{1}{4}-3\zeta\_{3}\right)+\beta\_{0}\left(\frac{1}{8}+\frac{\pi^{2}}{6}\right)\right] \\ &-C\_{F}^{2}\left\{\frac{1+z^{2}}{1-z}\left[2\ln(1-z)+\frac{3}{2}\right]\ln z + \frac{1+z}{2}\ln^{2}z + \frac{3+7z}{2}\ln z + 5(1-z)\right\} \\ &+C\_{A}C\_{F}\left[\frac{1}{2}\frac{1+z^{2}}{1-z}\ln^{2}z + (1+z)\ln z + 3(1-z)\right] \\ &+\beta\_{0}\left[\frac{1}{2}\frac{1+z^{2}}{1-z}\ln z + 1-z\right], \\ \mathcal{P}^{(2)}\_{q\bar{q}V}(z) = \left(2C\_{F}-C\_{A}\right)C\_{F}\left\{\frac{1+z^{2}}{1+z}\left[S\_{2}(z) + \frac{1}{2}\ln^{2}z\right] + (1+z)\ln z + 2(1-z)\right\}, \\ \mathcal{P}^{(2)}\_{q\bar{q}S}(z) = T\_{R}C\_{F}\left[-(1+z)\ln^{2}z + \left(1+5z+\frac{8}{3}z^{2}\right)\ln z + \frac{20}{9z}-2+6z-\frac{56}{9}z^{2}\right], \\ \end{aligned} \tag{6.45}$$

and

$$\begin{split} \mathcal{P}\_{qg}^{(2)}(z) &= C\_{F}T\_{R} \Big[ \left( z^{2} + (1-z)^{2} \right) \left( \ln^{2} \frac{1-z}{z} - 2 \ln \frac{1-z}{z} - \frac{\pi^{2}}{3} + 5 \right) + 2 \ln(1-z) \\ & \quad - \frac{1-2z}{2} \ln^{2} z - \frac{1-4z}{2} \ln z + 2 - \frac{9}{2}z \Big] \\ &+ C\_{A}T\_{R} \Big\{ \left( z^{2} + (1-z)^{2} \right) \left[ -\ln^{2}(1-z) + 2 \ln(1-z) + \frac{22}{3} \ln z - \frac{109}{9} + \frac{\pi^{2}}{6} \right] \\ &+ \left( z^{2} + (1+z)^{2} \right) S\_{2}(z) - 2 \ln(1-z) - (1+2z) \ln^{2} z \\ &+ \frac{68z-19}{3} \ln z + \frac{20}{9z} + \frac{91}{9} + \frac{7}{9}z \Big]. \end{split}$$

The solution of the PDF evolution equations given earlier, in the presence of these splitting functions, is obtained numerically. Some examples of the effects of DGLAP evolution on the PDFs will be given later in this chapter.

## 6.2 Fitting parton distribution functions

#### 6.2.1 Processes involved in global analysis fits

#### <IP_ADDRESS> Interlude: Deep-inelastic scattering

The parton distribution functions were introduced in Chapter 2 with reference to deep-inelastic scattering (DIS) of leptons and protons. These processes play an important role in the extraction of PDFs since, within the parton model, the PDFs are related in a fairly straightforward manner to cross-sections that can be measured experimentally. The kinematics of a DIS process can be described by the following variables: Q2 , the square of the four-momentum transferred to the proton in the exchange; x, the fraction of the proton's momentum carried by the struck quark, y, the fraction of the incident lepton's energy, in the proton's rest frame, that is lost in the

#### 412 Parton Distribution Functions

collision. In terms of these variables the doubly-differential DIS cross-section can be written as,

$$\frac{\mathrm{d}^2 \sigma}{\mathrm{d}x \mathrm{d}y} = \frac{2\pi \alpha^2}{xyQ^4} \left[ \left( 1 + (1 - y)^2 \right) F\_2 - \left( 1 - (1 - y)^2 \right) x F\_3 - y^2 F\_L \right] \,. \tag{6.47}$$

This introduces a parameterization of the cross-section in terms of the structure functions F2, F3, and FL. This formula is correct for neutral current (NC) processes where the exchanged particle is a photon or a Z-boson. For charged current (CC) interactions the particle exchanged is a W-boson and, for a beam of the appropriate helicity lepton to facilate the weak interaction, Eq. (6.47) must be multiplied by an additional coupling and propagator factor,

$$2\left(\frac{G\_F m\_W^2}{4\pi\alpha} \frac{Q^2}{Q^2 + m\_W^2}\right)^2. \tag{6.48}$$

In the quark-parton model the relationship between the structure functions and the PDFs is:

$$F\_2^{\rm NC} = x \sum\_q C\_q \left[ f\_q(x, Q^2) + f\_{\overline{q}}(x, Q^2) \right] \,, \tag{6.49}$$

$$F\_3^{\rm NC} = \sum\_q C\_q' \left[ f\_q(x, Q^2) - f\_{\bar{q}}(x, Q^2) \right] \,, \tag{6.50}$$

$$F\_2^{\text{CC}(-)} = 2x \left[ f\_u(x, Q^2) + f\_{\bar{d}}(x, Q^2) + f\_{\bar{s}}(x, Q^2) + f\_c(x, Q^2) + \dots \right], \quad \text{(6.51)}$$

$$F\_3^{\text{CC}(-)} = 2\left[f\_u(x, Q^2) - f\_{\bar{d}}(x, Q^2) - f\_{\bar{s}}(x, Q^2) + f\_c(x, Q^2) + \dots\right].\tag{6.52}$$

Here, Cq and C 0 q represent the combinations of couplings and propagators necessary to account for the separate photon and Z-boson contributions to the neutral current process. In particular, F3 would be zero for the case of pure photon exchange. The charged-current process corresponds to the case of an incident electron, so that a W− is exchanged, with the corresponding result for W+ obtained by interchanging d ↔ u and s ↔ c. Since the contribution of sea quarks and anti-quarks is equal, measurements of F NC 3 are particularly useful probes of the valence quark distributions.

#### <IP_ADDRESS> Interlude: End

Measurements of deep-inelastic scattering (DIS) structure functions (F2, F3), or of the related cross-sections, in lepton-hadron scattering and of lepton pair production crosssections in hadron-hadron collisions provide the main source of information on quark distributions fq/p(x, µ2 F ) inside hadrons. For simplicity, for the rest of this chapter, the scale µF will be replaced by the scale Q, representing both the renormalization and factorization scales.

At leading-order, the gluon distribution function fg/p(x, Q2 ) enters directly in hadron-hadron scattering processes with jet final states. Modern global parton distribution fits are carried out to NLO and NNLO, which allows αS(Q2 ), fq/p(x, Q2 ) and fg/p(x, Q2 ) to all mix and contribute in the theoretical formulae for all processes.

Nevertheless, the broad picture described above still holds to some degree in global PDF analyses. Again, Q here refers to a scale representing the hardness of the interaction, as for example the jet transverse momentum in a calculation of the inclusive jet cross-section.

A NLO (NNLO) global PDF fit requires thousands of iterations and thus thousands of estimates of NLO (NNLO) matrix elements. The NLO (NNLO) matrix elements require too much time for evaluation to be used directly in global fits. Previously, a Kfactor (NLO/LO or NNLO/LO) was calculated for each data point used in the global fit, and the LO matrix element (which can be calculated very quickly) was changed in the global fit (multiplied by the K-factor). Currently, a routine such as fastNLO [686] or Applgrid [329] is often used for fast evaluation of the NLO matrix element with the new iterated PDF. Practically speaking, both provide the same order of accuracy. Even when fastNLO or Applgrid is used at NLO, a K-factor approach (NNLO/NLO) is still needed at NNLO, at least until the fastNLO/Applgrid technique can be adapted for NNLO calculations (in progress at the completion of this book).

The data from DIS, DY, and jet processes utilized in PDF fits cover a wide range in x and Q2 . HERA data [77, 115] are predominantly at low x, while the fixed target DIS [169, 170, 217, 900] and DY [769, 874] data are at higher x. Collider jet data at both the TEVATRON and LHC [9, 61, 86, 103, 116, 117, 119, 126, 173, 360] cover a broad range in x and Q2 by themselves and are particularly important in the determination of the high x gluon distribution. Jet data from the LHC have now been used in global PDF fits, and their importance will increase as high statistics data, and their detailed systematic error information, are published. In addition, jet production data from HERA have been used in the HERAPDF global PDF fits [78, 79, 389, 390].

As an example, the kinematic coverage of the data used in the NNPDF2.3 fit [192] is shown in Fig. 6.2.

There is a tradeoff between the size and the consistency of a data set used in a global PDF fit, in that a wider data set contains more information, but information coming from different experiments may be partially inconsistent. Most of the fixed target data have been taken on nuclear targets and suffer from uncertainties in the nuclear corrections that must be made [693]. This is unfortunate as it is the neutrino fixed target data that provide most of the quark flavour differentiation, for example between up, down, and strange quarks. As LHC collider data become more copious, it may be possible to reduce the reliance on fixed target nuclear data. For example, the rapidity distributions for W+, W−, and Z production at the LHC (as well as the TEVATRON) are proving to be very useful in constraining u and d valence and sea quarks, as described in Chapter 9.

There is considerable overlap, however, for the kinematic coverage among the datasets with the degree of overlap increasing with time as the full statistics of the HERA experiments have been published. Parton distributions determined at a given x and Q2 'feed-down' or evolve to lower x values at higher Q2 values, as discussed in Chapter 2. DGLAP-based NLO and NNLO pQCD should provide an accurate description of the data (and of the evolution of the parton distributions) over the entire kinematic range present in current global fits. At very low x and Q2 , DGLAP evolution is believed to be no longer applicable and a BFKL [189, 517, 708, 709] description

Fig. 6.2 The kinematical coverage in x and Q 2 for the data sets included in the NNPDF2.3 global PDF fit. Reprinted with permission from Ref. [192].

should be used. No clear evidence of BFKL physics is seen in the current range of data; thus all global analyses use conventional DGLAP evolution of PDFs.

There is a remarkable consistency between the data in the PDF fits and the perturbative QCD theory fit to them. The CT,MMHT, and NNPDF groups use over 3000 data points in their global PDF analyses and the χ 2/DOF for the fit of theory to data is on the order of unity, for both the NLO and NNLO analyses. For most of the data points, the statistical errors are smaller than the systematic errors, so a proper treatment of the systematic errors and their bin-to-bin correlations is important. All modern day experiments provide the needed correlated systematic error information. The H1 and ZEUS experiments have combined the data from the two experiments from Run 1 at HERA (and now Run 2) in such a way as to reduce both the systematic and statistical errors, providing errors of both types of the order of a percent or less over much of the HERA kinematics [77]. In the Run 1 combination,for example, 1402 data points are combined to form 742 cross-section measurements (including both neutral current and charged current cross-sections). The combined data sets, with their small statistical and systematic errors, form a very strong constraint for all modern global PDF fits. Thus, it can be hard for other data sets, for example from the LHC, to match the statistical and systematic errors of the HERA data. The manner of using the systematic errors in a global fit will be discussed later in Section 6.3.

The accuracy of the extrapolation to higher Q2 depends on the accuracy of the

original measurement, any uncertainty on αS(Q2 ) and the accuracy of the evolution code. Most global PDF analyses are carried out at NLO and NNLO. Both the NLO and the NNLO evolution codes have now been benchmarked against each other and found to be consistent [190, 191, 251, 262, 470, 572]. Most processes of interest have been calculated to NLO and there is the possibility, as discussed previously, of including data from these processes in global fits. Fewer processes have been calculated at NNLO [134, 296]. The processes that have been calculated include DIS, DY, diphoton [340], tt¯ production [427], and inclusive W,Z, and Higgs boson + jet production [269, 271, 272, 274]. Late in the writing of this book, the complete NNLO inclusive jet production cross-section has been completed (a monumental feat), but the results are not yet in the form to be easily used in global PDF fits [422]. Typically, jet production has been included in global PDF fits using NLO matrix elements. Threshold corrections [445, 676] can be used to make an approximate NNLO prediction, but the corrections are valid only over a limited phase space at the LHC, thus greatly reducing the size and power of the jet data in the global fits. Thus, any of the NNLO global PDF analyses discussed here are still approximate for this reason, but in practice the approximation should work reasonably well. The NNLO corrections for the inclusive jet cross-section have been found to be small (and relatively constant for the LHC phase space), if a scale equal to the transverse momentum of the jet is used [335].1 The CT14, MMHT2014, and NNPDF3.0 PDFs follow different philosophies regarding the use of LHC jet data in NNLO fits. CT14 makes no cuts on the LHC jet data, MMHT2014 doesn't include the LHC jet data, and NNPDF3.0 uses only the jet data for which threshold resummation provides a reasonable prediction. Current evolution programmes should be able to carry out the evolution using NLO and NNLO DGLAP to an accuracy of a few percent over the hadron collider kinematic range, except perhaps at very large and very small x.

The kinematics appropriate for the production of a state of mass M and rapidity y at the LHC is shown in Fig. 6.3 [320]. For example, to produce a state of mass 100 GeV and rapidity 2 requires partons of x values 0.05 and 0.001 at a Q2 value of 1 × 104 GeV2 . Compare this figure to the scatterplot of the x and Q2 range included in the recent NNPDF2.3 fit and it is clear that an extrapolation to higher Q2 (M2 ) is required for predictions for many of the LHC processes of interest. As more Standard Model processes are included in global PDF fits, the need for extrapolation will be reduced.

#### 6.2.2 Parameterizations and schemes

A global PDF analysis carried out at NLO or NNLO needs to be performed in a specific renormalization and factorization scheme. The evolution kernels are calculated in a specific scheme and to maintain consistency, any hard scattering cross section calculations used for the input processes or utilizing the resulting PDFs need to have been implemented in that same renormalization scheme. As we saw earlier in Chapter 2, one needs to specify a scheme or convention in subtracting the divergent terms from the PDFs. These divergent terms result from collinear gluon emission from the initial

1The NNLO cross-section is higher than the NLO one at smaller jet transverse momenta if a scale equal to that of the largest transverse momentum jet in the event is used.

2

Fig. 6.3 A plot showing the x and Q values needed for the colliding partons to produce a final state with mass M and rapidity y at the LHC (14 TeV).

state partons. The collinear emissions result in pole terms of the form 1/, where is the dimensional regularization parameter. Basically the scheme definition specifies how much of the finite corrections to subtract along with the divergent pieces. Almost universally, the MS scheme is used; using dimensional regularization, in this scheme the pole terms and accompanying log 4π and Euler constant terms are subtracted.2 PDFs are also available in the DIS scheme (where the full order αs corrections for F2

2Within the MS scheme, PDFs can also be defined for a fixed number of flavours, which then have a validity over the kinematic range for which that number (and only that number) of flavours can be present in the proton.

Fig. 6.4 CTEQ6.5 up and down quark distributions normalized to those of CTEQ6.1, showing the impact of the heavy quark mass corrections. Reprinted with permission from Ref. [876].

are absorbed into the quark PDFs).

Basically all modern PDFs now incorporate a treatment of heavy quark effects in their fits, either via the ACOT general-mass (GM) variable flavour number scheme [131] (supplemented by by a unified treatment of both kinematical and dynamical effects using the S-ACOT [694] and ACOT-χ [698, 875] concepts), used by CTEQ/CT, or by the Thorne-Roberts scheme [872, 873], used by both MSTW and HERAPDF, and the FONLL scheme, used by NNPDF [302, 532].

Incorporation of the full heavy-quark mass effects in the general-mass formalism suppresses the heavy flavour contributions to the DIS structure functions, especially at low x and Q2 . In order for the theoretical calculations in the global fits to agree with the data in these kinematic regions, the contributions of the light quark and anti-quark PDFs must increase accordingly. This has a noticeable impact, especially on predictions for W and Z cross-sections at the LHC.

Fig. 6.4 shows the impact of the heavy quark mass corrections on the up and down quark distributions for CTEQ6.5, at a Q value of 2 GeV [876]. The CTEQ6.5 up and down quark distributions are normalized to the corresponding ones from CTEQ6.1 (which does not have the heavy quark mass corrections). The shaded areas indicate the CTEQ6.1 PDF uncertainty. The dashed curves represent slightly different parameterizations for the CTEQ6.5 PDFs. The heavy quark mass corrections have a strong effect (larger than the PDF uncertainty for CTEQ6.1) at low x, in a region sensitive to W and Z production at the LHC.

The impact of general-mass variable flavour number schemes (GM-VFNS) lies mostly in the low x and Q2 regions. Aside from modifications to the fits to the HERA

#### 418 Parton Distribution Functions

data, and the commensurate change in the fitted PDFs, there is basically no modification for predictions at high Q2 at the LHC.

It is also possible to use only leading-order matrix element calculations in the global fits which results in leading-order parton distribution functions, which have been made available, for example, by the CTEQ [489, 818], MSTW/MMHT [614, 755] and NNPDF [191, 194] groups. For many hard matrix elements for processes used in the global analysis, there exist K factors significantly different from unity. Thus, one expects there to be noticeable differences between the LO and NLO parton distributions (and indeed this is often the case, especially at low x and high x).

Global analyses have traditionally used a generic form for the parameterization of both the quark and gluon distributions at some reference value Q0: 3

$$F(x, Q\_0) = x^{A\_1} (1 - x)^{A\_2} P(x; A\_3, A\_4 \dots). \tag{6.53}$$

The reference value Q0 is usually chosen in the range of 1–2 GeV. The parameter A1 is associated with small-x Regge behaviour while A2 is associated with large-x valence counting rules. We expect A1 to be approximately -1 for gluons and anti-quarks, and of the order of 1/2 for valence quarks, from the Regge arguments mentioned in Chapter 2. Counting rule arguments tell us that the A2 parameter should be related to 2ns − 1, where ns is the minimum number of spectator quarks. So, for valence quarks in a proton, there are two spectator quarks, and we expect A2 = 3. For a gluon, there are three spectator quarks, and A2 = 5; for anti-quarks in a proton, there are four spectator quarks, and thus A2 = 7. Such arguments are useful, for example in telling us that the gluon distribution should fall more rapidly with x than quark distributions, but it is not clear exactly at what value of Q that the arguments made above are valid.

The first two factors, in general, are not sufficient to completely describe either quark or gluon distributions. The term P(x; A3, ...) is a suitably chosen smooth function, depending on one or more parameters, that adds more flexibility to the PDF parameterization. P(x; A3, ...) is chosen so as to tend towards a constant for x approaching either 0 or 1, so that the limiting behaviour is given by the first two terms.

In general, both the number of free parameters and the functional form can have an influence on the global fit. A too-limited parameterization not only can lead to a worse description of the data, but also to PDFs in different kinematic regions being tied together not by the physics, but by the limitations of the parameterization. Note that the parameterization forms shown here imply that PDFs are positive-definite. As they are not physical objects by themselves, it is possible for them to be negative, especially at low Q2 . Some PDF groups (such as CT) use a positive-definite form for the parameterization; others do not. For example, the MSTW2008 gluon distribution is negative for x < 0.0001, Q2 = 2GeV 2 . Evolution quickly brings the gluon into positive territory.

The CT14 fit uses 28 free parameters (many of the PDF parameters are either fixed at reasonable values, or are constrained by sum rules). There are a total of 8 free

3Recently, there has been a trend towards the use of more sophisticated forms of parameterization in global fits, but the physics arguments listed here are still valid. For example, in the CT14 global fit [489], P(x) is defined by a fourth-order polynomial in √ x; the polynomial is then re-expressed in terms of Bernstein polynomials in order to reduce correlations among the coefficients.

Fig. 6.5 The CT14 NNLO parton distribution functions evaluated at a Q 2 value of 2 GeV2 , with a linear (left) and logarithmic (right) scale.

parameters for the valence quarks, 5 for the gluon and 15 for the sea quarks.

The MMHT2014 fit uses 20 free parameters, while the NNPDF fits effectively has 259 free parameters. The NNPDF approach attempts to minimize the parameterization bias by exploring global fits using a large number of free parameters in a Monte Carlo approach. The general form for NNPDF can be written as fi(x, Qo) = ci(x)NNi(x), where NNi(x) is a neural network, and ci(x) is a "pre-processing function".

In the past, PDFs were often made available to the world in a form where the x and Q2 dependence was parameterized. Now, almost universally, the PDFs for a given x and Q2 range can be interpolated from a grid that is provided by the PDF groups, or the grid can be generated given the starting parameters for the PDFs (see the discussion on LHAPDF in Section 6.6). All techniques should provide an accuracy on the output PDF distributions on the order of a few percent or better.

The parton distributions from the CT14 NNLO PDFs are plotted in Fig. 6.5 at a Q2 value of 2 GeV 2 (near the starting point of evolution) and in Fig. 6.6 at a Q2 value of 10000 GeV 2 (more typical of LHC processes). At the lower Q2 value, the up quark and down quark distribution peak near x values of 1/3, the remnant of the spike in the primitive model described in Chapter 2. The charm PDF is very suppressed as it is produced entirely by evolution, and the Q2 value is near the starting scale for the evolution. There is no bottom quark distribution since it is below threshold. At higher Q2 values, the up and down quark peaks become shoulders due to the effects of evolution. At high Q2 , the gluon distribution is dominant at x values of less than 0.1 with the valence quark distributions dominant at higher x. One of the major influences of the HERA data has been to steepen the gluon distribution at low x.

The CT14 up quark, up-bar quark, b-quark, and gluon distributions are shown as a function of Q2 for x values of 0.001, 0.01, 0.1 and 0.3 in Figs. 6.7 and 6.8. At low x, the PDFs increase with Q2 , while at higher x, the PDFs decrease with Q2 . Both effects are due to DGLAP evolution, as discussed previously. An x value of approximately

Fig. 6.6 The CT14 NNLO parton distribution functions evaluated at a Q 2 of 10000 GeV2 .

Fig. 6.7 The CT14 NNLO up quark, up-bar quark, c quark, and gluon parton distribution functions evaluated as a function of Q 2 at x values of 0.001 (left) and 0.01 (right).

0.1 is the pivot − point for gluon evolution; at this x value, the gluon distribution changes little as the value of Q2 increases. It can also be seen that both the charm (and bottom) quark distributions are generated perturbatively from gluon splitting, and thus the distributions are zero below threshold for heavy quark pair production and rise rapidly thereafter with increasing Q2 . It is also possible for there to be intrinsic charm, where charm quarks are present below threshold. However, no strong evidence has been observed to date for the existence of intrinsic charm.

Fig. 6.8 The CT14 NNLO up quark, up-bar quark, c quark and gluon parton distribution functions evaluated as a function of Q 2 at an x value of 0.1 (left) and 0.3 (right).

Fig. 6.9 The momentum fractions carried by the CT14 NNLO quark and gluon distributions, as a function of Q. The gluon distribution in the right figure is shown without (solid) and with (dotted) the presence of a top quark PDF.

The average proton momentum carried by each parton species is shown in Fig. 6.9. As Q2 increases, the momentum carried by up and down valence quarks decreases, while the momentum carried by the gluon and by sea quarks increases. For typical LHC hard-scattering scales, the gluon carries slightly less than 50% of the parent proton's momentum. Note that a 5-flavour scheme is most commonly used, i.e. the charm and bottom (but not top) can appear as sea quarks once the Q value is sufficient to pair produce them from gluon splitting. It is also possible to allow top quarks in the sea in a 6-flavour scheme; even at the highest Q values, only about 1% of the proton's momentum is carried by top quarks. The momentum added to the top quark distribution comes primarily from the gluon distribution.

The photon is also a parton constituent of the proton, just as a quark or gluon is, and can be produced from QED radiation from quark lines [754]. This source of

#### 422 Parton Distribution Functions

photons is known as the inelastic component. There is also a (well-known) elastic component, which cannot be ignored, resulting from coherent electromagnetic radiation from the proton as a whole, leaving the proton intact [292, 587]. Both components contribute to photon-induced processes at the LHC. The inelastic component evolves with Q2 , while the elastic component is relatively constant (varying with the running of the QED coupling). To include the photon PDF, the QCD evolution of the partons in the proton now has to be expanded to a QCD+QED evolution, where the QED aspect involves the electromagnetic coupling α(Q2 ) instead of αs(Q2 ) and the corresponding splitting function is Abelian rather than non-Abelian. At the LHC, especially for the 13-14 TeV running, processes involving photons in the initial state will become increasingly important, c.f. γγ → WW, or photon-initiated production of W H. There is little hadronic data to directly constrain the photon PDF, but it has been known that less than 1% of the proton's momentum is carried by photons. The first attempt (MRST2004QED) to model the inelastic component just considered photon emission from quark lines, using the known quark distributions, and using either the current quark mass (few MeV) or the constituent quark mass (few hundred MeV) as a cutoff [754] (see also Refs. [587, 772]). The disparity in the cutoffs leads to a wide range in the possible size of the photon PDF. Other attempts to determine the photon distribution used Drell-Yan data from the LHC [193, 194, 275](NNPDF2.3qed followed by NNPDF3.0qed) to fit for the photon PDF.4

Another approach [80, 837] (CT14qed,CT14qed inc) used data from the scattering process ep → eγX measured by the ZEUS collaboration [391], leading to an upper constraint on the total inelastic photon PDF momentum, at a scale of 1.3 GeV, of approximately 0.14% (and a lower constraint of 0). This is to be compared to the total elastic photon PDF momentum fraction of 0.15%. The photon PDFs for the inelastic component only, assuming the maximum intrinsic momentum fraction of 0.14%, and the total photon PDF, including the elastic component as well, are shown in Fig. 6.10 for Q=1.3 GeV (left) and Q=85 GeV (right). The dominance of the inelastic component at high Q can be observed.

Recently, the photon PDF has been determined to a high precision (1-2%) using electron-proton scattering data, considering an equivalence between a template crosssection calculated either using proton structure functions, or a photon PDF [746]. The resulting PDF (LUXqed) is in good agreement with CT14qed inc and is at the lower edge of the NNPDF2.3qed photon PDF uncertainty band at high x, as shown in Fig. 6.11.

The charm quark distribution also has a dynamic component, generated through gluon splitting into a cc¯ pair.5 The photon/charm ratio increases with increasing Q and increasing x value. One reason for the variation in Q is that while αs decreases with Q, α remains approximately constant (actually, it rises slightly). At low x, the photon/charm ratio is of the order of 5-10%, due to the difference in coupling constants

4NNPDF3.0qed improves on the NNPDF2.3qed photon PDF with a correct treatment of α(αsL) n terms in the evolution. These resulted in a large uncertainty for the photon PDF, as the errors in the reasonably precise high-mass Drell-Yan data are still large compared to the expected contributions from photon-initiated processes.

5There may also be an intrinsic charm component at low Q, but the evidence is not convincing.

Fig. 6.10 The photon PDF from the inelastic component only, with a photon momentum fraction of 0.14%, and the total photon PDF, including the elastic component, at a Q value of 1.3 GeV (left) and 85 GeV (right).

Fig. 6.11 The ratio of various photon PDF sets to that of the LUXqed photon PDF set, all evaluated at a scale of 100 GeV. Note that the vertical axis scales are different for the sub-plots. Reprinted with permission from Ref. [746].

(αs vs. α) and the larger gluon than quark (primarily up quark) distribution at low x. (Remember that the photon does have a substantial elastic component at small x which is included in this ratio.) At high x, the dominance of the valence up quark over the gluon results in the photon distribution becoming larger than the charm distribution.

#### 424 Parton Distribution Functions

## 6.3 PDF uncertainties

In addition to having the best estimates for the values of the PDFs in a given kinematic range, it is also important to understand the allowed range of variation of the PDFs, i.e. their uncertainties. A prior conventional method of estimating parton distribution uncertainties had been to compare different published parton distributions. Although in some cases, this procedure does provide useful information, in general this is unreliable since most published sets of parton distributions adopt similar assumptions and the differences between the sets do not fully explore the uncertainties that actually exist. In addition, some PDF fits may use only a limited set of data in their fits, and they may result in PDFs that (a) have larger intrinsic uncertainties and (b) differ in their central values from PDF sets which use a more global set of input data. Some comparisons of PDFs and PDF predictions at the LHC will be given later in this chapter.

The sum of the quark distributions (Σfq/p(x, Q2 ) + fg/p(x, Q2 )) is, in general, well determined over a wide range of x and Q2 . As stated earlier, the quark distributions are predominantly determined by the DIS and DY data sets which have large statistics, and systematic errors in the few percent range (±3% for 10−4 < x < 0.75). Thus the sum of the quark distributions is basically known to a similar accuracy. The individual quark flavours, though, may have a greater uncertainty than the sum. This can be important, for example, in predicting distributions that depend on specific quark flavours, like the W lepton asymmetry distribution and the W and Z rapidity distributions.

The largest uncertainty of any parton distribution, however, is that on the gluon distribution. The gluon distribution can be determined indirectly at low x by measuring the scaling violations in the quark distributions, but a direct measurement is necessary at moderate to high x. About 40-50% of the momentum of the proton is carried by gluons, and most of that momentum is at relatively small x (16% of the momentum of the proton, for example, is carried by gluons in the x range from 0.01 to 0.1.) The best direct information on the gluon distribution at moderate to high x comes from jet production at the TEVATRON and the LHC, although new processes (photon production, top pair production, etc.) will increasingly contribute.

There has been a great deal of activity on the subject of PDF uncertainties. Two techniques in particular, the Lagrange Multiplier and Hessian techniques, have been used by CTEQ/CT, MSTW/MMHT, and HERAPDF to estimate PDF uncertainties [413, 753, 817, 819, 868]. The Lagrange Multiplier technique is useful for probing the PDF uncertainty of a given process, such as the Higgs boson cross-section, while the Hessian technique provides a more general framework for estimating the PDF uncertainty for any cross-section. In addition, the Hessian technique results in tools more accessible to the general user, such as error PDFs. Both techniques are described in the following. The Monte Carlo technique used by the NNPDF group for PDF uncertainty estimation will be described later in this section.

#### 6.3.1 The Hessian method

The Hessian method for the determination of a central PDF and/or determination of PDF uncertainties involves minimizing a suitable log-likelihood function. The χ 2

Fig. 6.12 A comparison of the unshifted and shifted HERA1 combined neutral current data (left) and the comparison of the NLO CT10 predictions to the shifted data (right). Reprinted with permission from Ref. [713].

function may contain the full set of correlated errors, or only a partial set. The correlated systematic errors may be accounted for using a covariance matrix, or as a shift to the data, adopting a χ 2 penalty proportional to the size of the shift divided by the systematic error. The two methods should be equivalent. Below we discuss the shift method.

In the following description, CT10 NLO PDFs [713] are used (similar considerations apply at NNLO), along with the HERA Run I combined (H1+ZEUS) neutral current (e +p) cross sections [77], to discuss the use of the Hessian formalism. A comparison of the HERA data and the NLO predictions using the CT10 PDFs is shown in Fig. 6.12. On the left, the data are presented in unshifted form, and on the right, optimal systematic error shifts have been applied in the manner detailed below. There is good agreement between the combined HERA Run I data and the CT10 NLO predictions with a global χ 2 of about 680 for the 579 data points, which is typical for the global fit PDFs.

The HERA Run I combined data have Nλ = 114 independent sources of experimental systematic uncertainty, with parameters λα that should obey a standard Gaussian or normal distribution. The contribution of the HERA dataset to the χ 2 can be written as

$$\chi^2(\{a\}, \{\lambda\}) = \sum\_{k=1}^{N} \frac{1}{s\_k^2} \left( D\_k - T\_k(\{a\}) - \sum\_{\alpha=1}^{N\_\lambda} \lambda\_\alpha \beta\_{k\alpha} \right)^2 + \sum\_{\alpha=1}^{N\_\lambda} \lambda\_\alpha^2,\tag{6.54}$$

where N is the total number of points and Tk(a) is the theory value for the kth data

Fig. 6.13 Distribution of systematic parameters λα of the combined HERA Run I data set [871] in the CT10 best fit (CT10.00). Reprinted with permission from Ref. [713].

point, dependent on the PDF parameters a. Furthermore, sk is the total uncorrelated error on the measurement Dk that is obtained by summing the statistical and uncorrelated systematic errors on Dk in quadrature,

$$s\_k = \sqrt{s\_{k, \text{stat}}^2 + s\_{k, \text{unccorr sys}}^2} \tag{6.55}$$

The χ 2 function is minimized with respect to the size of the systematic error shifts λα using the algebraic procedure described. It is also possible to add a penalty term to the χ 2 function that prevents relatively unconstrained PDF parameters from reaching values that might lead to unphysical predictions in regions where experimental data are sparse (e.g. very low x).

As expected, a better agreement of data with theory is observed when the systematic error shifts are allowed. It is important to check that the systematic error parameters λα(a) contribute to the χ 2 an amount on the order of the total number of systematic errors (114) and (b) that the sizes of the parameters follow a Gaussian distribution. For the case of CT10 and the HERA Run I data, the systematic error contribution to the total χ 2 is 65, or somewhat better than expected, and their distribution is approximately Gaussian-distributed, as shown in Fig. 6.13.

All systematic errors are not equally important, though, and it is also crucial to verify that no 'major' systematic error needs to be shifted by several sigma. Given the precision of the HERA Run I data, the size of the systematic error shifts required are relatively small. This need not be the case, for example, for the case of inclusive jet production, either at the TEVATRON or the LHC.

The Hessian method results in the production of a central (best fit) PDF, and a set of error PDFs. In this method, a large matrix (26 × 26 for CT10 and 28 × 28 for CT14), with dimension equal to the number of free parameters in the fit, has to

#### PDF uncertainties 427

Fig. 6.14 A schematic representation of the transformation from the PDF parameter basis to the orthonormal eigenvector basis. Reprinted with permission from Ref. [818].

be diagonalized.6 The result is 26 (28) orthonormal eigenvector directions for CT10 (CT14) which provide the basis for the determination of the PDF error for any crosssection.

This process is shown schematically in Fig. 6.14. The eigenvectors are now admixtures of the PDF parameters left free in the global fit. There is a broad range for the eigenvalues, over a factor of one million. The eigenvalues are distributed roughly linearly as log i , where i is the eigenvalue for the i-th direction. The larger eigenvalues correspond to directions which are well-determined; for example, eigenvectors 1 and 2 are sensitive primarily to the valence quark distributions at moderate x, a region where they are well-constrained. The theoretical uncertainty on the determination of the W mass at the TEVATRON depends primarily on these 2 eigenvector directions, as W production at the TEVATRON proceeds primarily through collisions of valence quarks. The most significant eigenvector directions for determination of the W mass at the LHC correspond to larger eigenvector numbers, which are primarily determined by sea quark distributions. In most cases, the eigenvector can not be directly tied to the behaviour of a particular PDF in a specific kinematic region. There are exceptions, such as eigenvector 15 in the CTEQ6.1 fit, discussed in the following.

In the past, one of the most controversial aspect of PDF uncertainties has been the determination of the ∆χ 2 excursion from the central fit that is representative of a reasonable error. Nominally, a ∆χ 2 = T 2 (tolerance) would correspond to a 1 − σ(68%CL) error. PDF fits performed with a limited number of experiments may be able to maintain that criterion. For example, HERAPDF uses a χ 2 excursion of 1 for a 1σ error.7 For general global fits, such as from CT and MMHT, however, a χ 2

6As more data is included, more PDF parameters in the global fit can be set free, resulting in a larger number of eigenvectors.

7But the total error also includes other sources of uncertainty, for example from possible parameterization bias.

#### 428 Parton Distribution Functions

excursion of 1 (for a 1σ error) is too low of a value in a global PDF fit. These global fits use data sets arising from a number of different processes and different experiments; there is a non-negligible tension between some of the different data sets. In addition, the finite number of PDF parameters used in the parameterizations (parameterization bias) also leads to the need for a larger tolerance. Thus, a larger variation in ∆χ 2 is required for a 68% CL. For example, CT10 uses a tolerance T=10 for a 90% CL error, corresponding to T=6.1 for a 68% CL error,8 while MSTW uses a dynamical tolerance (varying from 1 to 6.5) for each eigenvector.

The uncertainties for all predictions should be linearly dependent on the tolerance parameter used; thus, it should be reasonable to scale the uncertainty for an observable from the 90% CL limit provided by the CT error PDFs to a one-sigma error by dividing by a factor of 1.645. Such a scaling will be a better approximation for observables more dependent on the lower number eigenvectors, where the χ 2 function is closer to a quadratic form.

Even though, the data sets and definitions of tolerance are different among the different PDF groups, we will see later in this chapter that the PDF uncertainties at the LHC are fairly similar. Note that relying on the errors determined from a single PDF group may be an underestimate of the true PDF uncertainty, as the central results among the PDF groups can in some cases differ by an amount similar to this one-sigma error. (See the discussion later in this chapter regarding benchmarking comparisons of predictions and uncertainties for the LHC.)

Each error PDF results from an excursion along the "+" and "−" directions for each eigenvector. Consider a variable X; its value using the central PDF for an error set (say CT14) is given by X0. X + i is the value of that variable using the PDF corresponding to the "+" direction for eigenvector i and X − i the value for the variable using the PDF corresponding to the "−" direction. The excursions are symmetric for the larger eigenvalues, but may be asymmetric for the more poorly determined directions. In order to calculate the PDF error for an observable, a Master Equation should be used:

$$\Delta X\_{max}^{+} = \sqrt{\sum\_{i=1}^{N} [\max(X\_i^{+} - X\_0, X\_i^{-} - X\_0, 0)]^2}$$

$$\Delta X\_{max}^{-} = \sqrt{\sum\_{i=1}^{N} [\max(X\_0 - X\_i^{+}, X\_0 - X\_i^{-}, 0)]^2}. \tag{6.56}$$

∆X+ adds in quadrature the PDF error contributions that lead to an increase in the observable X and ∆X− the PDF error contributions that lead to a decrease. The addition in quadrature is justified by the eigenvectors forming an orthonormal basis. The sum is over all N eigenvector directions, or 20 in the case of CTEQ6.1 and 26 (28) in the case of CT10 (CT14). Ordinarily, X + i − X0 will be positive and X − i − X0

8A penalty is applied in the CT10 approach, if a disproportionate fraction of the increase in χ 2 along a particular eigenvector direction is concentrated in one, or a few, experiments; this results in the tolerance being effectively smaller for those directions.

will be negative (or vice versa), and thus it is trivial as to which term is to be included in each quadratic sum. For the higher number eigenvectors, however, the "+" and "−" contributions may be in the same direction (see for example eigenvector 17 in Fig. 6.15). In this case, only the most positive term will be included in the calculation of ∆X+ and the most negative in the calculation of ∆X−. Thus, there may be less than N terms for either the "+" or "−" directions. There are other versions of the Master Equation in current use but the version listed above is the "official" recommendation of the authors.

There are two things that can happen when new PDFs (eigenvector directions) are added: a new direction in parameter space can be opened to which some cross-sections will be sensitive to (an example of this is eigenvector 15 in the CTEQ6.1 error PDF set, which is sensitive to the high x gluon behaviour and thus influences the high pT jet cross-section at the TEVATRON and LHC). This particular eigenvector direction happens to be dominated by a parameter which affects mostly the large x behaviour of the gluon distribution.

In this case, a smaller parameter space is an underestimate of the true PDF error since it did not sample a direction important for some physics. In the second case, adding new eigenvectors does not appreciably open new parameter space and the new parameters should not contribute much PDF error to most physics processes (although the error may be redistributed somewhat among the new and old eigenvectors).

In Fig. 6.15, the PDF errors are shown in the "+" and "−" directions for the 20 CTEQ eigenvector directions for predictions for inclusive jet production at the TEVATRON from the CTEQ6.1 PDFs. The excursions are symmetric for the first 10 eigenvectors but can be asymmetric for the last 10, as they correspond to less welldetermined directions.

Either X0 and X ± i can be calculated separately in a matrix element/Monte Carlo program (requiring the program to be run 2N + 1 times) or X0 can be calculated with the program and at the same time the ratio of the PDF luminosities (the product of the two PDFs at the x values used in the generation of the event) for eigenvector i (±) to that of the central fit can be calculated and stored. This results in an effective sample with 2N +1 weights, but identical kinematics, requiring a substantially reduced amount of time to generate. PDF re-weighting will be discussed later in this chapter.

As an example of PDF uncertainties using the Hessian method, the CT10 and MSTW2008 NLO uncertainties for the up quark and gluon distributions are shown in Figs. 6.16 and 6.17. While the CT10 and MSTW2008 PDF distributions and uncertainties are reasonably close to each other, some differences are evident, especially at low and high x.