Top-quark pair production at the LHC: Fully differential QCD predictions at NNLO

TL;DR

This paper presents a new fully differential NNLO QCD calculation for top-quark pair production at the LHC, enabling precise predictions of various observables and comparison with experimental data.

Contribution

It introduces a novel NNLO QCD computational framework using the $q_T$ subtraction method within the Matrix platform for top-quark pair production.

Findings

NNLO predictions match well with CMS data

Provides detailed differential distributions at 13 TeV

Enhances precision in top-quark physics analyses

Abstract

We report on a new fully differential calculation of the next-to-next-to-leading-order (NNLO) QCD radiative corrections to the production of top-quark pairs at hadron colliders. The calculation is performed by using the $q_T$ subtraction formalism to handle and cancel infrared singularities in real and virtual contributions. The computation is implemented in the Matrix framework, thereby allowing us to efficiently compute arbitrary infrared-safe observables for stable top quarks. We present NNLO predictions for several single- and double-differential kinematical distributions in $pp$ collisions at the centre-of-mass energy $\sqrt{s}=13$ TeV, and we compare them with recent LHC data by the CMS collaboration.

Figures

Equations

$$d{\sigma }_{\mathrm{NNLO}}^{t\bar{t}}={\mathcal{H}}_{\mathrm{NNLO}}^{t\bar{t}}\otimes d{\sigma }_{\mathrm{LO}}^{t\bar{t}}+\left[d{\sigma }_{\mathrm{NLO}}^{t\bar{t}+\mathrm{jet}}-d{\sigma }_{\mathrm{NNLO}}^{t\bar{t},CT}\right],$$

$${\sigma }_{\text{LO}}^{t\bar{t}}=478.9(1{)}_{-21.4\%}^{+29.6\%}\text{ pb},{\sigma }_{\text{NLO}}^{t\bar{t}}=726.9(1{)}_{-11.9\%}^{+11.7\%}\text{ pb},{\sigma }_{\text{NNLO}}^{t\bar{t}}=794.0(8{)}_{-5.7\%}^{+3.5\%}\text{ pb}.$$

$$H_T=m_{T,t}+m_{T,\bar{t}},$$

$$m_{T,t(\bar{t})}=\sqrt{m_t^2+p_{T,t(\bar{t})}^2},$$

Full text

ZU-TH 31/19

**Top-quark pair production at the LHC:

Fully differential QCD predictions at NNLO **

Stefano Catani(a), Simone Devoto*(b), Massimiliano Grazzini(b)*,

Stefan Kallweit*(c)*** and Javier Mazzitelli(b)**

*(a)*INFN, Sezione di Firenze and Dipartimento di Fisica e Astronomia,

Università di Firenze, I-50019 Sesto Fiorentino, Florence, Italy

*(b)*Physik Institut, Universität Zürich, CH-8057 Zürich, Switzerland

*(c)*Dipartimento di Fisica, Università degli Studi di Milano-Bicocca and

INFN, Sezione di Milano-Bicocca, I-20126, Milan, Italy

Abstract

We report on a new fully differential calculation of the next-to-next-to-leading-order (NNLO) QCD radiative corrections to the production of top-quark pairs at hadron colliders. The calculation is performed by using the $q_{T}$ subtraction formalism to handle and cancel infrared singularities in real and virtual contributions. The computation is implemented in the Matrix framework, thereby allowing us to efficiently compute arbitrary infrared-safe observables for stable top quarks. We present NNLO predictions for several single- and double-differential kinematical distributions in $pp$ collisions at the centre-of-mass energy $\sqrt{s}=13$ TeV, and we compare them with recent LHC data by the CMS collaboration.

June 2019

1 Introduction

The production of top quarks at high-energy colliders is a process of utmost importance, both in testing the validity of the Standard Model (SM) and in the quest for new physics. Within the SM, the main source of top-quark events in hadronic collisions is top-quark pair ($t{\bar{t}}$) production. The large data set delivered by the CERN LHC enables precise measurements of the $t{\bar{t}}$ production cross section as a function of the $t{\bar{t}}$ kinematics (see e.g. Refs. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]), which can be compared with the SM predictions. At the same time these studies have a wider relevance, since $t{\bar{t}}$ production is a crucial background in many new-physics searches.

Next-to-leading-order (NLO) QCD corrections for this production process were obtained thirty years ago [11, 12, 13, 14, 15]. Beyond the on-shell approximation of $t{\bar{t}}$ production, first NLO QCD studies were carried out within the narrow-width approximation [16, 17, 18]. Such NLO studies were later performed by considering the complete $W^{+}W^{-}b\bar{b}$ final states, with off-shell leptonic [19, 20, 21, 22] and semi-leptonic [23] decays. In the leptonic channel the case of massive bottom quarks was investigated in Refs. [24, 25]. NLO QCD results for off-shell $t{\bar{t}}$ production in association with an additional jet were obtained in Refs. [26, 27]. NLO electroweak (EW) corrections for on-shell $t{\bar{t}}$ production were studied in Refs. [28, 29, 30, 17, 31, 32], and a merged calculation for $t{\bar{t}}+0,1\,$jets including EW corrections was presented in Ref. [33]. For the leptonic decay channel the complete NLO EW corrections to the production of the six-particle final state are known [34].

The calculation of the next-to-next-to-leading-order (NNLO) QCD corrections to the $t{\bar{t}}$ total cross section was completed a few years ago [35, 36, 37, 38]. NNLO results for some differential distributions were presented in Refs. [39, 40, 41]. This calculation was recently combined with NLO EW corrections [42]. The $t{\bar{t}}$ charge asymmetry is known at NLO [43] and NNLO [44] in QCD, and also including NLO EW corrections [45]. First NNLO QCD results including top-quark decays are starting to appear [46].

In the present paper we deal with on-shell $t{\bar{t}}$ production in NNLO QCD. The calculation of the $t{\bar{t}}$ production cross section at this perturbative order requires tree-level contributions with two additional final-state partons, one-loop contributions with one additional parton and purely virtual contributions. The required tree-level and one-loop scattering amplitudes are known. They enter the NLO calculation of the associated production of a $t{\bar{t}}$ pair and one jet [47, 48], but in the case of NNLO $t{\bar{t}}$ production they need to be accurately evaluated also in the infrared-singular regions where the jet becomes unresolved. The purely virtual contributions entail the square of one-loop scattering amplitudes and the two-loop scattering amplitudes. The squared one-loop amplitudes are known [49, 50, 51]. The complete computation of the two-loop amplitudes has been carried out numerically [52, 53]. Partial results for these amplitudes are available in analytic form [54, 55, 56, 57]. Recent progress in the computation of non-planar two-loop master integrals [58, 59] indicates that the analytic calculation can be completed in the near future.

