Next-to-next-to-leading order event generation for top-quark pair production

The production of top-quark pairs in hadronic collisions is among the most important reactions in modern particle physics phenomenology and constitutes an instrumental avenue to study the properties of the heaviest quark observed in nature. The analysis of this process at the Large Hadron Collider relies heavily on Monte Carlo simulations of the final state events, whose accuracy is challenged by the outstanding precision of experimental measurements. In this letter we present the first matched computation of top-quark pair production at next-to-next-to-leading order in QCD with all-order radiative corrections as implemented via parton-shower simulations. Besides its intrinsic relevance for LHC phenomenology, this work also establishes an important step towards the simulation of other hadronic processes with colour charges in the final state.

The production of top-quark pairs in hadronic collisions is among the most important reactions in modern particle physics phenomenology and constitutes an instrumental avenue to study the properties of the heaviest quark observed in nature. The analysis of this process at the Large Hadron Collider relies heavily on Monte Carlo simulations of the final state events, whose accuracy is challenged by the outstanding precision of experimental measurements. In this letter we present the first matched computation of top-quark pair production at next-to-next-to-leading order in QCD with all-order radiative corrections as implemented via parton-shower simulations. Besides its intrinsic relevance for LHC phenomenology, this work also establishes an important step towards the simulation of other hadronic processes with colour charges in the final state.
Top quarks are the heaviest elementary particles observed in nature and play a unique role in the Standard Model (SM) of particle physics. The large Yukawa coupling of the top quark to the Higgs boson establishes a special avenue in the exploration of the Higgs sector of the SM [1] and of new physics signals [2]. Moreover, the value of the top-quark mass enters precision Electro-Weak (EW) tests [3] and theoretical considerations on the stability of our universe [4].
At the Large Hadron Collider (LHC), top quarks are predominantly produced by strong interactions in association with their own anti-particle (tt). The large value of the top-quark mass is such that its production dynamics is safely inside the region of validity of QCD perturbation theory. Remarkable theoretical advancements in the past years have led to very accurate predictions for this process. Specifically, fixed-order computations that rely on a power expansion in the strong coupling constant α s are known up to next-to-next-to-leading order (NNLO) [5][6][7][8][9][10][11][12] (also including the top-quark decays [13]), and EW corrections have been computed up to next-to-leading order (NLO) [14][15][16][17][18]. In specific kinematic regimes, a reliable perturbative description requires the all-order resummation of large radiative corrections [19][20][21][22][23][24][25]. Some of the above calculations have been consistently combined to obtain the state-of-the-art predictions at the LHC [26]. The striking accuracy of experimental measurements of the top-quark mass requires pushing theoretical calculations to the edge of what can be achieved with perturbative methods (for recent reviews see [27,28]), and motivates new studies of non-perturbative aspects of topquark physics (see e.g. [29][30][31][32]).
The large number of top-quark pairs produced at the LHC has allowed for very satisfactory tests of the theory, both for the inclusive production cross section and for multi-differential distributions [33][34][35][36][37][38][39][40][41]. These tests have paved the way to the exploitation of top cross-section measurements for the extraction of fundamental parameters of the Standard Model, such as α s and Parton Den-sity Functions (PDFs), the top mass itself and the top Yukawa couplings (see e.g. [42][43][44][45][46][47][48]).
Experimental analyses involving tt production heavily rely upon its fully exclusive simulation. This is important not only for the study of the production dynamics itself, but, due to the complex final states that involve a combination of leptons, jets, b hadrons and missing energy, also for several SM processes and new physics searches for which the tt process acts as background. The needed simulations rely on event generators, which combine a prediction for the high-energy scattering, that produces the tt pair, initial-and final-state QCD radiation at all perturbative orders via parton shower (PS) algorithms and hadronization models (for a review see [49]). These event generators have been subject of considerable research, dealing with the matching of NLO QCD calculations to PS and a consistent description of the top-quark resonance [50][51][52][53][54][55]. Current research for the improvement of event generators focuses upon the development of more accurate PS [56][57][58][59][60][61][62], as well as a framework to combine NNLO computations with PS into a consistent event generator (NNLO+PS in the following).
Different frameworks for NNLO+PS computations have been developed in recent years in the context of colour-singlet production [63][64][65][66][67]. However, nearly a decade after these developments, a NNLO+PS method to deal with hadron-collider processes with colour charges in the final state (e.g. tt), which are considerably more complex, is still missing.
In this letter, we show how a generator of the same type as the ones developed for colour singlet production processes in Refs. [66,67], dubbed there MiNNLO PS , can be constructed for top-quark pair production. Our work constitutes the first computation of this type for reactions with coloured particles in the final state.
The MiNNLO PS procedure involves three steps. The first one (referred to as Step I in the following) corresponds to the generation of a tt pair plus one light parton (i.e. the underlying Born configuration) according to the POWHEG method [68][69][70][71], carried out at the NLO level, inclusive over the radiation of a second light parton.
The second step (Step II) characterizes the MiNNLO PS procedure, and it concerns the limit in which the light partons in the above calculation become unresolved (i.e. the underlying Born degenerates into a tt configuration without light jets). In this limit the calculation must be supplemented with an appropriate Sudakov form factor and higher-order terms, so as to guarantee that the simulation remains finite as well as NNLO accurate for inclusive tt production. Most of the novelties in this letter have to do with Step II and will be illustrated below.
In the third step (Step III), the kinematics of the second radiated parton, accounted for inclusively in Step I, is generated according to the standard POWHEG method, which guarantees that the NLO accuracy of tt+jet cross section is preserved. From this point on, subsequent radiation is included by the parton shower, with the constraint of having a transverse momentum softer than that of the last POWHEG real emission.
The starting point to achieve NNLO accuracy in Step II is the well known factorization theorem for tt pair production at small transverse momentum p T ≡ | p T | differential in the phase space of the tt pair dΦ tt ≡ dx 1 dx 2 [dΦ 2 ]. Herex 1,2 = m tt / √ s e ±y tt , with y tt being the rapidity of the tt system, m tt its invariant mass, [dΦ 2 ] denotes the Lorentz-invariant two-body phase space and √ s is the collider centre-of-mass energy. It reads [19][20][21][22] dσ S c is the Sudakov radiator which also enters the description of the production of a colour singlet system at small transverse momentum The first sum in Eq. (1) runs over all possible flavour configurations of the incoming partons p 1 of flavour c and p 2 of flavourc. The collinear coefficient functions describe the structure of constant terms related to the emission of collinear radiation, and the parton densities are denoted by f i and are evaluated at b 0 /b. The operation ⊗ denotes the standard convolution over the momentum fraction z carried by initial state radiation. The factor Tr( has different expressions for the qq and gg channels and has here a symbolic meaning. In particular, it has a rich Lorentz structure that we omit for simplicity in Eq. (1), which is a source of azimuthal correlations in the collinear limit [21,72]. All quantities in bold face denote operators in colour space, and the trace Tr(H c ∆) in Eq. (1) runs over the colour indices. This term can be expressed conveniently in the colour space formalism of Ref. [73], where the infrared-subtracted amplitude |M cc for the production of the tt system is a vector in colour space [12,21]. It reads where |M cc . The hard function H c = H c (Φ tt ; α s (m tt )) is obtained from the subtracted amplitudes and the ambiguity in its definition corresponds to using a specific resummation scheme [74]. We adopt here the definition of Ref. [21]. The operator ∆ encodes the structure of the quantum interference due to the exchange of soft radiation at large angle between the initial and final state, and within the final state. It is given by The symbol P denotes the path ordering (with increasing scales from left to right) of the exponential matrix with respect to the integration variable q 2 . Γ t is the anomalous dimension accounting for the effect of real soft radiation at large angles, and D = D(Φ tt , b; α s (b 0 /b)) encodes the azimuthal dependence of the corresponding constant terms, and is defined such that [D] φ = 1, where [· · · ] φ denotes the average over the azimuthal angle φ of p T . All of the above quantities admit a perturbative expansion in a series in α s (µ)/(2π), where µ is the scale indicated explicitly in the argument of each function. We generically denote these expansions as with {x} being any other set of arguments of F ≡ {A, B, C ij , H c , D, Γ t , |M cc }. We do not indicate explicitly the scale of the amplitude |M cc . In this case, the expansion (5) is in powers of α s (m tt ), and each of its perturbative coefficients |M includes an extra single power of α s (µ R ∼ m tt . The above coefficients F (i) up to two loops are given in Refs. [12,20,21,72,[75][76][77][78][79][80][81][82][83][84][85][86]. The two-loop coefficient D (2) is irrelevant at NNLO for observables averaged over φ, since [D] φ = 1, and therefore we do not include it here.
Expanding the second order term in the exponential of Eq. (4), one can write where we neglected N 3 LL corrections which do not contribute at NNLO in Eq. (1). In the following we will denote by V NLL the right hand side of Eq. (6) with Γ (2) t → 0. This enters the description of the p T spectrum at NLL accuracy.
To make contact with the procedure described in Ref. [66] we need to simplify further the structure of the term H c ∆, which encodes the difference with the colour singlet case. Our goal is to obtain a closed formula in p T space that retains NNLO accuracy. To this order, we observe that we can take the Γ (2) t term in Eq. (6) out of the path ordering symbol. We then perform a rotation in colour space to diagonalize Γ (1) t and evaluate the exponential matrix in Eq. (6). Eq. (1) can be reorganized using where the trace is to be interpreted as in Eq.
(3). The Sudakov radiatorŜ c is obtained from S c via the replacement The remainder term E(Φ tt , b) in Eq. (7) contributes at order α 2 s ln(m tt b), but it is irrelevant for our computation since it vanishes upon azimuthal integration (i.e. [E] φ = 0). For this reason, we will ignore it in the following. We then obtain This expression has almost the structure needed in order to carry out the same procedure used in the colour singlet case [66,67], except for the H c function which needs to be evaluated at the scale b 0 /b rather than m tt . To the relevant accuracy we can perform this scale change provided B (2) inŜ c is also modified as follows [63,66] (see e.g. Eq. (4.25) of Ref. [66]) In the colour basis in which Γ is a linear combination of complex exponential terms, each of which has the same factorized structure used as a starting point in the appendix of Ref. [66].
Using this observation, we finally integrate over b by expanding the integrand about b 0 /b ∼ p T . Noticing that the matrix element M To obtain Eq. (11), we introduced the quantitiesS c ,H c andC ij , which are obtained by applying the transformations given in Eq. (4.24) of Ref. [66] toŜ c , H c and C ij . The latter quantities are now evaluated at the scale p T . The azimuthally averaged term [· · · ] φ in Eq. (11) is taken from the NNLO computation of the tt cross section of Refs. [11,12,86] (see also Ref. [87] for more details). We also included the remainder R f = R f (p T ), which denotes the regular contribution to the p T distribution through O(α 4 s ), such that p T R f (p T ) vanishes in the p T → 0 limit. The integral of Eq. (11) over p T provides a NNLO accurate description of tt production differential in Φ tt .
In order to build a Monte Carlo algorithm for the gen-eration of events with NNLO accuracy, we have to modify the formula for the underlying born cross section of Step I in such a way that it matches Eq. (11) maintaining its NNLO accuracy. The procedure is described in detail in Refs. [63,66,67,88] and for this reason we omit it here. For the practical implementation, we use the NLO+PS code for tt+jet of Ref. [89], and apply the MiNNLO PS procedure for heavy-quark pair production given in this letter. The PS simulation is obtained with Pythia 8 [90], without the modelling of non-perturbative effects, and under the assumption of stable top quarks. We stress that up to Eq. (10) we retained also NLL accuracy in the p T spectrum, while Eq. (11) is strictly LL accurate.
Higher logarithmic accuracy could in principle be maintained in Eq. (11). However, this higher accuracy would be spoiled by the PS used here, which is limited to LL. On the other hand, Eq. (11) also preserves the class of NLL corrections associated with the coefficient A (2) in the Sudakov, that are traditionally included in PS algorithms [91]. The formulation of a (N)NLO matching to PS that preserves logarithmic accuracy beyond LL is still an open problem.
In the phenomenological study presented below, we consider LHC collisions with a center of mass energy of 13 TeV. The top-quark pole mass is set to 173.3 GeV and we consider five massless quark flavours using the corresponding NNLO set of the NNPDF31 [92] parton densities with α s (m Z ) = 0.118. The renormalization scale for the two powers of the strong coupling constant entering the Born cross section is set to µ (0) R = K R m tt /2. In the rest of Eq. (11), we implement the renormalization (µ R = K R µ 0 ) and factorization (µ F = K F µ 0 ) scale dependence as described in Ref. [66], with the central scale µ 0 = m tt /2 e −L (hence replacing the scales set to p T in Eq. (11)), where we defined L = ln Q/p T and Q = m tt /2. The logarithm L is turned off in the hard region of the p T spectrum so that the total derivative in Eq. (11) smoothly vanishes for p T Q as in Refs. [74,[93][94][95][96]. The dependence of Eq. (11) on µ R , µ F and Q is of order O(α 5 s ). At small p T the scale of the strong coupling and the parton densities is smoothly frozen around Q 0 = 2 GeV following the procedure of Ref. [67] to avoid the Landau singularity. To estimate the scale uncertainties we vary K R and K F by a factor of 2 around their central value, while keeping 1/2 ≤ K R /K F ≤ 2. Results with a different central scale choice are reported in Ref. [97].
For comparison, we consider results from the fixedorder NNLO calculation of Ref. [11,12] obtained with the Matrix code [98] using µ 0 = m tt /2. Furthermore, we also show MiNLO results, obtained with the NLO+PS generator for tt plus zero and one jet, constructed by turning off the NNLO corrections in Eq. (11). The latter constitutes a new calculation as well. Table I shows the total cross section for top-quark pair production for MiNLO , NNLO and MiNNLO PS . The central MiNLO result is about 10.3% (9.6%) smaller than the MiNNLO PS (NNLO) prediction and features much larger scale scale uncertainties. The MiNNLO PS result instead agrees with NNLO at the sub-percent level, well within the perturbative uncertainties. Small numerical differences are expected even for inclusive observables, since the MiNNLO PS and NNLO calculations differ by terms beyond accuracy.
In Fig. 1 we examine a set of differential distributions. To validate MiNNLO PS , we compare it to the NNLO prediction without fiducial cuts, which could lead to significant differences due to the PS. Experimental data from the CMS collaboration unfolded and extrapolated to the inclusive phase space [99], and divided by the appropriate branching fractions, are also shown. The topleft plot shows the rapidity difference between the tt system and the leading jet defined with p T, j1 ≥ 120 GeV. Both MiNLO and MiNNLO PS are formally NLO accurate in this case, and the agreement between them indicates that this accuracy is retained by the MiNNLO PS procedure. The same conclusion holds for other observables that require at least one resolved hard jet.
The distributions in the average top-quark rapidity (y tav ) and transverse momentum (p T,tav ) as well as in the invariant mass (m tt ) and rapidity (y tt ) of the tt pair shown in Fig. 1 are inclusive over QCD radiation. For such distributions MiNNLO PS is expected to be NNLO accurate. Indeed, MiNNLO PS and NNLO yield consistent results, with fully overlapping uncertainty bands. The small differences in the central value are once again due to the different treatment of terms beyond NNLO accuracy. The larger uncertainty bands of the MiNNLO PS predictions are expected, since additional scale dependent terms are included within the first term in the r.h.s. of Eq. (11) that are not present in the fixed-order calculation. In comparison to the MiNLO results the inclusion of NNLO corrections through MiNNLO PS has an impact of about 10%-20% on the differential distributions and substantially reduces the perturbative uncertainties. Also the agreement with data is quite remarkable. All data points are within one standard deviation from the MiNNLO PS prediction, with the exception of the very first bin in the m tt distribution that, on the other hand, is strongly affected by the finite width of the top, whose effects are not included here.
We finally discuss the transverse-momentum spectrum of the tt pair, denoted by p T, tt in the bottom-right panel of Fig. 1. At large transverse momenta, the three predictions considered are effectively NLO accurate. Indeed, MiNLO and MiNNLO PS are essentially indistinguishable in that region, and at the same time consistent with the spectrum at fixed order. The small differences with NNLO are due to the generation of further radiation by the PS. At small transverse momenta, MiNNLO PS induces O(10%) corrections with respect to MiNLO and significantly reduces the large scale dependence. In this region, it also differs in shape from the NNLO calculation, which diverges and becomes unphysical for vanishing transverse momenta. Within the relatively large experimental errors, MiNNLO PS slightly improves the description of the data in terms of shape compared to NNLO for this observable.
In this letter we have presented the matching of the NNLO computation for top-quark pair production at FIG. 1. Distribution in the rapidity difference between the tt pair and the leading jet (∆y tt,j 1 ), in the rapidity (yt av ) and the average transverse-momentum (pT,t av ) of the top and the anti-top, as well as in the rapidity (y tt ), in the invariant mass (m tt ) and in the transverse momentum (p T,tt ) of the tt system. Predictions are shown for MiNNLOPS (blue, solid), MiNLO (black, dashed) and at NNLO (red, dashed). The black data points represent the CMS measurement at 13 TeV of Ref. [99], where the yt av and pT,t av distributions have been obtained with leptonically decaying top quarks.
hadron colliders with parton showers. This result has been obtained by constructing the MiNNLO PS method for the production of heavy quarks, which constitutes the first NNLO+PS prediction for reactions with colour charges in the final state in hadronic collisions. The comparisons presented in Fig. 1 provide a numerical validation of MiNNLO PS for top-quark pair production, demonstrating its NNLO accuracy. The simulations presented here also allow for the inclusion of the top-quark decay, paving the way to an accurate event generation for tt production at the LHC which will enable precise comparisons of fiducial measurements to theory.

SUPPLEMENTAL MATERIAL
In this appendix we complement the results presented in the letter by comparing MiNNLO PS , MiNLO , and NNLO predictions with different scale settings. In particular, we set µ (0) R = K R m tt and µ 0 = m tt e −L for MiNNLO PS and MiNLO , and we use µ 0 = m tt for the NNLO fixed-order calculation. The other settings are as reported in the main text. Table II reports the total cross sections, while Fig. 2 shows the same distributions as shown in the letter, but with the updated scale setting. Despite the fact that higher-order differences are expected for all observables, we observe an excellent agreement between MiNNLO PS and NNLO predictions for this scale choice.