The implementation of the various scattering amplitudes in a (fully differential) NNLO calculation is definitely a non-trivial task because of the presence of infrared (IR) divergences at intermediate stages of the calculation. Various methods have been proposed and used to overcome these difficulties at the NNLO level (the interested reader can consult the list of references in Ref. [60]).

Using the antenna subtraction method [61, 62], partial results for $t{\bar{t}}$ production in the $q{\bar{q}}$ partonic channel were obtained by considering the complete fermionic contributions and evaluating the remaining contributions in the leading-colour approximation [63, 64, 65]. The complete NNLO computation of Refs. [35, 36, 37, 38, 44, 39, 40, 41] was performed by using the Stripper method [66, 67, 68].

In a recent paper [69] we have presented a new calculation of the inclusive $t{\bar{t}}$ production cross section in NNLO QCD. This calculation completes a previous computation that was limited to the flavour off-diagonal partonic channels [70]. The calculation uses the $q_{T}$ subtraction formalism [71] to handle and cancel IR-singular contributions in real and virtual corrections, and it is now completely integrated into the Matrix framework [72]. This allows us to perform fast and efficient computations of fiducial cross sections and (multi-)differential kinematical distributions for the production of on-shell top quarks.

In the present paper we extend the results of Ref. [69] in various respects. We present NNLO QCD predictions for several differential distributions in the transverse momenta and rapidities of the top quarks, as well as in the invariant mass and the rapidity of the $t{\bar{t}}$ system, and we discuss the results obtained by using different scale choices. We compare these results with CMS measurements in the lepton+jets channel at the centre-of-mass energy $\sqrt{s}=13$ TeV [9]. We then consider double-differential distributions and compare our results with the corresponding measurements by CMS [9].

The paper is organized as follows. In Section 2 we illustrate the framework in which the calculation is performed. In Section 3 we present results for single-differential and double-differential distributions, and we compare them with the experimental measurements. Finally, in Section 4 we present our conclusions. In Appendix A we present a quantitative comparison of our NNLO differential results with those available in the literature.

2 Calculation within the MATRIX framework

Our fully differential NNLO computation of $t{\bar{t}}$ production is carried out within the Matrix [72] framework. Matrix features a completely automated implementation of the $q_{T}$ subtraction formalism [71] to compute NNLO corrections, and it is thus applicable to the production of an arbitrary set of colourless final-state particles in hadronic collisions [73], as long as the two-loop virtual corrections to the corresponding leading-order (LO) process are provided. With appropriate modifications of the NNLO subtraction counterterm and the explicit computation of additional soft contributions (see below), Matrix can now also deal with the production of heavy-quark pairs.

According to the $q_{T}$ subtraction method, the NNLO differential cross section $d{\sigma}^{t{\bar{t}}}_{\rm NNLO}$ for the production process $pp\rightarrow t{\bar{t}}+X$ can be written as

[TABLE]

where $d{\sigma}^{t{\bar{t}}+\rm{jet}}_{\rm NLO}$ is the $t{\bar{t}}$+jet cross section at NLO accuracy.

The square bracket term of Eq. (1) is IR finite in the limit in which the transverse momentum of the $t{\bar{t}}$ pair, $q_{T}$, vanishes. However, the individual contributions $d{\sigma}^{t{\bar{t}}+\rm{jet}}_{\rm NLO}$ and $d{\sigma}^{t{\bar{t}},\,CT}_{\rm NNLO}$ are separately divergent. The contribution $d{\sigma}^{t{\bar{t}}+\rm{jet}}_{\rm NLO}$ can be evaluated with any available NLO method to handle and cancel IR divergences. The IR subtraction counterterm $d{\sigma}^{t{\bar{t}},\,CT}_{\rm NNLO}$ is obtained from the NNLO perturbative expansion (see e.g. Refs. [74, 75, 70]) of the resummation formula of the logarithmically-enhanced contributions to the $q_{T}$ distribution of the $t{\bar{t}}$ pair [76, 77, 78]: the explicit form of $d{\sigma}^{t{\bar{t}},\,CT}_{\rm NNLO}$ is fully known.

To complete the NNLO calculation, the second-order functions ${\cal H}^{t{\bar{t}}}_{\rm NNLO}$ in Eq. (1) are needed. These functions embody process-independent and process-dependent contributions. The process-independent contributions to ${\cal H}^{t{\bar{t}}}_{\rm NNLO}$ are analogous to those entering Higgs boson [71] and vector-boson [79] production, and they are explicitly known [80, 81, 73, 82, 83]. Since in $t{\bar{t}}$ production both the $gg$ and the $q{\bar{q}}$ partonic channels contribute at the same perturbative order, all these process-independent contributions are required. In the flavour off-diagonal channels the process-dependent contributions to ${\cal H}^{t{\bar{t}}}_{\rm NNLO}$ involve only amplitudes of the partonic processes $q{\bar{q}}\rightarrow t{\bar{t}}$ and $gg\rightarrow t{\bar{t}}$ up to the one-loop level, and the explicit results on the NLO azimuthal-correlation terms in the transverse-momentum resummation formalism [78]. The computation of ${\cal H}^{t{\bar{t}}}_{\rm NNLO}$ in the flavour diagonal $q{\bar{q}}$ and $gg$ channels additionally requires the two-loop amplitudes for $q{\bar{q}}\rightarrow t{\bar{t}}$ and $gg\rightarrow t{\bar{t}}$, and the evaluation of new contributions of purely soft origin. The two-loop amplitudes are available in a numerical form [53], and the corresponding grids have been implemented into Matrix through a suitable interpolation routine. The computation of the additional soft contributions has been completed by some of us [84]111An independent computation of these soft contributions is presented in Ref. [85]., and it has been implemented into Matrix as well.

The core of the Matrix framework is the Monte Carlo program Munich222Munich is the abbreviation of “MUlti-chaNnel Integrator at Swiss (CH) precision” — an automated parton-level NLO generator by S. Kallweit., which includes a fully automated implementation of the dipole-subtraction method for massless [86, 87] and massive [88] partons, and an efficient phase-space integration. All the required (spin- and colour-correlated) tree-level and one-loop (squared) amplitudes are obtained by using OpenLoops 2 [89, 90], except for the four-parton tree-level colour correlations that are based on an analytic implementation. OpenLoops 2 relies on its new on-the-fly tensor reduction [91] that guarantees stability all over the phase space, especially in the IR-singular regions, while scalar integrals from Collier [92, 93] are used. To the purpose of validating our results for the real–virtual corrections, we have also used the independent matrix-element generator Recola [94, 95], which employs tensor reduction and scalar integrals from Collier, and we find complete agreement.

The subtraction in the square brackets of Eq. (1) is not local, but the cross section is formally finite in the limit $q_{T}\rightarrow 0$. In practice, a technical cut on $q_{T}$ is introduced to render $d{\sigma}^{t{\bar{t}}+\mathrm{jet}}_{\mathrm{(N)LO}}$ and $d{\sigma}^{\mathrm{CT}}_{\mathrm{(N)NLO}}$ separately finite. Therefore, in our actual implementation, the $q_{T}$ subtraction method is very similar to a phase-space slicing method. It turns out that a cut, $r_{\mathrm{cut}}$, on the dimensionless quantity $r=q_{T}/m_{t{\bar{t}}}$ ($m_{t{\bar{t}}}$ denotes the invariant mass of the $t{\bar{t}}$ pair) is more convenient from a practical point of view. The absence of any residual logarithmic dependence on $r_{\mathrm{cut}}$ is a strong evidence of the correctness of the computation, since any mismatch between the contributions would result in a divergence of the cross section in the limit $r_{\mathrm{cut}}\rightarrow 0$. The remaining power-suppressed contributions vanish in that limit, and they can be controlled by monitoring the $r_{\mathrm{cut}}$ dependence of the cross section.

The $r_{\text{cut}}\rightarrow 0$ extrapolation for the total cross section is carried out by using the approach introduced in Ref. [72]. A quadratic least $\chi^{2}$ fit to the $r_{\text{cut}}$ dependent results is performed and repeated by varying the upper bound of the $r_{\text{cut}}$ interval. Finally, the result with the lowest $\chi^{2}/$degrees-of-freedom value is taken as the best fit, while the remaining results are used to estimate the extrapolation uncertainty. In addition to this analysis at the level of the total cross section, we have performed a similar bin-wise extrapolation in the computation of differential cross sections. We find that the results are in good agreement with those obtained by directly using a sufficiently low value of $r_{\text{cut}}$ ($r_{\rm cut}\lesssim 0.15\%$).

3 Results

To present our quantitative results, we consider $pp$ collisions at $\sqrt{s}=13$ TeV, and we fix the pole mass ${m_{t}}$ of the top quark to the value ${m_{t}}=173.3$ GeV. We consider $n_{F}=5$ massless quark flavours, and we use the corresponding NNPDF31 [96] sets of parton distribution functions (PDFs) with $\alpha_{\mathrm{S}}({m_{Z}})=0.118$. In particular, NnLO (with $n=0,1,2$) predictions are obtained by using PDFs at the corresponding perturbative order and the evolution of $\alpha_{\mathrm{S}}$ at $(n+1)$-loop order, as provided by the PDF set.

QCD scale uncertainties are estimated through the customary procedure of independently varying the renormalization ($\mu_{R}$) and factorization ($\mu_{F}$) scales by a factor of two around their common central value $\mu_{0}$ with the constraint $0.5\leq\mu_{F}/\mu_{R}\leq 2$, i.e. we use the standard 7-point scale variation.

Setting $\mu_{0}={m_{t}}$, the total cross sections and their corresponding scale uncertainties read

[TABLE]

We note that the LO and NLO results are not fully consistent within the corresponding uncertainties, indicating that, at least at LO, scale variations cannot be trusted as perturbative uncertainties. Similar features are shared by various other hard-scattering processes at hadron colliders. In contrast, the NLO and NNLO predictions are consistent, suggesting that scale variations can be used to estimate the size of perturbative contributions beyond NNLO.

The characteristic hard-scattering scale that controls the perturbative QCD behaviour of the total cross section $\sigma^{t{\bar{t}}}$ is ${m_{t}}$. In our calculation of $\sigma^{t{\bar{t}}}$, as reported in Eq. (2), we have used QCD scales ($\mu_{R}$ and $\mu_{F}$) at values of the order of ${m_{t}}$. Differential cross sections are controlled by corresponding characteristic hard scales, and we use QCD scales of that order in our computation of these observables. The characteristic hard scale specifically depends on the differential cross section under consideration.

Having at our disposal a fully differential calculation we can also use dynamical QCD scales. By dynamical we mean hard scales that refer to multi-differential cross sections eventually integrated over the phase space to obtain the specific differential cross section under consideration. The use of a dynamical scale produces practical simplifications since it allows us to compute several observables (e.g., differential cross sections) simultaneously, without changing the QCD scales on an observable-dependent basis. In practice, we use dynamical scales that are expected to be “effectively similar” to characteristic hard scales. Moreover, the study of dynamical scales is of interest independently of how we use them.

The default dynamical scale that we use throughout the paper is set to the central value $\mu_{0}=H_{T}/2$, where $H_{T}$ is the sum of the transverse masses of the top and antitop quarks,

[TABLE]

with

[TABLE]

and $p_{T,t}$ and $p_{T,{\bar{t}}}$ are the transverse momenta of the top and the antitop quark, respectively. We present differential cross sections that are obtained by using $\mu_{0}=H_{T}/2$ and values of $\mu_{0}$ of the order of the characteristic hard scale for that cross section. We also show results obtained by using central scales that are lowered by a factor of $1/2$. A reduced central scale, such as $H_{T}/4$, was considered in the studies of Ref. [97] on the basis of features of fastest perturbative convergence of some observables, and it was also already suggested in Ref. [21].

We have chosen the dynamical scale $H_{T}$ since it is expected to be parametrically of the same order as the characteristic hard scale of the observables that we examine in this paper. This a priori expectation is based on the kinematical features of these observables and on the general dynamical features of $t{\bar{t}}$ production. In the following paragraph we briefly comment about this. Independently of the expectation, throughout the paper we comment on the actual quantitative results that we obtain by using different QCD scales.

Owing to dynamics, the typical size of both $p_{T,t}$ and $p_{T,{\bar{t}}}$ is of the order of ${m_{t}}$ (see, e.g., Figs. 1–3 and 7). Therefore, in the case of observables that are inclusive over $p_{T,t}$ and $p_{T,{\bar{t}}}$, such as the total cross section and the pair rapidity distribution in Fig. 5, $H_{T}/2$ turns out to be of the same order as ${m_{t}}$, which is the characteristic hard scale for these observables. Analogously, since $p_{T,t}\sim p_{T,{\bar{t}}}$, $H_{T}/2$ turns out to be of the same order as the transverse masses, which are the characteristic hard scales for the differential cross sections in Figs. 1–3 and 7. The invariant mass $m_{t{\bar{t}}}$ of the $t{\bar{t}}$ pair is the characteristic hard scale in the case of the differential cross sections in Figs. 4, 6 and 8. The invariant mass is of the same order as $H_{T}$ with the exception of the kinematical subregions where the transverse momentum $p_{T,t{\bar{t}}}$ of the pair or the rapidity separation $|y_{t}-y_{\bar{t}}|$ between the top and the antitop quark are large. However, these subregions are dynamically suppressed, and therefore they give a minor contribution to the inclusive (over $p_{T,t{\bar{t}}}$ and $|y_{t}-y_{\bar{t}}|$) cross sections in Figs. 4, 6 and 8.

Our numerical results for differential cross sections are compared with the measurements of the CMS collaboration [9] (the data correspond to an integrated luminosity of $35.8~{}{\rm fb}^{-1}$) in the lepton+jets channel at parton level. The extrapolation from particle to parton level is carried out by the CMS collaboration in the inclusive phase space, and therefore no kinematical cuts are applied to obtain our theoretical predictions. To perform the comparison, our results are multiplied by the factor $0.292$, which corresponds to the value $0.438$ [98] of the semileptonic decay fraction of the $t{\bar{t}}$ pair, multiplied by a factor of $2/3$ since Ref. [9] considers only the decay into electrons and muons.

In Ref. [9] the CMS data for single- and double-differential distributions are compared to theoretical results obtained with the NLO Monte Carlo event generators POWHEG [99, 100, 101], interfaced either to PYTHIA8 [102] or to HERWIG++ [103], and MG5_aMC@NLO [104] interfaced to PYTHIA8 [102] (using the FxFx method [105] to deal with multijet merging). In addition, some of the measured parton-level single-differential distributions, namely the transverse-momentum and rapidity distributions of the leptonically and hadronically decaying top quark and the invariant-mass and rapidity distribution of the $t{\bar{t}}$ pair, are also compared to the NNLO QCD+NLO EW results of Ref. [42]. None of the double-differential distributions in Ref. [9] are compared to theoretical results beyond NLO QCD.

3.1 Single-differential distributions

In this section we present LO, NLO and NNLO results for a selection of single-differential distributions and compare them with the CMS measurements from Ref. [9]. At each perturbative order the scale-uncertainty bands in the figures are computed as explained at the beginning of Section 3.

We start the presentation by considering the transverse-momentum distributions of the top and antitop quarks. For each event we classify the transverse momenta according to their maximum and minimum values, $p_{T,t_{\text{high}}}$ and $p_{T,t_{\text{low}}}$.

In Figs. 1 and 2 (left) we show these distributions333NNLO results for these distributions have been recently presented in Ref. [106]. computed at our reference scale $\mu_{0}=H_{T}/2$. The characteristic hard scale of a transverse-momentum distribution is the transverse mass $m_{T}$. Therefore, in Fig. 1 (central and right) we also report the $p_{T,t_{\text{high}}}$ distribution for central scales $\mu_{0}=m_{T,t_{\text{high}}}$ and $\mu_{0}=m_{T,t_{\text{high}}}/2$, respectively, while in Fig. 2 (central and right) we consider $p_{T,t_{\text{low}}}$ for $\mu_{0}=m_{T,t_{\text{low}}}$ and $\mu_{0}=m_{T,t_{\text{low}}}/2$. The $p_{T,t_{\text{high}}}$ distribution is peaked at $p_{T,t_{\text{high}}}\sim 100$ GeV, while the $p_{T,t_{\text{low}}}$ distribution is peaked at $p_{T,t_{\text{low}}}\sim 60$ GeV.

We first discuss the $p_{T,t_{\text{high}}}$ distribution and focus on the $p_{T,t_{\text{high}}}\rightarrow 0$ region. If $p_{T,t_{\text{high}}}$ is small, both top quarks are forced to have small transverse momenta. As a consequence, this kinematical region corresponds to a small transverse momentum of the top-quark pair, $p_{T,t{\bar{t}}}$. The small-$p_{T,t{\bar{t}}}$ region exhibits Sudakov-type divergences [76, 77, 78, 107] at fixed order in perturbation theory, from the strong unbalance of real and virtual contributions due to soft-collinear emissions. In the computation of $p_{T,t_{\text{high}}}$, the unphysical fixed-order behaviour of $p_{T,t{\bar{t}}}$ is smeared due to the integration over $p_{T,t_{\text{low}}}$, and the Sudakov-type perturbative divergences disappear, by leaving (possibly large) residual effects. The amount of smearing is controlled by the shape of the $p_{T,t_{\text{high}}}$ distribution in the low-$p_{T}$ region at LO, which affects the unbalance between real and virtual contributions. The steeply rising LO distribution at low $p_{T}$ strongly suppresses real radiation, and the NLO radiative corrections to $p_{T,t_{\text{high}}}$ tend to be large and negative as $p_{T,t_{\text{high}}}$ decreases (a large and positive effect occurs at NNLO, and so forth). Accurate theoretical predictions of the detailed shape of the $p_{T,t_{\text{high}}}$ distribution at small $p_{T}$ require studies of all-order resummation effects of Sudakov type. However, in the case of large $p_{T}$ bins (as is the case in Fig. 1), reliable predictions can be obtained by considering perturbation theory at a sufficiently high order.

Comparing the results in Fig. 1 for the scales $\mu_{0}=H_{T}/2$ (left) and $\mu_{0}=m_{T,t_{\text{high}}}$ (central) we see that they are rather similar, and that the NNLO prediction agrees with the data. The scale $\mu_{0}=m_{T,t_{\text{high}}}/2$ also leads to good agreement with the data, but the corresponding NNLO uncertainty band is significantly narrower, especially in the intermediate region of transverse momenta, which is not observed for the corresponding band at NLO. This behaviour, namely the drastic shrinking of the scale-uncertainty band from NLO to NNLO, might indicate that for this choice scale variations cannot be trusted as an estimate of the perturbative uncertainties at NNLO. We also note that in the intermediate and large $p_{T}$ region the result obtained by using $\mu_{0}=m_{T,t_{\text{high}}}/2$ coincides with the upper bound of the corresponding NNLO band (i.e., the point $\mu_{0}=m_{T,t_{\text{high}}}/2$ corresponds to a local maximum of the NNLO cross section as a function of the scales).

We now discuss the $p_{T,t_{\text{low}}}$ distribution. In the region $p_{T,t_{\text{low}}}\rightarrow 0$, for LO kinematics both the top and the antitop quark are required to have small $p_{T}$. At NLO, real corrections open up a phase-space region where one top quark has a small $p_{T}$ and the other one has a relatively large $p_{T}$, thereby leading to large positive radiative contributions. The perturbative instability affecting the $p_{T,t_{\text{high}}}\rightarrow 0$ behaviour is now spread over the entire region of transverse momenta since, contrary to the low-$p_{T,t_{\text{high}}}$ region, small values of $p_{T,t_{\text{low}}}$ do not constrain $p_{T,t{\bar{t}}}$ to be small. The choices $\mu_{0}=H_{T}/2$ and $\mu_{0}=m_{T,t_{\text{low}}}$ lead to rather similar results. In both cases, at low and large $p_{T,t_{\text{low}}}$ NLO and NNLO bands overlap, whereas they do not in the intermediate region where the NLO band shrinks, showing that NLO perturbative uncertainties are underestimated. We note that the scale $\mu_{0}=m_{T,t_{\text{low}}}/2$ makes the perturbative convergence worse and the scale uncertainties larger at both NLO and NNLO.

We next consider the distribution in the transverse momentum of the hadronically decaying top or antitop quark, $p_{T,t_{\text{had}}}$. Since our calculation refers to stable top quarks, a prediction for the $p_{T,t_{\text{had}}}$ distribution can be obtained by computing the transverse-momentum spectra of the top and the antitop quark, and taking their average afterwards.444This is also the definition used in the data/theory comparison performed by the CMS collaboration [9]. Discussing our theoretical predictions, we refer to this as the $p_{T,t_{\mathrm{av}}}$ distribution. The corresponding LO, NLO and NNLO results are depicted in Fig. 3 for three different scale choices. We show the predictions for our default choice $\mu_{0}=H_{T}/2$ on the left. In the predictions for our natural scale choices, the top (antitop) $p_{T}$ distributions required to compute the average are evaluated at the corresponding transverse mass, $\mu_{0}=m_{T,t({\bar{t}})}$ (central) and $\mu_{0}=m_{T,t({\bar{t}})}/2$ (right). We denote these scale choices as $m_{T,t_{\rm av}}$ and $m_{T,t_{\rm av}}/2$, respectively. The $p_{T,t_{\text{had}}}$ distribution has a maximum at $p_{T,t_{\text{had}}}\sim 80~{}{\rm GeV}$. The LO and NLO scale-uncertainty bands do not overlap, except for $\mu_{0}=m_{T,t_{\rm av}}/2$. This is consistent with what happens for the corresponding total cross sections. The NLO and NNLO bands do overlap in the entire $p_{T,t_{\text{had}}}$ range, suggesting a good convergence of the perturbative expansion. In Fig. 3 we also observe that the scale choices $\mu_{0}=H_{T}/2$ (left) and $\mu_{0}=m_{T,t_{\rm av}}$ (central) give rather similar results. On the contrary, the choice $\mu_{0}=m_{T,t_{\rm av}}/2$ suggests a faster convergence of the perturbative expansion [97]. However, we also note that with this scale choice the NLO scale dependence is similar to what we obtain with $\mu_{0}=H_{T}/2$ and $\mu_{0}=m_{T,t_{\rm av}}$, whereas the NNLO scale dependence is significantly smaller than at NLO, thereby suggesting a possible underestimation of the perturbative uncertainty at NNLO. We note that $\mu_{0}=m_{T,t_{\rm av}}/2$ is also the scale used for the NNLO QCD+NLO EW prediction [42] to which the CMS data are compared in Ref. [9]. The data show that the measured $p_{T,t_{\text{had}}}$ distribution is slightly softer than the NNLO prediction. This is noticed also by the CMS collaboration [9] and in previous comparisons between NNLO results and LHC measurements [1, 2, 3, 4, 5, 6, 7, 8]. However Fig. 3 shows that the NNLO result and the data are consistent within the respective uncertainties. Our predictions for the $p_{T}$ spectrum of the leptonically decaying top quark are, of course, identical to those for $p_{T,t_{\text{had}}}$, and they are not shown here. The comparison with the data shows similar features.

We add a few comments on the perturbative behaviour of the $p_{T,t_{\text{high}}}$, $p_{T,t_{\text{low}}}$ and $p_{T,t_{\mathrm{av}}}$ distributions presented in Figs. 1, 2 and 3. The three distributions are identical at LO, but their behaviour beyond LO is clearly very different. As we can see from Fig. 3, the shape of the $p_{T,t_{\mathrm{av}}}$ distribution is almost unchanged with respect to the LO prediction. This feature is somehow expected, since the transverse-momentum spectrum of the top (antitop) quark at higher orders is affected by recoiling hard multijet radiation, which leads to a partly softer spectrum only at quite high $p_{T,t_{\mathrm{av}}}$ (beyond the $p_{T,t_{\text{had}}}$ range in Fig. 3), where hard multijet radiation is kinematically suppressed.

The physical shape of the $p_{T,t_{\text{high}}}$ and $p_{T,t_{\text{low}}}$ distributions is expected to be different from $p_{T,t_{\mathrm{av}}}$ at low and intermediate values of $p_{T}$. There we roughly have $p_{T,t_{\text{high}}}-p_{T,t_{\text{low}}}\sim p_{T,t{\bar{t}}}$. The $p_{T,t{\bar{t}}}$ distribution, which is confined to $p_{T,t{\bar{t}}}=0$ at LO, has an average value of about 50 GeV (which is already achieved at NLO), and it is localized in the small-$p_{T,t{\bar{t}}}$ region with a peak around 10 GeV [76, 107]. We thus expect that the physical shape of $p_{T,t_{\text{high}}}$ ($p_{T,t_{\text{low}}}$) is harder (softer) than its LO counterpart, with shape distortions of few tens of GeV as given by the size of $p_{T,t{\bar{t}}}$. Indeed, this is what we can observe from the comparison between the data and the LO prediction at small and intermediate values of $p_{T}$ in Figs. 1 and 2. This shape distortion has a smaller effect at high values of both $p_{T,t_{\text{high}}}$ and $p_{T,t_{\text{low}}}$.

In view of these physical expectations, it is not surprising that the shape of the $p_{T,t_{\text{high}}}$ and $p_{T,t_{\text{low}}}$ distributions is strongly affected by beyond-LO contributions. As discussed before, their fixed-order perturbative features are a smoothened version of the corresponding features of the $p_{T,t{\bar{t}}}$ distribution [107], the smoother behaviour being due to the smearing that is produced by the integration of $p_{T,t{\bar{t}}}$ over the respective unobserved $p_{T}$.

The invariant-mass distribution of the top-quark pair is reported in Fig. 4. The distribution is peaked at $m_{t{\bar{t}}}\sim 400~{}{\rm GeV}$. The characteristic hard-scattering scale for this distribution is of the order of $m_{t{\bar{t}}}$ itself. We use our default scale choice $\mu_{0}=H_{T}/2$ (left), and two other central values, namely $\mu_{0}=m_{t{\bar{t}}}/2$ (central) and $\mu_{0}=H_{T}/4$ (right).

We first comment on the convergence of the perturbative series for the three scale choices. In the cases $\mu_{0}=H_{T}/2$ and $\mu_{0}=m_{t{\bar{t}}}/2$ we see that LO and NLO bands do not overlap, analogously to what was previously observed for the $p_{T,t_{\text{had}}}$ distribution and for the total cross section in Eq. (2). In the case $\mu_{0}=H_{T}/2$, the NNLO corrections enhance the NLO result by about $10\%$ in the peak region, and their effect slightly increases with $m_{t{\bar{t}}}$ up to about $15\%$ in the highest-$m_{t{\bar{t}}}$ bin. Using $\mu_{0}=m_{t{\bar{t}}}/2$, the NNLO effect is similar in the peak region, but it increases to about $20\%$ at high $m_{t{\bar{t}}}$. The NLO and NNLO bands do overlap using both $\mu_{0}=H_{T}/2$ and $\mu_{0}=m_{t{\bar{t}}}/2$. As observed in Ref. [97], the choice $\mu_{0}=H_{T}/4$ leads to a faster convergence of the perturbative series. However, in the region where $m_{t{\bar{t}}}>360~{}{\rm GeV}$ we see that the size of the scale-variation band is very much reduced in going from the NLO to the NNLO result. This behaviour suggests that the central scale $\mu_{0}=H_{T}/4$ is (accidentally) quite close to a region of (local) minimal sensitivity [108] of the scale dependence of the NNLO result. In view of this feature, we think that the NNLO scale variation band with $\mu_{0}=H_{T}/4$ likely underestimates the perturbative uncertainty in this $m_{t{\bar{t}}}$ region and, especially, in the intermediate mass range 400~{}{\rm GeV}\raisebox{-1.72218pt}{\hbox to0.0pt{,\sim,\hss}}\raisebox{1.72218pt}{,<,}m_{t{\bar{t}}}\raisebox{-1.72218pt}{\hbox to0.0pt{,\sim,\hss}}\raisebox{1.72218pt}{,<,}1~{}{\rm TeV}.

We now comment on the comparison with the data. The first bin, $300\,{\rm GeV}<m_{t{\bar{t}}}<360\,{\rm GeV}$, deserves a separate discussion. The experimental result in this bin is significantly above the theoretical predictions, independently of the scale choice. This disagreement may have various origins. We first note that the result of our calculation for on-shell top quarks is non-vanishing only in the limited region where $2{m_{t}}=346.6~{}{\rm GeV}<m_{t{\bar{t}}}<360~{}{\rm GeV}$. In this region the NLO and NNLO effects are largest, regardless of the scale choice. Moreover, this region is particularly sensitive to the value of the top-quark mass. If the top-quark mass is smaller555Incidentally, we note that in Ref. [10] the CMS collaboration uses NLO QCD results obtained with NNPDF31 PDFs to extract the value ${m_{t}}=170.81\pm 0.68$ GeV from a fit of a measurement of the triple-differential cross section as a funtion of $m_{t{\bar{t}}}$, $y_{t{\bar{t}}}$ and the multiplicity of additional jets. by a couple of GeV, the predicted cross section in this bin becomes larger, without significantly affecting the NNLO result in each of the higher-$m_{t{\bar{t}}}$ bins. Another possible reason for the discrepancy may be the unfolding procedure [9] that is used to convert the data from particle to parton level. We expect such unfolding procedure to be particularly delicate in the threshold region: both the actual value of the top-quark mass used in the extrapolation and off-shell effects may have a significant impact. Further perturbative QCD effects beyond NNLO may also be relevant, but they are not expected to be large enough to explain the observed discrepancy with respect to the data. Overall, excluding the first bin, the NNLO result for the central scale $\mu_{0}=m_{t{\bar{t}}}/2$ leads to the best agreement with the data. The result for $\mu_{0}=H_{T}/2$ is also consistent with the data within uncertainties.

We finally consider the distribution in the rapidity of the top-quark pair, $y_{t{\bar{t}}}$. A natural scale choice for this distribution is the top-quark mass, as in the case of the computation of the total cross section. In Fig. 5 we report the rapidity distribution of the $t{\bar{t}}$ system for our default choice $\mu_{0}=H_{T}/2$ (left), for $\mu_{0}={m_{t}}$ (central) and for $\mu_{0}=H_{T}/4$ (right). We first observe that the choices $\mu_{0}=H_{T}/2$ and $\mu_{0}={m_{t}}$ lead to rather similar results. The LO and NLO bands marginally overlap, and the impact of NNLO corrections in the central rapidity region is about $10\%$, consistently with what we find for the total cross section. As previously observed for the $m_{t{\bar{t}}}$ distribution, the choice $\mu_{0}=H_{T}/4$ leads to a faster convergence of the perturbative expansion [97], and to a strong reduction of scale uncertainties at NNLO. The data nicely agree with the NNLO predictions for all the three scale choices.

We conclude this section with a few comments on the effects of EW corrections on the distributions considered so far. The effect of the complete EW corrections to top-quark pair production was studied in Ref. [32], and NNLO QCD+NLO EW predictions for some differential observables were presented in Ref. [42]. In the case of the $p_{T,t_{\text{had}}}$ distribution, the EW corrections are negative and lead to percent-level effects at low transverse momenta, which increase up to about $-5\%$ in the highest-$p_{T}$ bins in Fig. 3. Such effects are of the order of (or smaller than) the residual perturbative uncertainties of the NNLO result, and significantly smaller than the current experimental uncertainties. An analogous conclusion can be drawn for the invariant-mass distribution, which receives EW corrections of about $+2\%$ in the small-$m_{t{\bar{t}}}$ region and $-2\%$ in the highest-$m_{t{\bar{t}}}$ region in Fig. 4. The impact of EW corrections on the rapidity distributions is typically of ${\cal O}(1\%)$ or smaller. Therefore, we conclude that the impact of EW corrections is not expected to significantly affect our comparison with the data. Similar conclusions can be drawn for the double-differential distributions that we present in the following. More precise data, with higher reach in transverse momenta, will definitely call for the inclusion of EW corrections [42].

3.2 Double-differential distributions

In this section we present results for double-differential distributions and compare them with the CMS measurements from Ref. [9]. We consider three double-differential distributions, namely

• •

the $y_{t{\bar{t}}}$ distribution in $m_{t{\bar{t}}}$ intervals (see Fig. 6);

• •

the $p_{T,t_{\text{had}}}$ distribution in $y_{t_{\text{had}}}$ intervals (see Fig. 7);

• •

the $m_{t{\bar{t}}}$ distribution in $p_{T,t_{\text{had}}}$ intervals (see Fig. 8).

We present results for the default scale choice $\mu_{0}=H_{T}/2$ and at an additional scale of the order of the characteristic hard scale of the corresponding double-differential distribution.

The $y_{t{\bar{t}}}$ distribution in four $m_{t{\bar{t}}}$ intervals is presented in Fig. 6. The central scales are $\mu_{0}=H_{T}/2$ (upper) and $\mu_{0}=m_{t{\bar{t}}}/2$ (lower), and they lead to similar results, consistently with those for the single-differential cross sections in Figs. 4 and 5. The impact of the radiative corrections is relatively uniform in $y_{t{\bar{t}}}$ and in $m_{t{\bar{t}}}$. The scale $\mu_{0}=m_{t{\bar{t}}}/2$ leads to slightly larger radiative corrections in the highest-$m_{t{\bar{t}}}$ interval. In the first $m_{t{\bar{t}}}$ interval, the NNLO prediction slightly undershoots the data, consistently with what is observed in the low-$m_{t{\bar{t}}}$ region in Fig. 4. We note that, at variance with the first $m_{t{\bar{t}}}$ bin in Fig. 4, here the first $m_{t{\bar{t}}}$ interval extends up to $450$ GeV, thereby leading to a better agreement between the NNLO prediction and the data. At high values of $m_{t{\bar{t}}}$, the scale $\mu_{0}=m_{t{\bar{t}}}/2$ leads to a slightly better agreement with the data. In the highest-$y_{t{\bar{t}}}$ bin, for both central scales, the experimental result is below the NNLO prediction in all invariant-mass intervals except the first one.

In Fig. 7 we present the $p_{T,t_{\text{had}}}$ distribution in $y_{t_{\text{had}}}$ intervals, where $y_{t_{\text{had}}}$ is the rapidity of the hadronically decaying top or antitop quark with transverse momentum $p_{T,t_{\text{had}}}$. The QCD results use the scales $\mu_{0}=H_{T}/2$ (upper) and $\mu_{0}=m_{T,t_{\rm av}}$ (lower) as a natural scale, analogously to the $p_{T,t_{\text{had}}}$ distribution in Fig. 3. We observe that the two scale choices lead to very similar results. The impact of the radiative corrections is rather uniform in the considered $y_{t_{\text{had}}}$ intervals. Discussing the single-differential $p_{T,t_{\text{had}}}$ spectrum in Fig. 3, we observed that the NNLO prediction is slightly harder than the data. In Fig. 7 we observe that this feature holds in all the $y_{t_{\text{had}}}$ intervals. Nonetheless, the overall agreement between the NNLO prediction and the data is good.

We finally discuss the $m_{t{\bar{t}}}$ distribution in $p_{T,t_{\text{had}}}$ intervals, which is shown in Fig. 8. We use the central scales $\mu_{0}=H_{T}/2$ (upper) and $\mu_{0}=m_{t{\bar{t}}}/2$ (lower). In each $p_{T,t_{\text{had}}}$ interval we have $p_{T,{\rm min}}<p_{T,t_{\text{had}}}<p_{T,{\rm max}}$. The introduction of kinematical cuts on $p_{T,t_{\text{had}}}$ significantly affects the $m_{t{\bar{t}}}$ distribution. In the presence of a cut, $p_{T,{\rm max}}$, on the maximum $p_{T}$ of the hadronically decaying top quark, high values of $m_{t{\bar{t}}}$ can be reached by either increasing the transverse momentum of the top quark decaying to a leptonic final state, or increasing the rapidity separation $|y_{t}-y_{\bar{t}}|$. These kinematical regions are both dynamically suppressed for heavy-quark production and, therefore, $p_{T,{\rm max}}$ produces a faster suppression at high invariant masses. This effect can be observed by comparing the $m_{t{\bar{t}}}$ single-differential distribution in Fig. 4 with the present double-differential distribution.

The presence of a lower cut, $p_{T,{\rm min}}$, on the minimum $p_{T,t_{\text{had}}}$ has an even more drastic effect on the $m_{t{\bar{t}}}$ distribution, as it imposes the constraint $m_{t{\bar{t}}}>2\sqrt{{m^{2}_{t}}+p_{T,{\rm min}}^{2}}\equiv 2m_{T,{\rm min}}$ for LO kinematics. Below this unphysical threshold the LO result vanishes, and the NLO and NNLO results are effectively LO and NLO predictions, respectively. As a consequence, they suffer from larger theoretical uncertainties, which is reflected by the stronger scale dependence. Above this threshold, the LO distribution sharply increases up to a kinematical peak close to $2m_{T,{\rm min}}$. Owing to this LO behaviour, soft-collinear radiation produces shape instabilities [109] in this $m_{t{\bar{t}}}$ region at each subsequent perturbative order. The qualitative behaviour of these shape instabilities is completely analogous to that observed and discussed in Ref. [110] (see Figs. 10 and 20 and related comments therein) in the case of diphoton production in the presence of $p_{T}$ cuts. These perturbative instabilities are localized in a very narrow region around the LO threshold, and therefore their effect is smeared if a sufficiently large bin size is considered, as is the case for the differential distribution in Fig. 8.

The comparison with the data shows that in the first two $p_{T,t_{\text{had}}}$ intervals the NNLO prediction undershoots the data in the first $m_{t{\bar{t}}}$ bin. This discrepancy at low-$m_{t{\bar{t}}}$ is not resolved in the two highest-$p_{T,t_{\text{had}}}$ intervals, since the larger bin size ($300\text{ GeV}<m_{t{\bar{t}}}<430\text{ GeV}$) renders the distribution less sensitive to the behaviour close to the physical threshold. These observations are consistent with the expectations from the behaviour of the single-differential distribution in Fig. 4 at low $m_{t{\bar{t}}}$. Excluding the narrower bins at low $m_{t{\bar{t}}}$, the NNLO prediction in Fig. 8 is in very good agreement with the experimental measurements. The results with $\mu_{0}=H_{T}/2$ (upper) and $\mu_{0}=m_{t{\bar{t}}}/2$ (lower) turn out to be rather similar, consistently with our general expectation at the beginning of Section 3. The exception is the highest-$p_{T,t_{\text{had}}}$ interval, where the scale $\mu_{0}=m_{t{\bar{t}}}/2$ leads to a quite large NLO scale dependence at low $m_{t{\bar{t}}}$, which is drastically reduced at NNLO. In this region of low $m_{t{\bar{t}}}$ and high $p_{T,t_{\text{had}}}$, the invariant mass $m_{t{\bar{t}}}$ is not the characteristic hard scale anymore: the scale choice $\mu_{0}=m_{t{\bar{t}}}/2$ is not expected to be optimal, and the scale $\mu_{0}=H_{T}/2$ turns out to be more appropriate.

4 Summary and outlook

In this paper we have presented a new fully differential NNLO calculation of top-quark pair production at hadron colliders. The calculation is carried out by using the $q_{T}$ subtraction formalism to handle IR divergences from real and virtual contributions, and it is implemented in the Matrix framework. Our code enables fast and efficient calculations of fiducial cross sections and multi-differential distributions.

We have computed several single- and double-differential distributions of the top quarks, and we have compared our results with recent measurements performed by the CMS collaboration in the lepton+jets decay channel. We have considered several values of the renormalization and factorization scales to compute each of the distributions. We have used natural scales (i.e. ${m_{t}}$, $m_{t{\bar{t}}}/2$ and the relevant transverse masses $m_{T}$) of the order of the characteristic hard scale of the computed distribution, and we have shown that the corresponding results are similar to what is obtained with the overall choice $\mu_{0}=H_{T}/2$. We find that both the natural scales and $\mu_{0}=H_{T}/2$ lead to a reasonable perturbative behaviour for all the distributions that we have examined. The NNLO corrections substantially reduce the uncertainties of the theoretical predictions, and they improve the overall agreement with the experimental measurements. The largest deviation between data and the NNLO result occurs close to the $m_{t{\bar{t}}}$ threshold in single- and double-differential distributions. This discrepancy could be related to a variety of effects, including issues in the extrapolation of the data from particle to parton level, which is expected to be delicate in such threshold region. A lower value of the top-quark mass also has a significant impact close to the threshold.

The code that is used to perform these calculations is going to become public in a future Matrix release, providing a fast and flexible tool to compute (multi-)differential distributions with arbitrary cuts on the top-quark kinematical variables. The inclusion of NLO EW corrections and of top-quark decays is left to future work.

Acknowledgements

We are very grateful to Hayk Sargsyan for his contribution at early stages of this work. We are also indebted to Federico Buccioni, Jean-Nicolas Lang, Jonas Lindert and Stefano Pozzorini for their ongoing support with OpenLoops 2. We wish to thank Ben Kilminster and Florencia Canelli for useful discussions. This work is supported in part by the Swiss National Science Foundation (SNF) under contract 200020-169041. The work of SK is supported by the ERC Starting Grant 714788 REINVENT.

Appendix A Appendix

The calculation presented in this paper represents the first complete application of the $q_{T}$ subtraction formalism to the NNLO computation of (multi-)differential cross sections for the production of a colourful final state. As a validation of our results, in this appendix we carry out a detailed comparison of our predictions with those available for $t{\bar{t}}$ production at NNLO.

In Ref. [70], the $q_{T}$ subtraction results for the contributions of the flavour off-diagonal partonic channels to the $t{\bar{t}}$ total cross section were compared with those of the numerical program Top++ [111] by considering $p{\bar{p}}$ collisions at $\sqrt{s}=2$ TeV and $pp$ collisions at $\sqrt{s}=8$ TeV. In Ref. [69], the comparison with Top++ results at the LHC energies $\sqrt{s}=8$ and 13 TeV was performed by considering all the partonic channels, and we found agreement within the numerical uncertainties. We have successfully repeated this comparison in $p{\bar{p}}$ collisions at Tevatron energies and in $pp$ collisions with centre-of-mass energies $\sqrt{s}$ up to 100 TeV. We thus conclude that our calculation of the total cross section is in perfect agreement with the calculation of Refs. [35, 36, 37, 38]. The typical computing time needed to obtain ${\cal O}(0.1\%)$ precise inclusive cross sections (including scale uncertainties) with our program is about 1000 CPU days.

We point out that the only ingredient that is not computed independently in our implementation and in the calculation of Refs. [35, 36, 37, 38] is the finite part of the two-loop virtual amplitudes [52, 53]. This contribution, however, turns out to have a very small quantitative impact, namely ${\cal O}(0.1\%)$ of the NNLO total cross section at $\sqrt{s}=13$ TeV.

In the following we present a comparison of our NNLO differential results with the analogous results from Ref. [97]. To this purpose, we exactly follow the setup therein: we consider $pp$ collisions at $\sqrt{s}=13$ TeV, and we use CT14NNLO [112] PDFs. The QCD running of $\alpha_{\mathrm{S}}$ is evaluated at three-loop order with $\alpha_{\mathrm{S}}(m_{Z})=0.118$, and the pole mass of the top quark is fixed to ${m_{t}}=173.3$ GeV. We consider four differential distributions: invariant mass ($m_{t{\bar{t}}}$) and absolute rapidity ($|y_{t{\bar{t}}}|$) of the top-quark pair, and averages of the transverse momenta ($p_{T,t_{\mathrm{av}}}$) and of the absolute rapidities ($|y_{t_{\mathrm{av}}}|$) of the top and the antitop quark.

The values of the renormalization and factorization scales are fixed to different values for the different distributions. In case of the $m_{t{\bar{t}}}$, $|y_{t{\bar{t}}}|$ and $|y_{t_{\mathrm{av}}}|$ distributions the scales are set to $H_{T}/4$. In case of the $p_{T,t_{\mathrm{av}}}$ distribution, each of the $p_{T,t}$ and $p_{T,{\bar{t}}}$ distributions is calculated with both scales set to half of the corresponding transverse mass. This scale choice corresponds to $\mu_{0}=m_{T,t_{\rm av}}/2$ in the notation of Section 3.1.

The comparison of our results with those of Ref. [97] is illustrated in Fig. 9. In each case the upper panel shows the respective distribution, while the lower panel reports the ratio with respect to the corresponding result from Ref. [97] (CHM). Our results are stated with their numerical uncertainties. The computing time needed to obtain our results is about 2000 CPU days, mainly to achieve good statistical precision in the tails of the kinematical distributions. The results from Ref. [97] are quoted without an associated statistical uncertainty. However, Ref. [97] states that “the narrowest bins possible” are chosen in order to “keep the Monte Carlo integration error within about 1% in almost all bins”. Correspondingly, the approximate uncertainty of $\pm 1\%$ is reported in the plots for reference. We see that the agreement for all the considered distributions is excellent, also in the extreme kinematical regions, i.e. large rapidities, large transverse momenta and large invariant mass of the $t{\bar{t}}$ pair.